For this project I create a voronoi diagram on the map based on data points (or point of interests), voronoi diagram have applications in almost all areas of science and engineering. For geospatial use case, it is useful to tell us the closest point of interest (POI) 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.
Ok, let's start the project. As usual, for the first step let's import all necessary packages
import folium
import numpy as np
import geopandas as gpd
import osmnx as ox
import matplotlib.pyplot as plt
from shapely.ops import cascaded_union
from geovoronoi.plotting import subplot_for_map, plot_voronoi_polys_with_points_in_area
from geovoronoi import voronoi_regions_from_coords, points_to_coords
from shapely.geometry import Point, LineString, Polygon
Insert data points/POI into geodaframe
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)
gdf
geometry | |
---|---|
0 | POINT (106.64409 -6.30529) |
1 | POINT (106.65326 -6.30131) |
2 | POINT (106.63775 -6.28477) |
3 | POINT (106.66506 -6.28460) |
4 | POINT (106.62758 -6.28352) |
5 | POINT (106.64136 -6.27659) |
6 | POINT (106.62597 -6.30364) |
Determine the coverage area of the voronoi diagram & save it to geodaframe
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)
boundary
geometry | |
---|---|
0 | POLYGON ((106.61960 -6.30980, 106.61960 -6.275... |
Convert both dataframes to Web Mercator projection
gdf.crs = {'init' :'epsg:3395'}
boundary.crs = {'init' :'epsg:3395'}
/opt/anaconda3/lib/python3.7/site-packages/pyproj/crs/crs.py:53: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 return _prepare_from_string(" ".join(pjargs))
Convert the boundary geometry into a union of the polygon and POI data to NumPy array of coordinates.
boundary_shape = cascaded_union(boundary.geometry)
coords = points_to_coords(gdf.geometry)
boundary_shape
Calculate voronoi regions
poly_shapes, pts, poly_to_pt_assignments = voronoi_regions_from_coords(coords, boundary_shape)
fig, ax = subplot_for_map()
plot_voronoi_polys_with_points_in_area(ax, boundary_shape, poly_shapes, pts, poly_to_pt_assignments)
ax.set_title('Voronoi regions')
plt.tight_layout()
plt.show()
Create a graph from OSM within the boundaries of coverage area. Save the graph to geodataframe
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)
gdf_all_streets.head()
osmid | highway | oneway | length | geometry | service | access | name | lanes | est_width | bridge | junction | tunnel | maxspeed | width | u | v | key | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 816728678 | service | False | 76.900 | LINESTRING (106.63699 -6.28336, 106.63630 -6.2... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | 7627964418 | 7847970909 | 0 |
1 | 841101956 | service | True | 42.409 | LINESTRING (106.63699 -6.28336, 106.63684 -6.2... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | 7627964418 | 7847940435 | 0 |
2 | 841101955 | service | False | 70.538 | LINESTRING (106.63485 -6.28342, 106.63491 -6.2... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | 7627964419 | 7847970904 | 0 |
3 | [591697100, 695004245] | service | False | 305.493 | LINESTRING (106.63485 -6.28342, 106.63484 -6.2... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | 7627964419 | 6528394635 | 0 |
4 | 816728678 | service | False | 39.691 | LINESTRING (106.63485 -6.28342, 106.63502 -6.2... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | 7627964419 | 7627964413 | 0 |
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)
gdf_streets_by_region.head()
geometry | voronoi_region | |
---|---|---|
34 | LINESTRING (106.65206 -6.28108, 106.65205 -6.2... | 0 |
1126 | LINESTRING (106.64790 -6.28385, 106.64835 -6.2... | 0 |
1127 | LINESTRING (106.64790 -6.28385, 106.64761 -6.2... | 0 |
1128 | LINESTRING (106.64790 -6.28385, 106.64811 -6.2... | 0 |
1132 | LINESTRING (106.64854 -6.28395, 106.64838 -6.2... | 0 |
Visualize the voronoi diagram using folium
m = folium.Map([-6.304029, 106.645874], zoom_start=13, tiles='cartodbpositron')
#draw the voronoi diagram within coverage area
for x in range(len(poly_shapes)):
folium.GeoJson(poly_shapes[x]).add_to(m)
#draw the data points
points = np.array([[-6.305286, 106.644085], [-6.301309, 106.653261], [-6.284774, 106.637751], [-6.284598, 106.665062], [-6.283521, 106.627582],[-6.276593, 106.641365], [-6.303643, 106.625972]])
locs = points
for location in locs:
folium.CircleMarker(location=location,
color = "#4925a2", radius=0.01).add_to(m)
#draw street networks for each voronoi region
folium.GeoJson(gdf_streets_by_region[gdf_streets_by_region['voronoi_region'] == 0], style_function=lambda x: {'color': '#30e861', 'weight':'1'}).add_to(m)
folium.GeoJson(gdf_streets_by_region[gdf_streets_by_region['voronoi_region'] == 1], style_function=lambda x: {'color': '#ca1890', 'weight':'1'}).add_to(m)
folium.GeoJson(gdf_streets_by_region[gdf_streets_by_region['voronoi_region'] == 2], style_function=lambda x: {'color': '#d5e849', 'weight':'1'}).add_to(m)
folium.GeoJson(gdf_streets_by_region[gdf_streets_by_region['voronoi_region'] == 3], style_function=lambda x: {'color': '#1b12e4', 'weight':'1'}).add_to(m)
folium.GeoJson(gdf_streets_by_region[gdf_streets_by_region['voronoi_region'] == 4], style_function=lambda x: {'color': '#ee3d07', 'weight':'1'}).add_to(m)
folium.GeoJson(gdf_streets_by_region[gdf_streets_by_region['voronoi_region'] == 5], style_function=lambda x: {'color': '#eead07', 'weight':'1'}).add_to(m)
folium.GeoJson(gdf_streets_by_region[gdf_streets_by_region['voronoi_region'] == 6], style_function=lambda x: {'color': '#059b3a', 'weight':'1'}).add_to(m)
m