In [1]:
# GNU LESSER GENERAL PUBLIC LICENSE
# Version 3, 29 June 2007.
# Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>.
# Everyone is permitted to copy and distribute verbatim copies of this
# license document, but changing it is not allowed.
In [2]:
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
from scipy.spatial import Voronoi, voronoi_plot_2d
%matplotlib inline
mpl.rcParams['figure.figsize'] = 20,20 #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 [3]:
t1_shp_reads = time.time()
'''
# Waverly Hills streets and bounding box
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
'''
# Leon County streets and bounding box
streets = gpd.read_file(path+'Leon_County/LCSTSEG.shp')
streets.to_crs(epsg=2779, inplace=True) # NAD83(HARN) / Florida North
streets.to_file(path+'leon/leon.shp')
shp_Leon = ps.open(path+'leon/leon.shp')
b_box = shp_Leon.bbox
t2_shp_reads = time.time()-t1_shp_reads
print t2_shp_reads
In [4]:
# Individual street buffers in meters
intial_buffer = streets.buffer(700)
intial_buffer[:5]
Out[4]:
In [5]:
# Union of individual buffers
union_buffer = intial_buffer.unary_union
union_buffer = gpd.GeoSeries(union_buffer)
union_buffer.crs = streets.crs
union_buffer = gpd.GeoDataFrame(union_buffer, crs=streets.crs)
union_buffer.columns = ['geometry']
In [6]:
# Create convex hull buffer union
convex_hull_of_union = gpd.GeoDataFrame(union_buffer.convex_hull)
convex_hull_of_union.columns = ['geometry']
In [7]:
# Generate (n) simulated households
n = 100000
t1_allpoints = time.time()
x = np.random.uniform(b_box[0], b_box[2], n)
np.random.seed(850)
y = np.random.uniform(b_box[1], b_box[3], n)
simulated_coords = zip(x,y)
simulated_coord_points = [shapely.geometry.Point(xy_coord) for xy_coord in simulated_coords]
simulated_households = gpd.GeoDataFrame(simulated_coord_points, crs=streets.crs)
simulated_households.columns = ['geometry']
t2_allpoints = time.time()-t1_allpoints
print t2_allpoints, 'seconds to simulate households'
In [8]:
# Identify households within the buffer
t1_intersection = time.time()
intersection = [union_buffer['geometry'].intersection(hh) for hh in simulated_households['geometry']]
t2_intersection = time.time()-t1_intersection
print t2_intersection, 'seconds to identify households within the street buffer'
intersection = gpd.GeoDataFrame(intersection, crs=streets.crs)
intersection.columns = ['geometry']
# Ignore households not in the buffer --> shapely.geometry.collection.GeometryCollection
households_gdf = intersection[intersection.geom_type != 'GeometryCollection']
households_gdf[:5]
Out[8]:
In [9]:
# Create list of household coordinates to feed into Voronoi calculator
household_coords = []
for ids,coords in households_gdf['geometry'].iteritems():
household_coords.append([coords.x, coords.y])
In [26]:
household_coords = pd.DataFrame(household_coords)
household_coords.to_csv(path+'household_coords.csv')
In [27]:
# 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 [12]:
# Create geometries from the vertices and ridges of the Voronoi
'''
# Option 1: ignore all lines that go to infinity (-1)
polylines = [shapely.geometry.LineString(vor.vertices[line])
for line in vor.ridge_vertices
if -1 not in line]
# Option 2: use all lines
polylines = [shapely.geometry.LineString(vor.vertices[line])
for line in vor.ridge_vertices]
'''
# Option 3: remove -1 from ridges to make all voronoi polygons finite
vor_polys = []
for region in vor.regions:
if region:
if -1 in region:
region = region[:] # so we don't change the original voronoi output
region.remove(-1) # there can only be one -1 in the index list
vor_polys.append(shapely.geometry.MultiPoint(vor.vertices[region]).convex_hull)
In [13]:
# Create a dataframe of the complete set of Voronoi polygons
voronoi_all_df = gpd.GeoDataFrame(vor_polys, columns=['geometry'], crs=streets.crs)
In [14]:
# Create a dataframe of the complete set of Voronoi polygon centroids
voronoi_all_centroids = gpd.GeoDataFrame(voronoi_all_df.centroid, columns=['geometry'], crs=streets.crs)
In [15]:
# Create a dataframe of the set of Voronoi polygons within the convex hull of the street buffer union
poly_in_hull = [convex_hull_of_union['geometry'].intersection(hh_vor) for hh_vor in voronoi_all_df['geometry']]
voronoi_in_df = gpd.GeoDataFrame(poly_in_hull, crs=streets.crs)
voronoi_in_df.columns = ['geometry']
In [16]:
# Create a dataframe of the set of Voronoi polygon centroids within the bounding box
voronoi_in_centroids = gpd.GeoDataFrame(voronoi_in_df.centroid, crs=streets.crs)
voronoi_in_centroids.columns = ['geometry']
In [17]:
# Write out all dataframes to .shp
#households_gdf.to_file(path+'households/households.shp')
#convex_hull_of_union.to_file(path+'convex_hull/hull.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 [18]:
print len(vor_polys)
print len(households_gdf)
print len(voronoi_in_df)
print len(voronoi_all_df)
print 'The number of household points is equal to the number of voronoi polygons? --> ',\
len(vor_polys) == len(households_gdf) == len(voronoi_in_df) == len(voronoi_all_df)
In [19]:
# All finite voronoi polygons
voronoi_all_df.plot()
Out[19]:
In [20]:
# All households and all finite voronoi polygons trimmed down to a convex hull of the union buffer
voronoi_in_df.plot()
households_gdf.plot()
Out[20]:
In [21]:
convex_hull_of_union.plot()
union_buffer.plot()
Out[21]: