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 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.57it/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.70it/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 ((0.12722 -44.99993, 0.25444 -44.99972...
4	POLYGON ((-135.00000 -35.08447, -135.00000 -34...
1	POLYGON ((0.12722 44.99993, 0.25444 44.99972, ...
3	MULTIPOLYGON (((135.00000 0.17992, 135.00000 0...
2	POLYGON ((45.00000 0.08996, 45.00000 0.17992, ...
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 regions: 100%|██████████| 24/24 [00:06<00:00,  3.70it/s]

In [16]:

uk_result_gdf.head()

uk_result_gdf.head()

Out[16]:

	geometry
region_id
6	POLYGON ((-5.29365 50.99122, -5.15432 50.96423...
11	POLYGON ((0.13422 51.49328, 0.26825 51.45972, ...
20	POLYGON ((0.14326 52.26937, 0.28664 52.28930, ...
3	POLYGON ((2.01703 52.51474, 1.87217 52.49692, ...
16	MULTIPOLYGON (((-5.72215 57.24859, -5.76929 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
...	...	...	...
68186	46.169020	10.763700	IT
68187	47.369300	-1.177800	FR
68188	44.870900	-0.614000	FR
68189	41.224651	16.314623	IT
68190	45.234444	15.296389	HR

64048 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
31	POLYGON ((18.57741 50.16138, 18.87397 50.40943...
24	POLYGON ((20.20662 50.84976, 20.21139 50.84962...
22	POLYGON ((21.03307 51.08242, 21.83410 50.60255...
44	POLYGON ((18.58816 51.21469, 18.88671 50.71613...
43	POLYGON ((18.88671 50.71613, 19.46110 51.11189...
...	...
60	POLYGON ((19.19571 54.70446, 19.40174 54.60438...
14	POLYGON ((18.26959 54.43121, 17.94014 54.05019...
15	MULTIPOLYGON (((15.10605 53.71909, 15.01669 53...
21	POLYGON ((22.15032 53.92536, 21.94892 53.93306...
9	POLYGON ((18.44605 55.02332, 18.44710 55.02311...

63 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%|██████████| 63/63 [00:09<00:00,  6.59it/s]

In [32]:

pl_result_gdf

pl_result_gdf

Out[32]:

	geometry
region_id
31	MULTIPOLYGON (((20.45569 49.88224, 20.56089 49...
27	POLYGON ((18.50653 49.97701, 18.56751 50.05802...
13	MULTIPOLYGON (((15.09098 51.38437, 15.22261 51...
56	MULTIPOLYGON (((22.62949 50.25347, 22.75025 50...
21	POLYGON ((14.97739 51.35242, 14.97757 51.35305...
...	...
8	POLYGON ((20.74544 53.39756, 20.61375 53.32776...
16	POLYGON ((22.19171 54.09386, 22.32900 54.17035...
46	POLYGON ((21.63514 54.30617, 21.74064 54.23407...
17	POLYGON ((19.35177 54.57781, 19.19978 54.54693...
39	POLYGON ((22.60509 54.32286, 22.46679 54.24668...

63 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