Create Voronoi regions using Python
Overview
One of the most common spatial problems is to find the nearest point of interest (POI) from our current location. Let’s say someone will run out of gas soon and he/she needs to find the closest fuel station before it’s too late, what is the best solution for this problem? For sure the driver can check the map to find the nearest fuel station but it can be problematic if there are more than one station in the area and he/she needs to quickly decide which station is the closest. The best solution is by representing each POI with a dot inside a polygon shape. So within a polygon, the closest POI is definitely the dot inside the polygon. These polygons are called Voronoi regions.
Data Collection
For this project I create Voronoi regions on the map based on POI data. All POI data are chosen randomly while the street network data downloaded from OpenStreetMap with the help of OSMnx package.
Create Voronoi Regions
Currently the easiest way to build Voronoi regions using Python is by using geovoronoi package. Geovoronoi is a package to create and plot Voronoi regions inside geographic areas. As for the map visualization, I choose folium package.
First I start by creating random points around the map.
gdf = gpd.GeoDataFrame()
gdf = gdf.append({'geometry': Point(106.644085,-6.305286)}, ignore_index=True)
gdf = gdf.append({'geometry': Point(106.653261,-6.301309)}, ignore_index=True)
gdf = gdf.append({'geometry': Point(106.637751,-6.284774)}, ignore_index=True)
gdf = gdf.append({'geometry': Point(106.665062,-6.284598)}, ignore_index=True)
gdf = gdf.append({'geometry': Point(106.627582,-6.283521)}, ignore_index=True)
gdf = gdf.append({'geometry': Point(106.641365,-6.276593)}, ignore_index=True)
gdf = gdf.append({'geometry': Point(106.625972,-6.303643)}, ignore_index=True)Next step is to determine the coverage area of Voronoi regions & save it to the geodataframe.
area_max_lon = 106.670929
area_min_lon = 106.619602
area_max_lat = -6.275227
area_min_lat = -6.309795
lat_point_list = [area_min_lat, area_max_lat,area_max_lat,area_min_lat]
lon_point_list = [area_min_lon, area_min_lon, area_max_lon, area_max_lon]
polygon_geom = Polygon(zip(lon_point_list, lat_point_list))boundary = gpd.GeoDataFrame()
boundary = boundary.append({'geometry': polygon_geom}, ignore_index=True)
Don’t forget to convert both gdf and boundary dataframes to Web Mercator projection.
gdf.crs = {'init' :'epsg:3395'}
boundary.crs = {'init' :'epsg:3395'}Convert the boundary geometry dataframe into an union of the polygon and POI dataframe into an array of coordinates.
boundary_shape = cascaded_union(boundary.geometry)
coords = points_to_coords(gdf.geometry)Calculate Voronoi regions.
poly_shapes, pts, poly_to_pt_assignments = voronoi_regions_from_coords(coords, boundary_shape)Create a graph from OpenStreetMap within the boundaries of the coverage area. Using the graph to collect all street networks within the boundaries of the coverage area and save it to dataframe.
G = ox.graph_from_polygon(boundary_shape, network_type='all_private')
gdf_all_streets = ox.graph_to_gdfs(G, nodes=False, edges=True,node_geometry=False, fill_edge_geometry=True)Create new dataframe to collect street networks within each Voronoi regions
gdf_streets_by_region = gpd.GeoDataFrame()
for x in range(len(poly_shapes)):
gdf_streets = gpd.GeoDataFrame()
gdf_streets['geometry'] = gdf_all_streets.intersection(poly_shapes[x])
gdf_streets['voronoi_region'] = x
gdf_streets = gdf_streets[gdf_streets['geometry'].astype(str) != 'LINESTRING EMPTY']
gdf_streets_by_region = gdf_streets_by_region.append(gdf_streets)Below is the visualization of the Voronoi regions on the map.
Conclusion
The map looks great! Unfortunately it is not yet ready to be used on real life application, the problem is these Voronoi regions created by using euclidean distance instead of network distance.
Hope you enjoyed reading this article. For more details of the code, you can just click this link.
