Voronoi regionalizer
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()
Generating spherical polygons: 0%| | 0/6 [00:00<?, ?it/s]
Generating spherical polygons: 17%|████████████████████████▌ | 1/6 [00:00<00:00, 7.51it/s]
Generating spherical polygons: 33%|█████████████████████████████████████████████████ | 2/6 [00:00<00:00, 7.79it/s]
Generating spherical polygons: 50%|█████████████████████████████████████████████████████████████████████████▌ | 3/6 [00:00<00:00, 8.07it/s]
Generating spherical polygons: 67%|██████████████████████████████████████████████████████████████████████████████████████████████████ | 4/6 [00:00<00:00, 8.15it/s]
Generating spherical polygons: 83%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▌ | 5/6 [00:00<00:00, 8.16it/s]
Generating spherical polygons: 100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 6/6 [00:00<00:00, 8.14it/s]
Generating spherical polygons: 100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 6/6 [00:00<00:00, 8.06it/s]
Interpolating edges: 0%| | 0/36 [00:00<?, ?it/s]
Interpolating edges: 3%|████▎ | 1/36 [00:00<00:16, 2.16it/s]
Interpolating edges: 100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 36/36 [00:00<00:00, 68.95it/s]
Generating polygons: 0%| | 0/6 [00:00<?, ?it/s]
Generating polygons: 83%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 5/6 [00:00<00:00, 45.19it/s]
Generating polygons: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 6/6 [00:00<00:00, 45.21it/s]
In [6]:
Copied!
result_gdf
result_gdf
Out[6]:
| geometry | |
|---|---|
| region_id | |
| 6 | POLYGON ((45.19105 -35.35420, 45.28673 -35.398... |
| 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 ((134.80895 35.35420, 134.71327 35.398... |
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)
Generating spherical polygons: 0%| | 0/30 [00:00<?, ?it/s]
Generating spherical polygons: 3%|████▊ | 1/30 [00:00<00:03, 7.80it/s]
Generating spherical polygons: 20%|█████████████████████████████▏ | 6/30 [00:00<00:00, 29.69it/s]
Generating spherical polygons: 33%|████████████████████████████████████████████████▎ | 10/30 [00:00<00:00, 21.02it/s]
Generating spherical polygons: 43%|██████████████████████████████████████████████████████████████▊ | 13/30 [00:00<00:00, 17.11it/s]
Generating spherical polygons: 53%|█████████████████████████████████████████████████████████████████████████████▎ | 16/30 [00:00<00:00, 17.02it/s]
Generating spherical polygons: 60%|███████████████████████████████████████████████████████████████████████████████████████ | 18/30 [00:01<00:00, 13.95it/s]
Generating spherical polygons: 70%|█████████████████████████████████████████████████████████████████████████████████████████████████████▌ | 21/30 [00:01<00:00, 16.02it/s]
Generating spherical polygons: 77%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████▏ | 23/30 [00:01<00:00, 10.69it/s]
Generating spherical polygons: 87%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▋ | 26/30 [00:01<00:00, 10.88it/s]
Generating spherical polygons: 97%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▏ | 29/30 [00:02<00:00, 12.66it/s]
Generating spherical polygons: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 30/30 [00:02<00:00, 13.71it/s]
Interpolating edges: 0%| | 0/162 [00:00<?, ?it/s]
Interpolating edges: 48%|█████████████████████████████████████████████████████████████████████████▋ | 78/162 [00:00<00:00, 776.17it/s]
Interpolating edges: 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 162/162 [00:00<00:00, 815.78it/s]
Generating polygons: 0%| | 0/30 [00:00<?, ?it/s]
Generating polygons: 70%|███████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 21/30 [00:00<00:00, 175.08it/s]
Generating polygons: 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 30/30 [00:00<00:00, 164.41it/s]
In [16]:
Copied!
uk_result_gdf.head()
uk_result_gdf.head()
Out[16]:
| geometry | |
|---|---|
| region_id | |
| 20 | POLYGON ((-2.41138 51.18386, -2.26822 51.19562... |
| 6 | POLYGON ((-4.50696 50.16602, -4.49602 50.25580... |
| 21 | POLYGON ((-5.75825 51.76373, -5.63100 51.71710... |
| 16 | POLYGON ((0.10761 51.24626, 0.21489 51.17704, ... |
| 29 | POLYGON ((0.66406 50.71132, 0.62114 50.79716, ... |
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 |
| ... | ... | ... | ... |
| 74667 | 41.779016 | 15.445426 | IT |
| 74668 | 48.446400 | 1.476700 | FR |
| 74669 | 43.848670 | -1.358814 | FR |
| 74670 | 48.779036 | -3.238814 | FR |
| 74671 | 43.789600 | -1.395040 | FR |
70517 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()
Generating polygons: 0%| | 0/13617 [00:00<?, ?it/s]
Generating polygons: 3%|███▊ | 349/13617 [00:00<00:03, 3484.07it/s]
Generating polygons: 5%|███████▋ | 698/13617 [00:00<00:03, 3272.65it/s]
Generating polygons: 8%|███████████▏ | 1027/13617 [00:00<00:04, 2869.38it/s]
Generating polygons: 10%|███████████████ | 1381/13617 [00:00<00:03, 3106.54it/s]
Generating polygons: 13%|██████████████████▊ | 1730/13617 [00:00<00:03, 3234.17it/s]
Generating polygons: 15%|██████████████████████▋ | 2082/13617 [00:00<00:03, 3325.06it/s]
Generating polygons: 18%|██████████████████████████▎ | 2425/13617 [00:00<00:03, 3356.73it/s]
Generating polygons: 20%|██████████████████████████████ | 2764/13617 [00:00<00:03, 3333.21it/s]
Generating polygons: 23%|█████████████████████████████████▋ | 3102/13617 [00:00<00:03, 3344.39it/s]
Generating polygons: 25%|█████████████████████████████████████▎ | 3438/13617 [00:01<00:03, 3130.48it/s]
Generating polygons: 28%|████████████████████████████████████████▊ | 3755/13617 [00:01<00:03, 2819.26it/s]
Generating polygons: 30%|███████████████████████████████████████████▉ | 4044/13617 [00:01<00:03, 2800.59it/s]
Generating polygons: 32%|███████████████████████████████████████████████▋ | 4386/13617 [00:01<00:03, 2969.99it/s]
Generating polygons: 35%|███████████████████████████████████████████████████▍ | 4736/13617 [00:01<00:02, 3117.89it/s]
Generating polygons: 37%|███████████████████████████████████████████████████████▏ | 5077/13617 [00:01<00:02, 3200.44it/s]
Generating polygons: 40%|██████████████████████████████████████████████████████████▉ | 5423/13617 [00:01<00:02, 3273.47it/s]
Generating polygons: 42%|██████████████████████████████████████████████████████████████▋ | 5767/13617 [00:01<00:02, 3319.86it/s]
Generating polygons: 45%|██████████████████████████████████████████████████████████████████▍ | 6116/13617 [00:01<00:02, 3369.91it/s]
Generating polygons: 47%|██████████████████████████████████████████████████████████████████████▏ | 6455/13617 [00:02<00:02, 3272.38it/s]
Generating polygons: 50%|█████████████████████████████████████████████████████████████████████████▋ | 6784/13617 [00:02<00:02, 3069.55it/s]
Generating polygons: 52%|█████████████████████████████████████████████████████████████████████████████▏ | 7097/13617 [00:02<00:02, 3083.51it/s]
Generating polygons: 54%|████████████████████████████████████████████████████████████████████████████████▌ | 7408/13617 [00:02<00:02, 2608.89it/s]
Generating polygons: 57%|████████████████████████████████████████████████████████████████████████████████████▎ | 7756/13617 [00:02<00:02, 2831.98it/s]
Generating polygons: 59%|████████████████████████████████████████████████████████████████████████████████████████ | 8101/13617 [00:02<00:01, 2996.56it/s]
Generating polygons: 62%|███████████████████████████████████████████████████████████████████████████████████████████▋ | 8434/13617 [00:02<00:01, 3087.14it/s]
Generating polygons: 65%|███████████████████████████████████████████████████████████████████████████████████████████████▌ | 8791/13617 [00:02<00:01, 3222.44it/s]
Generating polygons: 67%|███████████████████████████████████████████████████████████████████████████████████████████████████▏ | 9121/13617 [00:02<00:01, 3157.98it/s]
Generating polygons: 69%|██████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 9463/13617 [00:03<00:01, 3232.82it/s]
Generating polygons: 72%|██████████████████████████████████████████████████████████████████████████████████████████████████████████▌ | 9809/13617 [00:03<00:01, 3298.43it/s]
Generating polygons: 75%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 10173/13617 [00:03<00:01, 3396.97it/s]
Generating polygons: 77%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████▋ | 10527/13617 [00:03<00:00, 3439.06it/s]
Generating polygons: 80%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▍ | 10873/13617 [00:03<00:00, 3421.92it/s]
Generating polygons: 82%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▏ | 11226/13617 [00:03<00:00, 3451.06it/s]
Generating polygons: 85%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▉ | 11572/13617 [00:03<00:00, 3424.81it/s]
Generating polygons: 88%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▋ | 11916/13617 [00:03<00:00, 3085.78it/s]
Generating polygons: 90%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▎ | 12258/13617 [00:03<00:00, 3177.38it/s]
Generating polygons: 92%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 12582/13617 [00:03<00:00, 3124.04it/s]
Generating polygons: 95%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▎ | 12904/13617 [00:04<00:00, 3150.84it/s]
Generating polygons: 97%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████ | 13258/13617 [00:04<00:00, 3187.62it/s]
Generating polygons: 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▋| 13589/13617 [00:04<00:00, 3191.34it/s]
Generating polygons: 100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 13617/13617 [00:04<00:00, 3148.13it/s]
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.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]:
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 | |
|---|---|
| 62 | POLYGON ((17.72716 50.20484, 17.72725 50.20479... |
| 10 | POLYGON ((18.86746 50.59649, 18.89623 50.65140... |
| 11 | POLYGON ((17.71513 50.91523, 17.94013 51.29233... |
| 27 | POLYGON ((18.76317 51.55247, 19.43041 51.29638... |
| 50 | POLYGON ((16.90944 51.58201, 17.71513 50.91523... |
| ... | ... |
| 18 | POLYGON ((14.88752 53.53729, 14.81359 53.42514... |
| 4 | POLYGON ((15.67858 54.09567, 15.47711 53.72779... |
| 17 | POLYGON ((17.60918 54.43117, 17.29406 54.59986... |
| 22 | POLYGON ((22.09392 53.75768, 20.81480 54.26650... |
| 0 | POLYGON ((19.29098 54.56342, 18.93009 54.18890... |
68 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)
Generating spherical polygons: 0%| | 0/68 [00:00<?, ?it/s]
Generating spherical polygons: 4%|██████▍ | 3/68 [00:00<00:03, 18.10it/s]
Generating spherical polygons: 7%|██████████▋ | 5/68 [00:00<00:04, 14.51it/s]
Generating spherical polygons: 10%|███████████████ | 7/68 [00:00<00:05, 11.51it/s]
Generating spherical polygons: 15%|█████████████████████▎ | 10/68 [00:00<00:05, 10.91it/s]
Generating spherical polygons: 21%|█████████████████████████████▊ | 14/68 [00:01<00:03, 14.01it/s]
Generating spherical polygons: 25%|████████████████████████████████████▎ | 17/68 [00:01<00:03, 14.56it/s]
Generating spherical polygons: 32%|██████████████████████████████████████████████▉ | 22/68 [00:01<00:02, 17.99it/s]
Generating spherical polygons: 41%|███████████████████████████████████████████████████████████▋ | 28/68 [00:01<00:02, 15.94it/s]
Generating spherical polygons: 44%|███████████████████████████████████████████████████████████████▉ | 30/68 [00:02<00:02, 15.06it/s]
Generating spherical polygons: 51%|██████████████████████████████████████████████████████████████████████████▋ | 35/68 [00:02<00:01, 19.88it/s]
Generating spherical polygons: 65%|█████████████████████████████████████████████████████████████████████████████████████████████▊ | 44/68 [00:02<00:00, 29.69it/s]
Generating spherical polygons: 71%|██████████████████████████████████████████████████████████████████████████████████████████████████████▎ | 48/68 [00:02<00:00, 27.92it/s]
Generating spherical polygons: 91%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▏ | 62/68 [00:02<00:00, 44.99it/s]
Generating spherical polygons: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 68/68 [00:02<00:00, 25.87it/s]
Interpolating edges: 0%| | 0/290 [00:00<?, ?it/s]
Interpolating edges: 38%|█████████████████████████████████████████████████████████▊ | 111/290 [00:00<00:00, 1106.29it/s]
Interpolating edges: 77%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████▌ | 222/290 [00:00<00:00, 1038.08it/s]
Interpolating edges: 100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 290/290 [00:00<00:00, 1031.66it/s]
Generating polygons: 0%| | 0/68 [00:00<?, ?it/s]
Generating polygons: 32%|█████████████████████████████████████████████████▊ | 22/68 [00:00<00:00, 202.12it/s]
Generating polygons: 91%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▍ | 62/68 [00:00<00:00, 305.40it/s]
Generating polygons: 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 68/68 [00:00<00:00, 309.08it/s]
In [32]:
Copied!
pl_result_gdf
pl_result_gdf
Out[32]:
| geometry | |
|---|---|
| region_id | |
| 27 | MULTIPOLYGON (((20.93369 49.32803, 21.05148 49... |
| 21 | MULTIPOLYGON (((17.67506 50.71827, 17.70398 50... |
| 2 | POLYGON ((14.86920 51.64570, 15.00040 51.68406... |
| 34 | POLYGON ((14.68473 52.30921, 14.84763 52.29609... |
| 61 | POLYGON ((23.43262 51.49468, 23.57614 51.51055... |
| ... | ... |
| 0 | POLYGON ((19.17348 54.43600, 19.27014 54.36134... |
| 40 | POLYGON ((20.60519 54.30506, 20.57908 54.20877... |
| 15 | MULTIPOLYGON (((14.78411 53.70234, 14.69366 53... |
| 20 | POLYGON ((16.04185 53.87793, 15.87085 53.93648... |
| 67 | POLYGON ((16.36318 54.04798, 16.29492 54.13301... |
68 rows × 1 columns
In [33]:
Copied!
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
