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: 100%|██████████| 6/6 [00:00<00:00, 11.61it/s]
Interpolating edges: 100%|██████████| 36/36 [00:00<00:00, 89.05it/s]
Generating polygons: 100%|██████████| 6/6 [00:00<00:00, 63.51it/s]

In [6]:

result_gdf

result_gdf

Out[6]:

	geometry
region_id
6	POLYGON ((0.12722 -44.99993, 0.25444 -44.99972...
4	POLYGON ((-45.19105 -35.35420, -45.28673 -35.3...
1	POLYGON ((45.00000 35.08447, 45.00000 34.99451...
3	MULTIPOLYGON (((180.00000 44.82000, 180.00000 ...
2	POLYGON ((45.00000 0.17992, 45.00000 0.26988, ...
5	POLYGON ((180.00000 90.00000, 180.00000 89.910...

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%|██████████| 34/34 [00:01<00:00, 25.51it/s]
Interpolating edges: 100%|██████████| 156/156 [00:00<00:00, 1332.08it/s]
Generating polygons: 100%|██████████| 34/34 [00:00<00:00, 291.08it/s]

In [16]:

uk_result_gdf.head()

uk_result_gdf.head()

Out[16]:

	geometry
region_id
4	POLYGON ((-5.26763 50.77545, -5.12872 50.79540...
19	POLYGON ((0.46236 50.64227, 0.45521 50.73250, ...
7	POLYGON ((-3.11695 59.54729, -3.11499 59.54800...
0	POLYGON ((-2.03762 51.75725, -1.91914 51.69067...
26	POLYGON ((0.42633 51.09342, 0.43359 51.00319, ...

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
...	...	...	...
74630	37.332730	-5.964503	ES
74631	44.506141	11.343406	IT
74632	45.545538	10.168659	IT
74633	46.174047	8.819575	CH
74634	37.501551	15.068982	IT

70492 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%|██████████| 13616/13616 [00:03<00:00, 4508.77it/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
47	MULTIPOLYGON (((19.97574 49.97716, 20.37316 49...
13	POLYGON ((19.95916 50.51218, 20.20137 50.74388...
29	POLYGON ((19.06819 50.14904, 19.35182 50.27707...
64	POLYGON ((19.07297 51.05189, 19.35182 50.27707...
0	POLYGON ((18.61961 51.20290, 18.88098 51.28519...
...	...
15	POLYGON ((19.98286 54.11059, 19.41376 54.44465...
67	POLYGON ((17.73715 54.58652, 17.93197 54.43082...
57	POLYGON ((16.44775 54.35495, 15.94471 53.87937...
52	POLYGON ((16.86776 54.84052, 16.87075 54.84154...
39	POLYGON ((16.86375 54.54525, 16.44775 54.35495...

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:01<00:00, 42.90it/s]
Interpolating edges: 100%|██████████| 289/289 [00:00<00:00, 1553.75it/s]
Generating polygons: 100%|██████████| 71/71 [00:00<00:00, 248.23it/s]

In [32]:

pl_result_gdf

pl_result_gdf

Out[32]:

	geometry
region_id
49	POLYGON ((22.48128 49.79983, 22.60510 49.75586...
40	POLYGON ((18.60553 49.82075, 18.72246 49.87509...
38	POLYGON ((22.69840 49.92818, 22.74531 50.01431...
27	POLYGON ((15.57958 50.81010, 15.69494 50.86457...
0	POLYGON ((23.43301 50.27011, 23.34113 50.33833...
...	...
50	POLYGON ((17.90225 54.35682, 17.73849 54.32159...
15	POLYGON ((17.43927 54.88357, 17.52817 54.79916...
3	POLYGON ((20.02075 53.64228, 20.01318 53.73756...
39	POLYGON ((20.59690 53.69564, 20.61576 53.78841...
62	POLYGON ((22.79659 54.32714, 22.77997 54.23720...

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