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: 33%|███▎ | 2/6 [00:00<00:00, 11.34it/s]
Generating spherical polygons: 67%|██████▋ | 4/6 [00:00<00:00, 11.66it/s]
Generating spherical polygons: 100%|██████████| 6/6 [00:00<00:00, 11.65it/s]
Generating spherical polygons: 100%|██████████| 6/6 [00:00<00:00, 11.61it/s]
Interpolating edges: 0%| | 0/36 [00:00<?, ?it/s]
Interpolating edges: 3%|▎ | 1/36 [00:00<00:13, 2.65it/s]
Interpolating edges: 100%|██████████| 36/36 [00:00<00:00, 85.42it/s]
Generating polygons: 0%| | 0/6 [00:00<?, ?it/s]
Generating polygons: 100%|██████████| 6/6 [00:00<00:00, 63.35it/s]
In [6]:
Copied!
result_gdf
result_gdf
Out[6]:
geometry | |
---|---|
region_id | |
6 | POLYGON ((180 -45.18, 180 -45.27, 180 -45.36, ... |
4 | POLYGON ((-45 -0.08996, -45 -0.17992, -45 -0.2... |
1 | POLYGON ((45 35.08447, 45 34.99451, 45 34.9045... |
3 | MULTIPOLYGON (((135 0.08996, 135 0.17992, 135 ... |
2 | POLYGON ((135 35.08447, 135 34.99451, 135 34.9... |
5 | POLYGON ((44.80895 35.3542, 44.71327 35.39899,... |
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/32 [00:00<?, ?it/s]
Generating spherical polygons: 12%|█▎ | 4/32 [00:00<00:01, 14.44it/s]
Generating spherical polygons: 22%|██▏ | 7/32 [00:00<00:01, 16.16it/s]
Generating spherical polygons: 28%|██▊ | 9/32 [00:00<00:01, 14.08it/s]
Generating spherical polygons: 34%|███▍ | 11/32 [00:00<00:01, 14.33it/s]
Generating spherical polygons: 41%|████ | 13/32 [00:00<00:01, 15.52it/s]
Generating spherical polygons: 66%|██████▌ | 21/32 [00:01<00:00, 28.67it/s]
Generating spherical polygons: 84%|████████▍ | 27/32 [00:01<00:00, 33.89it/s]
Generating spherical polygons: 97%|█████████▋| 31/32 [00:01<00:00, 34.18it/s]
Generating spherical polygons: 100%|██████████| 32/32 [00:01<00:00, 25.49it/s]
Interpolating edges: 0%| | 0/149 [00:00<?, ?it/s]
Interpolating edges: 93%|█████████▎| 139/149 [00:00<00:00, 1388.41it/s]
Interpolating edges: 100%|██████████| 149/149 [00:00<00:00, 1363.45it/s]
Generating polygons: 0%| | 0/32 [00:00<?, ?it/s]
Generating polygons: 44%|████▍ | 14/32 [00:00<00:00, 133.50it/s]
Generating polygons: 100%|██████████| 32/32 [00:00<00:00, 283.15it/s]
In [16]:
Copied!
uk_result_gdf.head()
uk_result_gdf.head()
Out[16]:
geometry | |
---|---|
region_id | |
9 | POLYGON ((0.02206 51.09168, 0.08839 51.07032, ... |
0 | POLYGON ((0.16366 50.62267, 0.14872 50.7122, 0... |
7 | POLYGON ((-5.47797 52.04071, -5.32917 52.01829... |
13 | POLYGON ((1.6716 51.9353, 1.63184 52.02361, 1.... |
6 | POLYGON ((-13.9945 57.65907, -13.98267 57.6735... |
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 |
... | ... | ... | ... |
74683 | 47.874740 | 7.028350 | FR |
74684 | 48.086210 | 7.278890 | FR |
74685 | 42.044830 | -6.564990 | ES |
74686 | 42.067773 | 13.897856 | IT |
74687 | 46.515031 | 6.559780 | CH |
70494 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%|▎ | 460/13617 [00:00<00:02, 4591.42it/s]
Generating polygons: 7%|▋ | 920/13617 [00:00<00:02, 4564.71it/s]
Generating polygons: 10%|█ | 1377/13617 [00:00<00:02, 4453.58it/s]
Generating polygons: 14%|█▍ | 1874/13617 [00:00<00:02, 4653.10it/s]
Generating polygons: 17%|█▋ | 2358/13617 [00:00<00:02, 4719.40it/s]
Generating polygons: 21%|██ | 2831/13617 [00:00<00:02, 4692.03it/s]
Generating polygons: 24%|██▍ | 3324/13617 [00:00<00:02, 4768.47it/s]
Generating polygons: 28%|██▊ | 3811/13617 [00:00<00:02, 4800.52it/s]
Generating polygons: 32%|███▏ | 4304/13617 [00:00<00:01, 4838.72it/s]
Generating polygons: 35%|███▌ | 4789/13617 [00:01<00:01, 4808.95it/s]
Generating polygons: 39%|███▊ | 5271/13617 [00:01<00:01, 4805.36it/s]
Generating polygons: 42%|████▏ | 5764/13617 [00:01<00:01, 4841.26it/s]
Generating polygons: 46%|████▌ | 6249/13617 [00:01<00:01, 4836.03it/s]
Generating polygons: 49%|████▉ | 6733/13617 [00:01<00:01, 4466.06it/s]
Generating polygons: 53%|█████▎ | 7186/13617 [00:01<00:01, 4261.90it/s]
Generating polygons: 56%|█████▌ | 7618/13617 [00:01<00:01, 4079.80it/s]
Generating polygons: 60%|█████▉ | 8112/13617 [00:01<00:01, 4313.88it/s]
Generating polygons: 63%|██████▎ | 8583/13617 [00:01<00:01, 4424.27it/s]
Generating polygons: 66%|██████▋ | 9050/13617 [00:01<00:01, 4492.09it/s]
Generating polygons: 70%|███████ | 9544/13617 [00:02<00:00, 4621.47it/s]
Generating polygons: 74%|███████▍ | 10043/13617 [00:02<00:00, 4727.24it/s]
Generating polygons: 77%|███████▋ | 10542/13617 [00:02<00:00, 4804.11it/s]
Generating polygons: 81%|████████ | 11041/13617 [00:02<00:00, 4857.25it/s]
Generating polygons: 85%|████████▍ | 11533/13617 [00:02<00:00, 4875.63it/s]
Generating polygons: 88%|████████▊ | 12028/13617 [00:02<00:00, 4896.34it/s]
Generating polygons: 92%|█████████▏| 12519/13617 [00:02<00:00, 4649.95it/s]
Generating polygons: 95%|█████████▌| 12987/13617 [00:02<00:00, 4643.69it/s]
Generating polygons: 99%|█████████▉| 13454/13617 [00:02<00:00, 4646.30it/s]
Generating polygons: 100%|██████████| 13617/13617 [00:02<00:00, 4607.52it/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