Voronoi regionalizer

In [1]:

import geopandas as gpd
import numpy as np
import plotly.express as px
from shapely.geometry import Point

from srai.regionalizers import VoronoiRegionalizer, geocode_to_region_gdf
from srai.constants import WGS84_CRS
from srai.plotting.folium_wrapper import plot_regions

import geopandas as gpd import numpy as np import plotly.express as px from shapely.geometry import Point from srai.regionalizers import VoronoiRegionalizer, geocode_to_region_gdf from srai.constants import WGS84_CRS from srai.plotting.folium_wrapper import plot_regions

Regionalizer whole Earth¶

Basic usage of VoronoiRegionalizer to cover whole Earth using 6 poles.

In [2]:

# 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]:

seeds_gdf.plot()

seeds_gdf.plot()

Out[3]:

<Axes: >

In [4]:

vr = VoronoiRegionalizer(seeds=seeds_gdf)

vr = VoronoiRegionalizer(seeds=seeds_gdf)

In [5]:

result_gdf = vr.transform()

result_gdf = vr.transform()

Generating spherical polygons:   0%|          | 0/6 [00:00<?, ?it/s]/opt/hostedtoolcache/Python/3.10.13/x64/lib/python3.10/site-packages/spherical_geometry/great_circle_arc.py:365: RuntimeWarning: invalid value encountered in divide
  return P / l
/opt/hostedtoolcache/Python/3.10.13/x64/lib/python3.10/site-packages/spherical_geometry/great_circle_arc.py:261: RuntimeWarning: invalid value encountered in intersects
  return math_util.intersects(A, B, C, D)
Generating spherical polygons: 100%|██████████| 6/6 [00:00<00:00,  7.86it/s]
Interpolating edges: 100%|██████████| 36/36 [00:00<00:00, 50.89it/s]
Generating polygons: 100%|██████████| 6/6 [00:00<00:00, 36.37it/s]

In [6]:

result_gdf

result_gdf

Out[6]:

	geometry
region_id
6	POLYGON ((180.00000 -45.18000, 180.00000 -45.2...
4	POLYGON ((-45.19105 -35.35420, -45.28673 -35.3...
1	POLYGON ((0.12722 44.99993, 0.25444 44.99972, ...
3	MULTIPOLYGON (((135.00000 0.17992, 135.00000 0...
2	POLYGON ((45.19105 35.35420, 45.28673 35.39899...
5	POLYGON ((180.00000 89.82000, 180.00000 89.730...

Globe view¶

In [7]:

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

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

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

uk_gdf

uk_gdf

Out[10]:

	geometry
region_id
United Kingdom	MULTIPOLYGON (((-14.01552 57.60263, -14.01459 ...

In [11]:

def generate_random_points(shape, n_points=100):
    minx, miny, maxx, maxy = shape.bounds
    pts = []

    while len(pts) < 4:
        randx = np.random.uniform(minx, maxx, n_points)
        randy = np.random.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): minx, miny, maxx, maxy = shape.bounds pts = [] while len(pts) < 4: randx = np.random.uniform(minx, maxx, n_points) randy = np.random.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]:

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

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

vr_uk = VoronoiRegionalizer(seeds=uk_seeds_gdf)

vr_uk = VoronoiRegionalizer(seeds=uk_seeds_gdf)

In [15]:

uk_result_gdf = vr_uk.transform(gdf=uk_gdf)

uk_result_gdf = vr_uk.transform(gdf=uk_gdf)

Generating spherical polygons: 100%|██████████| 29/29 [00:02<00:00, 13.68it/s]
Interpolating edges: 100%|██████████| 155/155 [00:00<00:00, 578.76it/s]
Generating polygons: 100%|██████████| 29/29 [00:00<00:00, 113.35it/s]

In [16]:

uk_result_gdf.head()

uk_result_gdf.head()

Out[16]:

	geometry
region_id
18	POLYGON ((-3.78107 51.14743, -3.63704 51.15530...
15	POLYGON ((-5.27678 51.18429, -5.17750 51.24928...
16	POLYGON ((-1.49198 52.49609, -1.36040 52.54685...
4	POLYGON ((-0.09859 52.30405, -0.19755 52.38666...
14	POLYGON ((0.11737 52.91652, 0.23505 52.97238, ...

Generated regions on a map¶

In [17]:

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

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
...	...	...	...
68224	40.484794	-3.662886	ES
68225	40.482000	-3.616292	ES
68226	40.598525	-3.715582	ES
68227	40.341551	-3.711546	ES
68228	40.381184	-3.780274	ES

64086 rows × 3 columns

In [19]:

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 not pos 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 not pos 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]:

vr_rail = VoronoiRegionalizer(seeds=stations_gdf)

vr_rail = VoronoiRegionalizer(seeds=stations_gdf)

In [21]:

rail_result_gdf = vr_rail.transform()

rail_result_gdf = vr_rail.transform()

Generating polygons: 100%|██████████| 7459/7459 [00:03<00:00, 2245.61it/s]

Germany view¶

In [22]:

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

# 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]:

from shapely.ops import voronoi_diagram
from plotly.subplots import make_subplots

from shapely.ops import voronoi_diagram from plotly.subplots import make_subplots

In [25]:

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

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

region_polygons = list(
    voronoi_diagram(pl_seeds_gdf.geometry.unary_union, envelope=pl_gdf_shape).normalize().geoms
)

region_polygons = list( voronoi_diagram(pl_seeds_gdf.geometry.unary_union, envelope=pl_gdf_shape).normalize().geoms )

In [28]:

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

pl_regions_2d_gdf

pl_regions_2d_gdf

Out[29]:

	geometry
64	POLYGON ((20.23346 50.10103, 20.67588 49.41200...
14	POLYGON ((19.33705 50.05861, 20.09478 50.35449...
66	POLYGON ((17.99117 50.72687, 18.86129 50.77224...
11	POLYGON ((20.17265 50.62586, 20.35645 50.85685...
26	POLYGON ((19.15913 51.30435, 19.24059 51.30713...
...	...
24	POLYGON ((20.38738 54.17046, 20.11512 54.08874...
47	POLYGON ((19.10351 54.44130, 19.62761 54.28687...
43	POLYGON ((19.73403 54.44730, 19.79734 54.44286...
63	POLYGON ((19.19583 54.70440, 19.40174 54.60438...
2	POLYGON ((19.10351 54.44130, 18.96475 54.27570...

71 rows × 1 columns

In [30]:

vr_pl = VoronoiRegionalizer(seeds=pl_seeds_gdf)

vr_pl = VoronoiRegionalizer(seeds=pl_seeds_gdf)

In [31]:

pl_result_gdf = vr_pl.transform(gdf=pl_gdf)

pl_result_gdf = vr_pl.transform(gdf=pl_gdf)

Generating spherical polygons: 100%|██████████| 71/71 [00:02<00:00, 26.74it/s]
Interpolating edges: 100%|██████████| 299/299 [00:00<00:00, 737.49it/s]
Generating polygons: 100%|██████████| 71/71 [00:00<00:00, 254.41it/s]

In [32]:

pl_result_gdf

pl_result_gdf

Out[32]:

	geometry
region_id
11	POLYGON ((21.84737 49.46957, 21.99174 49.47968...
25	POLYGON ((20.49043 49.72716, 20.59549 49.65862...
31	POLYGON ((16.06226 51.05129, 16.22355 51.04601...
12	POLYGON ((14.78358 51.67186, 14.92180 51.70414...
36	POLYGON ((15.22810 52.48024, 15.13921 52.40585...
...	...
41	POLYGON ((16.96257 54.43394, 17.07769 54.49409...
59	POLYGON ((19.11424 54.54621, 19.19622 54.62466...
24	POLYGON ((23.41292 54.12841, 23.35620 54.03984...
35	POLYGON ((20.17361 54.19695, 20.04841 54.26897...
1	POLYGON ((21.47911 54.30048, 21.43507 54.21408...

71 rows × 1 columns

In [33]:

choropleth_1 = px.choropleth(
    pl_result_gdf,
    geojson=pl_result_gdf.geometry,
    locations=pl_result_gdf.index,
    color=pl_result_gdf.index,
    color_continuous_scale=px.colors.qualitative.Plotly,
)

choropleth_2 = px.choropleth(
    pl_regions_2d_gdf,
    geojson=pl_regions_2d_gdf.geometry,
    locations=pl_regions_2d_gdf.index,
    color=pl_regions_2d_gdf.index,
    color_continuous_scale=px.colors.qualitative.Plotly,
)

points_plot = px.scatter_geo(pl_seeds_gdf, lat=pl_seeds_gdf.geometry.y, lon=pl_seeds_gdf.geometry.x)

fig = make_subplots(
    rows=2,
    cols=2,
    specs=[
        [{"type": "scattergeo"}, {"type": "scattergeo"}],
        [{"type": "scattergeo"}, {"type": "scattergeo"}],
    ],
    horizontal_spacing=0.01,
    vertical_spacing=0.0,
)

fig.add_trace(choropleth_1["data"][0], row=1, col=1)
fig.add_trace(choropleth_1["data"][0], row=2, col=1)
fig.add_trace(choropleth_2["data"][0], row=1, col=2)
fig.add_trace(choropleth_2["data"][0], row=2, col=2)
for r in [1, 2]:
    for c in [1, 2]:
        fig.add_trace(points_plot.data[0], row=r, col=c)

minx, miny, maxx, maxy = pl_gdf_shape.bounds

fig.update_traces(marker_color="black", marker_size=6, selector=dict(type="scattergeo"))
fig.update_layout(coloraxis_showscale=False)
fig.update_geos(
    projection_type="natural earth",
    lataxis_range=[miny - 1, maxy + 1],
    lonaxis_range=[minx - 1, maxx + 1],
    resolution=50,
    row=1,
    showlakes=False,
)

fig.update_geos(
    projection_type="natural earth",
    lataxis_range=[miny + 1, maxy - 1],
    lonaxis_range=[minx + 2, maxx - 2],
    resolution=50,
    row=2,
    showlakes=False,
)

fig.update_traces(marker={"opacity": 0.6}, selector=dict(type="choropleth"), row=1)
fig.update_traces(marker={"opacity": 0.25}, selector=dict(type="choropleth"), row=2)

fig.update_layout(
    height=800,
    width=800,
    title_text="Side By Side Subplots (Left: Spherical voronoi, Right: 2D voronoi)",
    margin={"r": 5, "t": 50, "l": 5, "b": 0},
)
fig.show(renderer="png")  # replace with fig.show() to allow interactivity

choropleth_1 = px.choropleth( pl_result_gdf, geojson=pl_result_gdf.geometry, locations=pl_result_gdf.index, color=pl_result_gdf.index, color_continuous_scale=px.colors.qualitative.Plotly, ) choropleth_2 = px.choropleth( pl_regions_2d_gdf, geojson=pl_regions_2d_gdf.geometry, locations=pl_regions_2d_gdf.index, color=pl_regions_2d_gdf.index, color_continuous_scale=px.colors.qualitative.Plotly, ) points_plot = px.scatter_geo(pl_seeds_gdf, lat=pl_seeds_gdf.geometry.y, lon=pl_seeds_gdf.geometry.x) fig = make_subplots( rows=2, cols=2, specs=[ [{"type": "scattergeo"}, {"type": "scattergeo"}], [{"type": "scattergeo"}, {"type": "scattergeo"}], ], horizontal_spacing=0.01, vertical_spacing=0.0, ) fig.add_trace(choropleth_1["data"][0], row=1, col=1) fig.add_trace(choropleth_1["data"][0], row=2, col=1) fig.add_trace(choropleth_2["data"][0], row=1, col=2) fig.add_trace(choropleth_2["data"][0], row=2, col=2) for r in [1, 2]: for c in [1, 2]: fig.add_trace(points_plot.data[0], row=r, col=c) minx, miny, maxx, maxy = pl_gdf_shape.bounds fig.update_traces(marker_color="black", marker_size=6, selector=dict(type="scattergeo")) fig.update_layout(coloraxis_showscale=False) fig.update_geos( projection_type="natural earth", lataxis_range=[miny - 1, maxy + 1], lonaxis_range=[minx - 1, maxx + 1], resolution=50, row=1, showlakes=False, ) fig.update_geos( projection_type="natural earth", lataxis_range=[miny + 1, maxy - 1], lonaxis_range=[minx + 2, maxx - 2], resolution=50, row=2, showlakes=False, ) fig.update_traces(marker={"opacity": 0.6}, selector=dict(type="choropleth"), row=1) fig.update_traces(marker={"opacity": 0.25}, selector=dict(type="choropleth"), row=2) fig.update_layout( height=800, width=800, title_text="Side By Side Subplots (Left: Spherical voronoi, Right: 2D voronoi)", margin={"r": 5, "t": 50, "l": 5, "b": 0}, ) fig.show(renderer="png") # replace with fig.show() to allow interactivity