In [1]:
import numpy as np
import pysal as ps
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import matplotlib as mpl
import shapely
import time
import sys
from scipy.spatial import Voronoi, voronoi_plot_2d
%matplotlib inline
mpl.rcParams['figure.figsize'] = 25,25 #set the default map size
mpl.rcParams['patch.linewidth'] = 0.5 #set default polygon line width
np.random.seed(352)
# Local path on user's machine
path = '/Users/jgaboardi/Algorithms_etc./Data/'
In [2]:
# Waverly Hills
streets = gpd.read_file(path+'Waverly_Trim/Waverly.shp')
streets.to_crs(epsg=2779, inplace=True) # NAD83(HARN) / Florida North
streets.to_file(path+'waverly/waverly.shp')
shp_waverly = ps.open(path+'waverly/waverly.shp')
b_box = shp_waverly.bbox
In [3]:
intial_buffer = streets.buffer(50) #Buffer in meters
intial_buffer[:5]
Out[3]:
In [4]:
union_buffer = intial_buffer.unary_union #Buffer Union
union_buffer = gpd.GeoSeries(union_buffer)
union_buffer.crs = streets.crs
final_buffer = gpd.GeoDataFrame(union_buffer, crs=streets.crs)
final_buffer.columns = ['geometry']
In [5]:
t1_allpoints = time.time()
np.random.seed(352)
x = np.random.uniform(b_box[0], b_box[2], 100)
np.random.seed(850)
y = np.random.uniform(b_box[1], b_box[3], 100)
coords = zip(x,y)
coords = [shapely.geometry.Point(i) for i in coords]
simulated_hh = gpd.GeoDataFrame(coords)
simulated_hh.crs = streets.crs
simulated_hh.columns = ['geometry']
t2_allpoints = time.time()-t1_allpoints
print t2_allpoints, 'seconds to simulate households in Waverly'
In [6]:
intersection = [final_buffer['geometry'].intersection(hh) for hh in simulated_hh['geometry']]
intersection = gpd.GeoDataFrame(intersection, crs=streets.crs)
intersection.columns = ['geometry']
intersection[:5]
Out[6]:
In [7]:
# Remove shapely.geometry.collection.GeometryCollection
t1_allpoints = time.time()
hh_in = []
for hh in intersection['geometry']:
if type(hh) == shapely.geometry.point.Point:
hh_in.append(hh)
In [8]:
households_gdf = gpd.GeoDataFrame(hh_in, crs=streets.crs)
households_gdf.columns = ['geometry']
households_gdf[:5]
Out[8]:
In [9]:
household_coords = []
for i,j in households_gdf['geometry'].iteritems():
household_coords.append([j.x, j.y])
In [10]:
def in_box(households, bounding_box):
return np.logical_and(np.logical_and(bounding_box[0] >= households[:, 0],
households[:, 0] >= bounding_box[1]),
np.logical_and(bounding_box[2] >= households[:, 1],
households[:, 1] >= bounding_box[3]))
In [11]:
def voronoi(households, bounding_box):
# Select households inside the bounding box
i = True
print i
# Mirror points
points_center = households[:]
print points_center
points_left = np.copy(points_center)
points_left[:, 0] = bounding_box[0] - (points_left[:, 0] - bounding_box[0])
points_right = np.copy(points_center)
points_right[:, 0] = bounding_box[1] + (bounding_box[1] - points_right[:, 0])
points_down = np.copy(points_center)
points_down[:, 1] = bounding_box[2] - (points_down[:, 1] - bounding_box[2])
points_up = np.copy(points_center)
points_up[:, 1] = bounding_box[3] + (bounding_box[3] - points_up[:, 1])
points = np.append(points_center,
np.append(np.append(points_left,
points_right,
axis=0),
np.append(points_down,
points_up,
axis=0),
axis=0),
axis=0)
# Compute Voronoi
vor = Voronoi(points)
# Filter regions
regions = []
for region in vor.regions:
flag = True
for index in region:
if index == -1:
flag = False
break
else:
x = vor.vertices[index, 0]
y = vor.vertices[index, 1]
if not(bounding_box[0] - eps <= x and x <= bounding_box[1] + eps and
bounding_box[2] - eps <= y and y <= bounding_box[3] + eps):
flag = False
break
if region != [] and flag:
regions.append(region)
vor.filtered_points = points_center
vor.filtered_regions = regions
return vor
In [12]:
def centroid_region(vertices):
# Polygon's signed area
A = 0
# Centroid's x
C_x = 0
# Centroid's y
C_y = 0
for i in range(0, len(vertices) - 1):
s = (vertices[i, 0] * vertices[i + 1, 1] - vertices[i + 1, 0] * vertices[i, 1])
A = A + s
C_x = C_x + (vertices[i, 0] + vertices[i + 1, 0]) * s
C_y = C_y + (vertices[i, 1] + vertices[i + 1, 1]) * s
A = 0.5 * A
C_x = (1.0 / (6.0 * A)) * C_x
C_y = (1.0 / (6.0 * A)) * C_y
return np.array([[C_x, C_y]])
In [13]:
#b_box_flip = [162927.08355677937, 621380.1358320742, 166307.62055677705, 623680.8311320741]
In [14]:
eps = sys.float_info.epsilon
households = np.array(household_coords)
bounding_box = b_box # [x_min, x_max, y_min, y_max]
vor = voronoi(households, bounding_box)
In [15]:
fig = plt.figure()
ax = fig.gca()
# Plot initial points
ax.plot(vor.filtered_points[:, 0], vor.filtered_points[:, 1], 'b.')
# Plot ridges points
for region in vor.filtered_regions:
vertices = vor.vertices[region, :]
ax.plot(vertices[:, 0], vertices[:, 1], 'go')
# Plot ridges
for region in vor.filtered_regions:
vertices = vor.vertices[region + [region[0]], :]
ax.plot(vertices[:, 0], vertices[:, 1], 'k-')
# Compute and plot centroids
centroids = []
for region in vor.filtered_regions:
vertices = vor.vertices[region + [region[0]], :]
centroid = centroid_region(vertices)
centroids.append(list(centroid[0, :]))
ax.plot(centroid[:, 0], centroid[:, 1], 'r.')
#ax.set_xlim([-0.1, 10.1])
#ax.set_ylim([-0.1, 10.1])
#plt.savefig("bounded_voronoi.png")
voronoi_plot_2d(vor)
#plt.savefig("voronoi.png")
Out[15]:
In [16]:
# Instantiate and Plot the Voronoi
t1_voronoi = time.time()
vor = Voronoi(household_coords)
voronoi_plot_2d(vor)
plt.show()
t2_voronoi = time.time()-t1_voronoi
print t2_voronoi, 'seconds to instantiate the Voronoi Tessellation and plot the diagram'
print vor.npoints, 'points'
print len(vor.regions), 'regions'
In [ ]:
#polylines = [shapely.geometry.LineString(vor.vertices[line])
# for line in vor.ridge_vertices
# if -1 not in line]
polylines = [shapely.geometry.LineString(vor.vertices[line])
for line in vor.ridge_vertices]
In [ ]:
voronoi_all_df = gpd.GeoSeries([p for p in shapely.ops.polygonize(polylines)])
voronoi_all_df = gpd.GeoDataFrame(voronoi_all_df, crs=streets.crs)
voronoi_all_df.columns = ['geometry']
In [ ]:
voronoi_all_centroids = gpd.GeoDataFrame(voronoi_all_df.centroid, crs=streets.crs)
voronoi_all_centroids.columns = ['geometry']
In [ ]:
bounding_box = gpd.read_file(path+'bounding_box/bounding_box.shp')
In [ ]:
poly_in_bbox = [bounding_box['geometry'].intersection(hh) for hh in voronoi_all_df['geometry']]
voronoi_in_df = gpd.GeoDataFrame(poly_in_bbox, crs=streets.crs)
voronoi_in_df.columns = ['geometry']
In [ ]:
voronoi_in_centroids = gpd.GeoDataFrame(voronoi_in_df.centroid, crs=streets.crs)
voronoi_in_centroids.columns = ['geometry']
In [ ]:
households_gdf.to_file(path+'households/households.shp')
voronoi_all_centroids.to_file(path+'voronoi/centroids_all.shp')
voronoi_all_df.to_file(path+'voronoi/voronoi_all.shp')
voronoi_in_centroids.to_file(path+'voronoi/centroids_in.shp')
voronoi_in_df.to_file(path+'voronoi/voronoi_in.shp')
In [ ]:
In [ ]: