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
from srai.constants import WGS84_CRS
from srai.plotting.folium_wrapper import plot_regions
from srai.utils import geocode_to_region_gdf

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

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 regions:   0%|          | 0/6 [00:00<?, ?it/s]/opt/hostedtoolcache/Python/3.10.12/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.12/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 regions:  50%|█████     | 3/6 [00:02<00:02,  1.38it/s]/opt/hostedtoolcache/Python/3.10.12/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 regions:  67%|██████▋   | 4/6 [00:02<00:01,  1.59it/s]/opt/hostedtoolcache/Python/3.10.12/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 regions:  83%|████████▎ | 5/6 [00:03<00:00,  1.73it/s]/opt/hostedtoolcache/Python/3.10.12/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 regions: 100%|██████████| 6/6 [00:03<00:00,  1.58it/s]

In [6]:

result_gdf

result_gdf

Out[6]:

	geometry
region_id
6	POLYGON ((135.19105 -35.35420, 135.28673 -35.3...
4	POLYGON ((-45.19105 -35.35420, -45.28673 -35.3...
1	POLYGON ((0.25444 44.99972, 0.38167 44.99936, ...
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 regions: 100%|██████████| 34/34 [00:06<00:00,  5.44it/s]

In [16]:

uk_result_gdf.head()

uk_result_gdf.head()

Out[16]:

	geometry
region_id
28	POLYGON ((-6.53427 49.72894, -6.40939 49.76881...
18	POLYGON ((1.62505 52.34082, 1.77094 52.35344, ...
7	POLYGON ((0.31510 52.61285, 0.33914 52.70554, ...
11	MULTIPOLYGON (((-3.29934 59.40946, -3.28784 59...
1	MULTIPOLYGON (((-1.53635 57.56700, -1.70691 57...

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
...	...	...	...
68177	44.674312	10.478726	IT
68178	44.661871	10.466608	IT
68179	44.694544	10.498249	IT
68180	43.125000	5.972400	FR
68181	40.726753	13.874011	IT

64039 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()

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
23	POLYGON ((19.58295 49.71063, 19.82752 49.99465...
5	POLYGON ((22.03071 50.17879, 23.45705 51.33318...
19	POLYGON ((20.88561 50.27163, 21.10690 50.71567...
21	POLYGON ((19.77365 50.71436, 19.80760 50.74745...
47	POLYGON ((20.13700 51.19459, 20.91696 50.94306...
...	...
3	POLYGON ((17.03751 54.67573, 16.77948 54.47242...
56	POLYGON ((17.12476 54.91886, 17.12565 54.91906...
52	POLYGON ((18.25607 55.03479, 18.27736 55.03561...
54	POLYGON ((18.25230 54.13957, 17.80361 54.00549...
11	POLYGON ((19.03971 54.02020, 19.03385 54.01985...

70 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 regions: 100%|██████████| 70/70 [00:08<00:00,  8.36it/s]

In [32]:

pl_result_gdf

pl_result_gdf

Out[32]:

	geometry
region_id
30	POLYGON ((20.16292 49.84101, 20.30190 49.81972...
41	MULTIPOLYGON (((18.78655 49.71630, 18.90205 49...
37	POLYGON ((21.93921 49.35217, 21.94965 49.44653...
61	POLYGON ((14.80331 51.73169, 14.95400 51.72188...
4	POLYGON ((14.36944 53.65261, 14.51893 53.63680...
...	...
32	POLYGON ((17.42129 54.52855, 17.49741 54.44699...
29	POLYGON ((19.14425 54.01089, 18.97912 53.93991...
22	POLYGON ((17.72397 54.20204, 17.64875 54.28374...
65	POLYGON ((20.43823 53.98743, 20.27657 53.99113...
13	POLYGON ((19.50587 54.48265, 19.37990 54.42330...

70 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