Adjacency neighbourhood
In [1]:
Copied!
import geopandas as gpd
import numpy as np
from shapely.geometry import Point
from srai.constants import WGS84_CRS
from srai.neighbourhoods import AdjacencyNeighbourhood
from srai.plotting.folium_wrapper import plot_all_neighbourhood, plot_neighbours, plot_regions
from srai.regionalizers import (
AdministrativeBoundaryRegionalizer,
VoronoiRegionalizer,
geocode_to_region_gdf,
)
import geopandas as gpd
import numpy as np
from shapely.geometry import Point
from srai.constants import WGS84_CRS
from srai.neighbourhoods import AdjacencyNeighbourhood
from srai.plotting.folium_wrapper import plot_all_neighbourhood, plot_neighbours, plot_regions
from srai.regionalizers import (
AdministrativeBoundaryRegionalizer,
VoronoiRegionalizer,
geocode_to_region_gdf,
)
Adjacency Neighbourhood¶
It can generate neighbourhoods for all geodataframes with touching geometries.
Real boundaries example - Italy¶
In [2]:
Copied!
it_gdf = geocode_to_region_gdf(query=["R365331"], by_osmid=True)
plot_regions(it_gdf)
it_gdf = geocode_to_region_gdf(query=["R365331"], by_osmid=True)
plot_regions(it_gdf)
Out[2]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [3]:
Copied!
regionalizer = AdministrativeBoundaryRegionalizer(admin_level=4)
it_regions_gdf = regionalizer.transform(it_gdf)
regionalizer = AdministrativeBoundaryRegionalizer(admin_level=4)
it_regions_gdf = regionalizer.transform(it_gdf)
Loading boundaries: 0: 0%| | 0/6 [00:00<?, ?it/s]
Loading boundaries: 1: 0%| | 0/6 [00:05<?, ?it/s]
Loading boundaries: 2: 0%| | 0/6 [00:06<?, ?it/s]
Loading boundaries: 3: 0%| | 0/6 [00:08<?, ?it/s]
Loading boundaries: 4: 0%| | 0/6 [00:09<?, ?it/s]
Loading boundaries: 5: 0%| | 0/6 [00:11<?, ?it/s]
Loading boundaries: 6: 0%| | 0/6 [00:13<?, ?it/s]
Loading boundaries: 7: 0%| | 0/6 [00:14<?, ?it/s]
Loading boundaries: 8: 0%| | 0/6 [00:16<?, ?it/s]
Loading boundaries: 9: 0%| | 0/6 [00:17<?, ?it/s]
Loading boundaries: 10: 0%| | 0/6 [00:18<?, ?it/s]
Loading boundaries: 11: 0%| | 0/6 [00:20<?, ?it/s]
Loading boundaries: 12: 0%| | 0/6 [00:22<?, ?it/s]
Loading boundaries: 13: 0%| | 0/6 [00:23<?, ?it/s]
Loading boundaries: 14: 0%| | 0/6 [00:25<?, ?it/s]
Loading boundaries: 15: 0%| | 0/6 [00:27<?, ?it/s]
Loading boundaries: 16: 0%| | 0/6 [00:28<?, ?it/s]
Loading boundaries: 17: 0%| | 0/6 [00:30<?, ?it/s]
Loading boundaries: 18: 0%| | 0/6 [00:31<?, ?it/s]
Loading boundaries: 19: 0%| | 0/6 [00:33<?, ?it/s]
Loading boundaries: 20: 0%| | 0/6 [00:35<?, ?it/s]
Loading boundaries: 21: 0%| | 0/6 [00:36<?, ?it/s]
Loading boundaries: 22: 0%| | 0/6 [00:38<?, ?it/s]
Loading boundaries: 23: 0%| | 0/6 [00:39<?, ?it/s]
Loading boundaries: 24: 0%| | 0/6 [00:41<?, ?it/s]
Loading boundaries: 25: 0%| | 0/6 [00:43<?, ?it/s]
Loading boundaries: 26: 0%| | 0/6 [00:45<?, ?it/s]
Loading boundaries: 27: 0%| | 0/6 [00:47<?, ?it/s]
Loading boundaries: 28: 0%| | 0/6 [00:49<?, ?it/s]
Loading boundaries: 29: 0%| | 0/6 [00:50<?, ?it/s]
Loading boundaries: 30: 0%| | 0/6 [00:52<?, ?it/s]
Loading boundaries: 31: 0%| | 0/6 [00:53<?, ?it/s]
Loading boundaries: 32: 0%| | 0/6 [00:55<?, ?it/s]
Loading boundaries: 33: 0%| | 0/6 [00:56<?, ?it/s]
Loading boundaries: 34: 0%| | 0/6 [00:58<?, ?it/s]
Loading boundaries: 35: 0%| | 0/6 [01:00<?, ?it/s]
Loading boundaries: 36: 0%| | 0/6 [01:01<?, ?it/s]
Loading boundaries: 37: 0%| | 0/6 [01:03<?, ?it/s]
Loading boundaries: 38: 0%| | 0/6 [01:04<?, ?it/s]
Loading boundaries: 39: 0%| | 0/6 [01:06<?, ?it/s]
Loading boundaries: 40: 0%| | 0/6 [01:07<?, ?it/s]
Loading boundaries: 41: 0%| | 0/6 [01:09<?, ?it/s]
Loading boundaries: 42: 0%| | 0/6 [01:10<?, ?it/s]
Loading boundaries: 43: 0%| | 0/6 [01:12<?, ?it/s]
Loading boundaries: 44: 0%| | 0/6 [01:14<?, ?it/s]
Loading boundaries: 45: 0%| | 0/6 [01:15<?, ?it/s]
Loading boundaries: 46: 0%| | 0/6 [01:17<?, ?it/s]
Loading boundaries: 47: 0%| | 0/6 [01:18<?, ?it/s]
Loading boundaries: 48: 0%| | 0/6 [01:20<?, ?it/s]
Loading boundaries: 49: 0%| | 0/6 [01:21<?, ?it/s]
Loading boundaries: 50: 0%| | 0/6 [01:23<?, ?it/s]
Loading boundaries: 51: 0%| | 0/6 [01:25<?, ?it/s]
Loading boundaries: 52: 0%| | 0/6 [01:26<?, ?it/s]
Loading boundaries: 53: 0%| | 0/6 [01:28<?, ?it/s]
Loading boundaries: 54: 0%| | 0/6 [01:29<?, ?it/s]
Loading boundaries: 55: 0%| | 0/6 [01:31<?, ?it/s]
Loading boundaries: 56: 0%| | 0/6 [01:33<?, ?it/s]
Loading boundaries: 57: 0%| | 0/6 [01:34<?, ?it/s]
Loading boundaries: 58: 0%| | 0/6 [01:36<?, ?it/s]
Loading boundaries: 59: 0%| | 0/6 [01:37<?, ?it/s]
Loading boundaries: 60: 0%| | 0/6 [01:39<?, ?it/s]
Loading boundaries: 61: 0%| | 0/6 [01:40<?, ?it/s]
Loading boundaries: 62: 0%| | 0/6 [01:42<?, ?it/s]
Loading boundaries: 63: 0%| | 0/6 [01:44<?, ?it/s]
Loading boundaries: 64: 0%| | 0/6 [01:46<?, ?it/s]
Loading boundaries: 65: 0%| | 0/6 [01:48<?, ?it/s]
Loading boundaries: 66: 0%| | 0/6 [01:49<?, ?it/s]
Loading boundaries: 67: 0%| | 0/6 [01:51<?, ?it/s]
Loading boundaries: 68: 0%| | 0/6 [01:54<?, ?it/s]
Loading boundaries: 68: 17%|█▋ | 1/6 [01:54<09:30, 114.12s/it]
Loading boundaries: 68: 33%|███▎ | 2/6 [01:57<03:15, 48.89s/it]
Loading boundaries: 68: 50%|█████ | 3/6 [01:57<01:20, 26.70s/it]
Loading boundaries: 68: 67%|██████▋ | 4/6 [01:58<00:33, 16.56s/it]
Loading boundaries: 68: 83%|████████▎ | 5/6 [01:59<00:10, 10.84s/it]
Loading boundaries: 68: 100%|██████████| 6/6 [02:00<00:00, 7.51s/it]
Loading boundaries: 68: 100%|██████████| 6/6 [02:00<00:00, 20.07s/it]
/opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/srai/regionalizers/administrative_boundary_regionalizer.py:169: FutureWarning: `unary_union` returned None due to all-None GeoSeries. In future, `unary_union` will return 'GEOMETRYCOLLECTION EMPTY' instead. ].geometry.unary_union
In [4]:
Copied!
plot_regions(it_regions_gdf)
plot_regions(it_regions_gdf)
Out[4]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [5]:
Copied!
neighbourhood = AdjacencyNeighbourhood(it_regions_gdf)
neighbourhood = AdjacencyNeighbourhood(it_regions_gdf)
Nearest neighbours¶
In [6]:
Copied!
region_id = "Lazio"
neighbours = neighbourhood.get_neighbours(region_id)
neighbours
region_id = "Lazio"
neighbours = neighbourhood.get_neighbours(region_id)
neighbours
Out[6]:
{'Abruzzo', 'Campania', 'Marche', 'Molise', 'Tuscany', 'Umbria'}
In [7]:
Copied!
plot_neighbours(it_regions_gdf, region_id, neighbours)
plot_neighbours(it_regions_gdf, region_id, neighbours)
Out[7]:
Make this Notebook Trusted to load map: File -> Trust Notebook
Neighbours at a distance¶
In [8]:
Copied!
region_id = "Basilicata"
neighbours = neighbourhood.get_neighbours_at_distance(region_id, 2)
neighbours
region_id = "Basilicata"
neighbours = neighbourhood.get_neighbours_at_distance(region_id, 2)
neighbours
Out[8]:
{'Lazio', 'Molise'}
In [9]:
Copied!
plot_neighbours(it_regions_gdf, region_id, neighbours)
plot_neighbours(it_regions_gdf, region_id, neighbours)
Out[9]:
Make this Notebook Trusted to load map: File -> Trust Notebook
Regions without neighbours¶
In [10]:
Copied!
region_id = "Sardinia"
neighbours = neighbourhood.get_neighbours(region_id)
neighbours
region_id = "Sardinia"
neighbours = neighbourhood.get_neighbours(region_id)
neighbours
Out[10]:
set()
In [11]:
Copied!
plot_neighbours(it_regions_gdf, region_id, neighbours)
plot_neighbours(it_regions_gdf, region_id, neighbours)
Out[11]:
Make this Notebook Trusted to load map: File -> Trust Notebook
Plotting all neighbourhood¶
In [12]:
Copied!
region_id = "Campania"
plot_all_neighbourhood(it_regions_gdf, region_id, neighbourhood)
region_id = "Campania"
plot_all_neighbourhood(it_regions_gdf, region_id, neighbourhood)
Out[12]:
Make this Notebook Trusted to load map: File -> Trust Notebook
Voronoi example - Australia¶
In [13]:
Copied!
au_gdf = geocode_to_region_gdf(query=["R80500"], by_osmid=True)
plot_regions(au_gdf)
au_gdf = geocode_to_region_gdf(query=["R80500"], by_osmid=True)
plot_regions(au_gdf)
Out[13]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [14]:
Copied!
def generate_random_points(shape, n_points=500):
"""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=500):
"""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 [15]:
Copied!
pts = generate_random_points(au_gdf.geometry[0])
au_seeds_gdf = gpd.GeoDataFrame(
{"geometry": pts},
index=list(range(len(pts))),
crs=WGS84_CRS,
)
pts = generate_random_points(au_gdf.geometry[0])
au_seeds_gdf = gpd.GeoDataFrame(
{"geometry": pts},
index=list(range(len(pts))),
crs=WGS84_CRS,
)
/opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/geopandas/geoseries.py:641: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` val = getattr(super(), mtd)(*args, **kwargs)
In [16]:
Copied!
vr = VoronoiRegionalizer(seeds=au_seeds_gdf)
au_result_gdf = vr.transform(gdf=au_gdf)
vr = VoronoiRegionalizer(seeds=au_seeds_gdf)
au_result_gdf = vr.transform(gdf=au_gdf)
Generating spherical polygons: 0%| | 0/91 [00:00<?, ?it/s]
Generating spherical polygons: 12%|█▏ | 11/91 [00:00<00:00, 95.30it/s]
Generating spherical polygons: 26%|██▋ | 24/91 [00:00<00:00, 81.95it/s]
Generating spherical polygons: 40%|███▉ | 36/91 [00:00<00:00, 80.82it/s]
Generating spherical polygons: 57%|█████▋ | 52/91 [00:00<00:00, 77.20it/s]
Generating spherical polygons: 66%|██████▌ | 60/91 [00:01<00:00, 46.46it/s]
Generating spherical polygons: 73%|███████▎ | 66/91 [00:01<00:00, 40.85it/s]
Generating spherical polygons: 78%|███████▊ | 71/91 [00:01<00:00, 42.01it/s]
Generating spherical polygons: 86%|████████▌ | 78/91 [00:01<00:00, 43.42it/s]
Generating spherical polygons: 93%|█████████▎| 85/91 [00:01<00:00, 46.35it/s]
Generating spherical polygons: 100%|██████████| 91/91 [00:01<00:00, 55.66it/s]
Interpolating edges: 0%| | 0/344 [00:00<?, ?it/s]
Interpolating edges: 0%| | 1/344 [00:00<00:52, 6.48it/s]
Interpolating edges: 46%|████▌ | 157/344 [00:00<00:00, 752.58it/s]
Interpolating edges: 91%|█████████ | 312/344 [00:00<00:00, 1075.17it/s]
Interpolating edges: 100%|██████████| 344/344 [00:00<00:00, 919.20it/s]
Generating polygons: 0%| | 0/91 [00:00<?, ?it/s]
Generating polygons: 69%|██████▉ | 63/91 [00:00<00:00, 543.86it/s]
Generating polygons: 100%|██████████| 91/91 [00:00<00:00, 569.97it/s]
In [17]:
Copied!
folium_map = plot_regions(au_result_gdf, tiles_style="CartoDB positron")
au_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(au_result_gdf, tiles_style="CartoDB positron")
au_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
In [18]:
Copied!
neighbourhood = AdjacencyNeighbourhood(regions_gdf=au_result_gdf)
neighbourhood = AdjacencyNeighbourhood(regions_gdf=au_result_gdf)
Nearest neighbours¶
In [19]:
Copied!
region_id = 0
neighbours = neighbourhood.get_neighbours(region_id)
neighbours
region_id = 0
neighbours = neighbourhood.get_neighbours(region_id)
neighbours
Out[19]:
{21, 22, 36, 41, 63, 71, 73}
In [20]:
Copied!
plot_neighbours(au_result_gdf, region_id, neighbours)
plot_neighbours(au_result_gdf, region_id, neighbours)
Out[20]:
Make this Notebook Trusted to load map: File -> Trust Notebook
Neighbours at a distance¶
In [21]:
Copied!
region_id = 0
neighbours = neighbourhood.get_neighbours_at_distance(region_id, 3)
neighbours
region_id = 0
neighbours = neighbourhood.get_neighbours_at_distance(region_id, 3)
neighbours
Out[21]:
{5, 6, 9, 11, 13, 14, 15, 18, 31, 34, 47, 49, 52, 55, 66, 78, 80, 88, 90}
In [22]:
Copied!
plot_neighbours(au_result_gdf, region_id, neighbours)
plot_neighbours(au_result_gdf, region_id, neighbours)
Out[22]:
Make this Notebook Trusted to load map: File -> Trust Notebook
Plotting all neighbourhood¶
In [23]:
Copied!
region_id = 0
plot_all_neighbourhood(au_result_gdf, region_id, neighbourhood)
region_id = 0
plot_all_neighbourhood(au_result_gdf, region_id, neighbourhood)
Out[23]:
Make this Notebook Trusted to load map: File -> Trust Notebook