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.

Creating Voronoi Tesselations

Automated Simulation in GIS


James Gaboardi, 2016


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/'


/Users/jgaboardi/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.
  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')

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


31.3784589767

In [4]:
# Individual street buffers in meters
intial_buffer = streets.buffer(700)  
intial_buffer[:5]


Out[4]:
0    POLYGON ((625158.2342299456 162820.9589747752,...
1    POLYGON ((632572.7322698359 164252.7532209178,...
2    POLYGON ((622877.6070010761 158529.3292612745,...
3    POLYGON ((622819.3539096914 159626.6340049335,...
4    POLYGON ((625503.3645330973 164138.5211709534,...
dtype: object

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'


1.98290491104 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]


2775.11395979 seconds to identify households within the street buffer
Out[8]:
geometry
0 POINT (641008.2878558711 186443.6193693635)
6 POINT (588503.8617661271 143879.9120721709)
9 POINT (608431.7504617096 146625.5249416636)
10 POINT (632118.1956706066 173602.2221251963)
11 POINT (617652.3112889249 147384.488869649)

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'


14.1305770874 seconds to instantiate the Voronoi Tessellation and plot the diagram
45617 points
45618 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)


45617
45617
45617
45617
The number of household points is equal to the number of voronoi polygons? -->  True

In [19]:
# All finite voronoi polygons
voronoi_all_df.plot()


Out[19]:
<matplotlib.axes._subplots.AxesSubplot at 0x139dc8a10>

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]:
<matplotlib.axes._subplots.AxesSubplot at 0x1dfe0cad0>

In [21]:
convex_hull_of_union.plot()
union_buffer.plot()


Out[21]:
<matplotlib.axes._subplots.AxesSubplot at 0x1ef5ca7d0>