Voronoi regionalizer
Run this notebook in Google Colab:
Remember to install the srai library before running the notebook:
%pip install srai[all]
In [1]:
Copied!
import geopandas as gpd
import numpy as np
import plotly.express as px
from shapely.geometry import Point
from srai.constants import WGS84_CRS
from srai.plotting.folium_wrapper import plot_regions
from srai.regionalizers import VoronoiRegionalizer, geocode_to_region_gdf
import geopandas as gpd
import numpy as np
import plotly.express as px
from shapely.geometry import Point
from srai.constants import WGS84_CRS
from srai.plotting.folium_wrapper import plot_regions
from srai.regionalizers import VoronoiRegionalizer, geocode_to_region_gdf
Regionalizer whole Earth¶
Basic usage of VoronoiRegionalizer to cover whole Earth using 6 poles.
In [2]:
Copied!
# 6 poles of the Earth
seeds_gdf = gpd.GeoDataFrame(
{
"geometry": [
Point(0, 0),
Point(90, 0),
Point(180, 0),
Point(-90, 0),
Point(0, 90),
Point(0, -90),
]
},
index=[1, 2, 3, 4, 5, 6],
crs=WGS84_CRS,
)
# 6 poles of the Earth
seeds_gdf = gpd.GeoDataFrame(
{
"geometry": [
Point(0, 0),
Point(90, 0),
Point(180, 0),
Point(-90, 0),
Point(0, 90),
Point(0, -90),
]
},
index=[1, 2, 3, 4, 5, 6],
crs=WGS84_CRS,
)
In [3]:
Copied!
seeds_gdf.plot()
seeds_gdf.plot()
Out[3]:
<Axes: >
In [4]:
Copied!
vr = VoronoiRegionalizer(seeds=seeds_gdf)
vr = VoronoiRegionalizer(seeds=seeds_gdf)
In [5]:
Copied!
result_gdf = vr.transform()
result_gdf = vr.transform()
In [6]:
Copied!
result_gdf
result_gdf
Out[6]:
| geometry | |
|---|---|
| region_id | |
| 6 | POLYGON ((45.19105 -35.3542, 45.28673 -35.3989... |
| 4 | POLYGON ((-45.19105 -35.3542, -45.28673 -35.39... |
| 1 | POLYGON ((45 35.08447, 45 34.99451, 45 34.9045... |
| 3 | MULTIPOLYGON (((180 44.82, 180 44.73, 180 44.6... |
| 2 | POLYGON ((135 35.08447, 135 34.99451, 135 34.9... |
| 5 | POLYGON ((134.80895 35.3542, 134.71327 35.3989... |
Globe view¶
In [7]:
Copied!
fig = px.choropleth(
result_gdf,
geojson=result_gdf.geometry,
locations=result_gdf.index,
color=result_gdf.index,
color_continuous_scale=px.colors.sequential.Viridis,
)
fig2 = px.scatter_geo(seeds_gdf, lat=seeds_gdf.geometry.y, lon=seeds_gdf.geometry.x)
fig.update_traces(marker={"opacity": 0.6}, selector=dict(type="choropleth"))
fig.add_trace(fig2.data[0])
fig.update_traces(marker_color="white", marker_size=10, selector=dict(type="scattergeo"))
fig.update_layout(coloraxis_showscale=False)
fig.update_geos(
projection_type="orthographic",
projection_rotation_lon=20,
projection_rotation_lat=30,
showlakes=False,
)
fig.update_layout(height=800, width=800, margin={"r": 0, "t": 0, "l": 0, "b": 0})
fig.show(renderer="png") # replace with fig.show() to allow interactivity
fig = px.choropleth(
result_gdf,
geojson=result_gdf.geometry,
locations=result_gdf.index,
color=result_gdf.index,
color_continuous_scale=px.colors.sequential.Viridis,
)
fig2 = px.scatter_geo(seeds_gdf, lat=seeds_gdf.geometry.y, lon=seeds_gdf.geometry.x)
fig.update_traces(marker={"opacity": 0.6}, selector=dict(type="choropleth"))
fig.add_trace(fig2.data[0])
fig.update_traces(marker_color="white", marker_size=10, selector=dict(type="scattergeo"))
fig.update_layout(coloraxis_showscale=False)
fig.update_geos(
projection_type="orthographic",
projection_rotation_lon=20,
projection_rotation_lat=30,
showlakes=False,
)
fig.update_layout(height=800, width=800, margin={"r": 0, "t": 0, "l": 0, "b": 0})
fig.show(renderer="png") # replace with fig.show() to allow interactivity
2D OSM View¶
In [8]:
Copied!
folium_map = plot_regions(result_gdf)
seeds_gdf.explore(
m=folium_map,
style_kwds=dict(color="#444", opacity=1, fillColor="#f2f2f2", fillOpacity=1),
marker_kwds=dict(radius=3),
)
folium_map = plot_regions(result_gdf)
seeds_gdf.explore(
m=folium_map,
style_kwds=dict(color="#444", opacity=1, fillColor="#f2f2f2", fillOpacity=1),
marker_kwds=dict(radius=3),
)
Out[8]:
Make this Notebook Trusted to load map: File -> Trust Notebook
Regionalize single country¶
Drawing a list of arbitrary points inside of the country boundary and using them for regionalization of the same geometry.
In [9]:
Copied!
uk_gdf = geocode_to_region_gdf(query=["R62149"], by_osmid=True)
uk_shape = uk_gdf.iloc[0].geometry # get the Polygon
uk_gdf = geocode_to_region_gdf(query=["R62149"], by_osmid=True)
uk_shape = uk_gdf.iloc[0].geometry # get the Polygon
In [10]:
Copied!
uk_gdf
uk_gdf
Out[10]:
| geometry | |
|---|---|
| region_id | |
| United Kingdom | MULTIPOLYGON (((-14.01552 57.60263, -14.01459 ... |
In [11]:
Copied!
def generate_random_points(shape, n_points=100):
"""Generates random points."""
minx, miny, maxx, maxy = shape.bounds
pts = []
rng = np.random.default_rng()
while len(pts) < 4:
randx = rng.uniform(minx, maxx, n_points)
randy = rng.uniform(miny, maxy, n_points)
coords = np.vstack((randx, randy)).T
# use only the points inside the geographic area
pts = [p for p in list(map(Point, coords)) if p.within(shape)]
del coords # not used any more
return pts
def generate_random_points(shape, n_points=100):
"""Generates random points."""
minx, miny, maxx, maxy = shape.bounds
pts = []
rng = np.random.default_rng()
while len(pts) < 4:
randx = rng.uniform(minx, maxx, n_points)
randy = rng.uniform(miny, maxy, n_points)
coords = np.vstack((randx, randy)).T
# use only the points inside the geographic area
pts = [p for p in list(map(Point, coords)) if p.within(shape)]
del coords # not used any more
return pts
In [12]:
Copied!
pts = generate_random_points(uk_shape)
uk_seeds_gdf = gpd.GeoDataFrame(
{"geometry": pts},
index=list(range(len(pts))),
crs=WGS84_CRS,
)
pts = generate_random_points(uk_shape)
uk_seeds_gdf = gpd.GeoDataFrame(
{"geometry": pts},
index=list(range(len(pts))),
crs=WGS84_CRS,
)
Random points on a map¶
In [13]:
Copied!
folium_map = plot_regions(uk_gdf, tiles_style="CartoDB positron")
uk_seeds_gdf.explore(
m=folium_map,
style_kwds=dict(color="#444", opacity=1, fillColor="#f2f2f2", fillOpacity=1),
marker_kwds=dict(radius=3),
)
folium_map = plot_regions(uk_gdf, tiles_style="CartoDB positron")
uk_seeds_gdf.explore(
m=folium_map,
style_kwds=dict(color="#444", opacity=1, fillColor="#f2f2f2", fillOpacity=1),
marker_kwds=dict(radius=3),
)
Out[13]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [14]:
Copied!
vr_uk = VoronoiRegionalizer(seeds=uk_seeds_gdf)
vr_uk = VoronoiRegionalizer(seeds=uk_seeds_gdf)
In [15]:
Copied!
uk_result_gdf = vr_uk.transform(gdf=uk_gdf)
uk_result_gdf = vr_uk.transform(gdf=uk_gdf)
In [16]:
Copied!
uk_result_gdf.head()
uk_result_gdf.head()
Out[16]:
| geometry | |
|---|---|
| region_id | |
| 22 | POLYGON ((-4.38215 52.43082, -4.23211 52.41257... |
| 19 | POLYGON ((-0.68758 51.72575, -0.51772 51.72312... |
| 10 | POLYGON ((0.15739 52.58628, 0.29963 52.54778, ... |
| 16 | POLYGON ((1.92352 52.3226, 1.8143 52.38719, 1.... |
| 2 | MULTIPOLYGON (((-2.28864 59.0367, -2.46599 59.... |
Generated regions on a map¶
In [17]:
Copied!
folium_map = plot_regions(uk_result_gdf, tiles_style="CartoDB positron")
uk_seeds_gdf.explore(
m=folium_map,
style_kwds=dict(color="#444", opacity=1, fillColor="#f2f2f2", fillOpacity=1),
marker_kwds=dict(radius=3),
)
folium_map = plot_regions(uk_result_gdf, tiles_style="CartoDB positron")
uk_seeds_gdf.explore(
m=folium_map,
style_kwds=dict(color="#444", opacity=1, fillColor="#f2f2f2", fillOpacity=1),
marker_kwds=dict(radius=3),
)
Out[17]:
Make this Notebook Trusted to load map: File -> Trust Notebook
Higher amount of points¶
Example of railway stations in Germany (5000+ seeds) with multiprocessing.
In [18]:
Copied!
stations_csv = gpd.pd.read_csv(
"https://raw.githubusercontent.com/trainline-eu/stations/master/stations.csv",
sep=";",
index_col="id",
usecols=["id", "latitude", "longitude", "country"],
)
stations_csv
stations_csv = gpd.pd.read_csv(
"https://raw.githubusercontent.com/trainline-eu/stations/master/stations.csv",
sep=";",
index_col="id",
usecols=["id", "latitude", "longitude", "country"],
)
stations_csv
Out[18]:
| latitude | longitude | country | |
|---|---|---|---|
| id | |||
| 1 | 44.081790 | 6.001625 | FR |
| 2 | 44.061565 | 5.997373 | FR |
| 3 | 44.063863 | 6.011248 | FR |
| 4 | 44.350000 | 6.350000 | FR |
| 6 | 44.088710 | 6.222982 | FR |
| ... | ... | ... | ... |
| 74704 | 48.265370 | 7.235370 | FR |
| 74705 | 48.262740 | 7.225320 | FR |
| 74706 | 48.794560 | 7.589270 | FR |
| 74707 | 44.360410 | 5.138650 | FR |
| 74708 | 45.604780 | 0.101240 | FR |
70497 rows × 3 columns
In [19]:
Copied!
stations = []
positions = set()
for idx, r in stations_csv.iterrows():
if r.country != "DE" or gpd.pd.isna(r.latitude) or gpd.pd.isna(r.longitude):
continue
pos = round(r.longitude, 5), round(r.latitude, 5)
if pos not in positions:
stations.append({"id": idx, "geometry": Point(*pos)})
positions.add(pos)
stations_gdf = gpd.GeoDataFrame(data=stations, crs=WGS84_CRS).set_index("id")
del stations_csv
del stations
del positions
stations_gdf.head()
stations = []
positions = set()
for idx, r in stations_csv.iterrows():
if r.country != "DE" or gpd.pd.isna(r.latitude) or gpd.pd.isna(r.longitude):
continue
pos = round(r.longitude, 5), round(r.latitude, 5)
if pos not in positions:
stations.append({"id": idx, "geometry": Point(*pos)})
positions.add(pos)
stations_gdf = gpd.GeoDataFrame(data=stations, crs=WGS84_CRS).set_index("id")
del stations_csv
del stations
del positions
stations_gdf.head()
Out[19]:
| geometry | |
|---|---|
| id | |
| 6691 | POINT (9.03081 51.73032) |
| 6692 | POINT (11.80785 49.09284) |
| 6693 | POINT (10.48471 53.14212) |
| 6695 | POINT (6.59351 50.44187) |
| 6696 | POINT (8.31165 49.93009) |
In [20]:
Copied!
vr_rail = VoronoiRegionalizer(seeds=stations_gdf)
vr_rail = VoronoiRegionalizer(seeds=stations_gdf)
In [21]:
Copied!
rail_result_gdf = vr_rail.transform()
rail_result_gdf = vr_rail.transform()
Germany view¶
In [22]:
Copied!
folium_map = plot_regions(rail_result_gdf, tiles_style="CartoDB positron")
stations_gdf.explore(
m=folium_map,
style_kwds=dict(color="#444", opacity=1, fillColor="#f2f2f2", fillOpacity=1),
marker_kwds=dict(radius=1),
)
folium_map.fit_bounds([(54.98310, 5.98865), (47.30248, 15.01699)])
folium_map
folium_map = plot_regions(rail_result_gdf, tiles_style="CartoDB positron")
stations_gdf.explore(
m=folium_map,
style_kwds=dict(color="#444", opacity=1, fillColor="#f2f2f2", fillOpacity=1),
marker_kwds=dict(radius=1),
)
folium_map.fit_bounds([(54.98310, 5.98865), (47.30248, 15.01699)])
folium_map
Out[22]:
Make this Notebook Trusted to load map: File -> Trust Notebook
Berlin view¶
In [23]:
Copied!
# Berlin
folium_map.fit_bounds([(52.51637 + 0.1, 13.40665 + 0.1), (52.51637 - 0.1, 13.40665 - 0.1)])
folium_map
# Berlin
folium_map.fit_bounds([(52.51637 + 0.1, 13.40665 + 0.1), (52.51637 - 0.1, 13.40665 - 0.1)])
folium_map
Out[23]:
Make this Notebook Trusted to load map: File -> Trust Notebook
Difference between spherical voronoi and 2d voronoi¶
Showing the difference between working on the sphere and projected 2D plane.
Uses shapely.voronoi_polygons function as an example.
shapely.voronoi_diagram function allows for a quick division of the Earth using list of seeds on a projected 2d plane.
This results in straight lines with angles distorted and polygons differences
might be substantial during comparisons or any area calculations.
In [24]:
Copied!
from plotly.subplots import make_subplots
from shapely.ops import voronoi_diagram
from plotly.subplots import make_subplots
from shapely.ops import voronoi_diagram
In [25]:
Copied!
pl_gdf = geocode_to_region_gdf(query=["R49715"], by_osmid=True)
pl_gdf_shape = pl_gdf.iloc[0].geometry # get the Polygon
pl_gdf = geocode_to_region_gdf(query=["R49715"], by_osmid=True)
pl_gdf_shape = pl_gdf.iloc[0].geometry # get the Polygon
In [26]:
Copied!
pts = generate_random_points(pl_gdf_shape)
pl_seeds_gdf = gpd.GeoDataFrame(
{"geometry": pts},
index=list(range(len(pts))),
crs=WGS84_CRS,
)
pts = generate_random_points(pl_gdf_shape)
pl_seeds_gdf = gpd.GeoDataFrame(
{"geometry": pts},
index=list(range(len(pts))),
crs=WGS84_CRS,
)
In [27]:
Copied!
region_polygons = list(
voronoi_diagram(pl_seeds_gdf.geometry.union_all(), envelope=pl_gdf_shape).normalize().geoms
)
region_polygons = list(
voronoi_diagram(pl_seeds_gdf.geometry.union_all(), envelope=pl_gdf_shape).normalize().geoms
)
In [28]:
Copied!
pl_regions_2d_gdf = gpd.GeoDataFrame(
{"geometry": [polygon for polygon in region_polygons]},
index=list(range(len(region_polygons))),
crs=WGS84_CRS,
).clip(pl_gdf_shape)
pl_regions_2d_gdf = gpd.GeoDataFrame(
{"geometry": [polygon for polygon in region_polygons]},
index=list(range(len(region_polygons))),
crs=WGS84_CRS,
).clip(pl_gdf_shape)
In [29]:
Copied!
pl_regions_2d_gdf
pl_regions_2d_gdf
Out[29]:
| geometry | |
|---|---|
| 1 | POLYGON ((18.39685 50.03854, 18.77187 50.2928,... |
| 35 | POLYGON ((15.7816 50.95991, 15.92487 51.09869,... |
| 57 | POLYGON ((15.33228 51.78947, 15.7816 50.95991,... |
| 52 | POLYGON ((17.1352 51.30507, 17.3558 51.06539, ... |
| 18 | POLYGON ((16.27108 51.85955, 16.54277 51.84197... |
| ... | ... |
| 27 | POLYGON ((19.36285 54.50397, 19.10472 54.59581... |
| 36 | MULTIPOLYGON (((16.02584 53.92859, 15.96071 53... |
| 67 | POLYGON ((16.21444 54.14815, 16.02584 53.92859... |
| 54 | POLYGON ((16.62739 54.22438, 16.21444 54.14815... |
| 19 | POLYGON ((18.09808 54.85238, 17.82177 54.21821... |
70 rows × 1 columns
In [30]:
Copied!
vr_pl = VoronoiRegionalizer(seeds=pl_seeds_gdf)
vr_pl = VoronoiRegionalizer(seeds=pl_seeds_gdf)
In [31]:
Copied!
pl_result_gdf = vr_pl.transform(gdf=pl_gdf)
pl_result_gdf = vr_pl.transform(gdf=pl_gdf)
In [32]:
Copied!
pl_result_gdf
pl_result_gdf
Out[32]:
| geometry | |
|---|---|
| region_id | |
| 55 | POLYGON ((21.48492 49.42813, 21.54008 49.51431... |
| 15 | POLYGON ((18.63742 50.09765, 18.77838 50.06944... |
| 24 | POLYGON ((15.51322 51.6003, 15.62515 51.52025,... |
| 66 | POLYGON ((16.75356 50.94624, 16.78543 50.85662... |
| 47 | POLYGON ((23.55606 50.31367, 23.46682 50.38764... |
| ... | ... |
| 9 | POLYGON ((19.01825 54.19165, 18.87413 54.24989... |
| 31 | POLYGON ((16.7588 54.17437, 16.60083 54.20292,... |
| 3 | POLYGON ((18.49371 54.99914, 18.61772 54.93081... |
| 43 | POLYGON ((21.25024 53.71767, 21.09711 53.6766,... |
| 59 | POLYGON ((19.30525 54.07465, 19.31939 54.17365... |
70 rows × 1 columns