Exploration of tessellation¶
In search of an optimal spatial unit, this notebook explores several options of tessellation.
import os
import geopandas as gpd
from sqlalchemy import create_engine
user = os.environ.get('DB_USER')
pwd = os.environ.get('DB_PWD')
host = os.environ.get('DB_HOST')
port = os.environ.get('DB_PORT')
db_connection_url = f"postgres+psycopg2://{user}:{pwd}@{host}:{port}/built_env"
engine = create_engine(db_connection_url)
engine.begin()
<sqlalchemy.engine.base.Engine._trans_ctx at 0x7f81ec0569d0>
x, y = 352125.32, 492802.86 # coordinates in epsg 27700
buffer = 5000 # radius in [m]
sql = f'SELECT * FROM openroads_200803_topological WHERE ST_DWithin(geometry, ST_SetSRID(ST_Point({x}, {y}), 27700), {buffer})'
roads = gpd.read_postgis(sql, engine, geom_col='geometry')
roads.plot()
<AxesSubplot:>
sql = f'SELECT * FROM openmap_buildings_200814 WHERE ST_DWithin(geometry, ST_SetSRID(ST_Point({x}, {y}), 27700), {buffer})'
buildings = gpd.read_postgis(sql, engine, geom_col='geometry')
buildings.plot()
<AxesSubplot:>
buildings.shape
(7834, 3)
from shapely.geometry import Point
limit = Point(x, y).buffer(buffer)
import momepy as mm
buildings
| gml_id | featureCode | geometry | |
|---|---|---|---|
| 0 | idD589F9B4-E82E-47EC-956C-D28DD3D2F706 | 15014 | POLYGON ((350715.880 497538.310, 350725.390 49... |
| 1 | id1A1D8810-34E8-4BCD-BD81-C1FAB20C0944 | 15014 | POLYGON ((351065.450 497619.740, 351051.940 49... |
| 2 | idC9145C5D-BF42-4CCA-8CAD-601FAD44E3C4 | 15014 | POLYGON ((351062.650 497630.330, 351055.710 49... |
| 3 | idA0A87301-8873-41D3-B9EE-A5163061DDC1 | 15014 | POLYGON ((351085.090 497646.280, 351086.240 49... |
| 4 | id540C742C-8BA8-4F16-8760-3EDB6FFB9F6B | 15014 | POLYGON ((351096.110 497573.180, 351096.080 49... |
| ... | ... | ... | ... |
| 7829 | id0620DFB7-966B-45F2-B2F6-86F9CBBAD225 | 15014 | POLYGON ((350879.770 490856.520, 350891.250 49... |
| 7830 | idA38FAB2E-CAAC-4F21-85B3-D1F3D7BB4FAF | 15014 | POLYGON ((350880.310 490761.290, 350876.970 49... |
| 7831 | idC9A8E7D1-179D-4447-B113-23385DF6B2D5 | 15014 | POLYGON ((350881.340 490881.890, 350887.600 49... |
| 7832 | idA55AB039-E66B-4097-AD1C-CFE92EF2E87B | 15014 | POLYGON ((347493.590 491128.970, 347493.870 49... |
| 7833 | idD95BE56F-3F70-464C-8A04-F623541B514D | 15014 | POLYGON ((347554.800 491127.770, 347557.100 49... |
7834 rows × 3 columns
roads['uID'] = range(len(roads))
tes_rd = mm.Tessellation(roads, 'uID', limit, segment=10)
Inward offset...
Discretization...
33%|███▎ | 758/2271 [00:00<00:00, 7524.43it/s]Generating input point array...
100%|██████████| 2271/2271 [00:00<00:00, 6093.07it/s]
Generating Voronoi diagram...
Vertices to Polygons: 4%|▍ | 1511/36736 [00:00<00:02, 15086.56it/s]Generating GeoDataFrame...
Vertices to Polygons: 100%|██████████| 36736/36736 [00:01<00:00, 26861.00it/s]
Dissolving Voronoi polygons...
Preparing limit for edge resolving...
Building R-tree...
Identifying edge cells...
Cutting...
roads_tessellation = tes_rd.tessellation
ax = roads_tessellation.plot(figsize=(12, 12))
roads.plot(ax=ax, color='k', linewidth=.2)
<AxesSubplot:>
buildings['uID'] = range(len(buildings))
tes_blg = mm.Tessellation(buildings, 'uID', limit, segment=2, shrink=0)
Inward offset...
Discretization...
7%|▋ | 540/7834 [00:00<00:01, 5394.51it/s]Generating input point array...
100%|██████████| 7834/7834 [00:02<00:00, 3899.12it/s]
Generating Voronoi diagram...
Generating GeoDataFrame...
Vertices to Polygons: 100%|██████████| 267225/267225 [00:08<00:00, 32991.76it/s]
Dissolving Voronoi polygons...
Preparing limit for edge resolving...
Building R-tree...
Identifying edge cells...
Cutting...
blg_tessellation = tes_blg.tessellation
ax = blg_tessellation.plot(figsize=(16, 16))
buildings.plot(ax=ax, color='k')
<AxesSubplot:>
union = buildings.buffer(100).unary_union
ind = roads.sindex.query(union, predicate='contains')
roads[~roads.index.isin(ind)].plot()
<AxesSubplot:>
country_roads = roads[~roads.index.isin(ind)].copy()
country_roads.uID = country_roads.uID + 10000
mixed = country_roads[['uID', 'geometry']].append(buildings[['uID', 'geometry']])
tes_mix = mm.Tessellation(mixed, 'uID', limit, segment=2, shrink=0)
Inward offset...
Discretization...
1%| | 96/8192 [00:00<00:08, 926.40it/s]Generating input point array...
100%|██████████| 8192/8192 [00:02<00:00, 3284.89it/s]
Generating Voronoi diagram...
Generating GeoDataFrame...
Vertices to Polygons: 100%|██████████| 340897/340897 [00:12<00:00, 27511.88it/s]
Dissolving Voronoi polygons...
100%|██████████| 7/7 [00:00<00:00, 568.84it/s]Preparing limit for edge resolving...
Building R-tree...
Identifying edge cells...
Cutting...
mixed_tessellation = tes_mix.tessellation
ax = mixed_tessellation.plot(figsize=(16, 16))
mixed.plot(ax=ax, color='k')
<AxesSubplot:>
total = buildings.total_bounds
import numpy as np
pts = []
for x in np.linspace(total[0], total[2], 100):
for y in np.linspace(total[1], total[3], 100):
pts.append(Point(x, y))
pts = gpd.GeoSeries(pts)
grid = pts.drop(pts.sindex.query(union, predicate='intersects'))
ax = grid.plot(markersize=.1, figsize=(12, 12))
buildings.plot(ax=ax)
<AxesSubplot:>
grid = gpd.GeoDataFrame(geometry=grid, crs=buildings.crs)
grid['uID'] = range(len(grid))
grid['uID'] = grid['uID'] + 100000
# fix for momepy
grid.geometry = grid.buffer(.1)
buildings_grid = buildings[['uID', 'geometry']].append(grid).reset_index(drop=True)
buildings_grid
| uID | geometry | |
|---|---|---|
| 0 | 0 | POLYGON ((350715.880 497538.310, 350725.390 49... |
| 1 | 1 | POLYGON ((351065.450 497619.740, 351051.940 49... |
| 2 | 2 | POLYGON ((351062.650 497630.330, 351055.710 49... |
| 3 | 3 | POLYGON ((351085.090 497646.280, 351086.240 49... |
| 4 | 4 | POLYGON ((351096.110 497573.180, 351096.080 49... |
| ... | ... | ... |
| 15185 | 107351 | POLYGON ((357117.710 497363.963, 357117.710 49... |
| 15186 | 107352 | POLYGON ((357117.710 497464.178, 357117.710 49... |
| 15187 | 107353 | POLYGON ((357117.710 497564.392, 357117.710 49... |
| 15188 | 107354 | POLYGON ((357117.710 497664.606, 357117.710 49... |
| 15189 | 107355 | POLYGON ((357117.710 497764.820, 357117.710 49... |
15190 rows × 2 columns
tes_grid = mm.Tessellation(gpd.clip(buildings_grid, limit), 'uID', limit, segment=2, shrink=0)
Inward offset...
Discretization...
4%|▍ | 526/12953 [00:00<00:02, 5258.59it/s]Generating input point array...
100%|██████████| 12953/12953 [00:03<00:00, 3282.79it/s]
Generating Voronoi diagram...
Generating GeoDataFrame...
Vertices to Polygons: 100%|██████████| 594931/594931 [00:16<00:00, 35040.58it/s]
Dissolving Voronoi polygons...
Preparing limit for edge resolving...
Building R-tree...
Identifying edge cells...
Cutting...
grid_tessellation = tes_grid.tessellation
ax = grid_tessellation.plot(figsize=(20, 20), edgecolor='w', linewidth=.2)
buildings_grid.plot(ax=ax, color='k')
<AxesSubplot:>
buildings_grid.to_file('units.gpkg', layer='buildings_grid', driver='GPKG')
grid_tessellation.to_file('units.gpkg', layer='grid_tessellation', driver='GPKG')
buildings.to_file('units.gpkg', layer='buildings', driver='GPKG')
roads.to_file('units.gpkg', layer='roads', driver='GPKG')
mixed.to_file('units.gpkg', layer='mixed', driver='GPKG')
blg_tessellation.to_file('units.gpkg', layer='blg_tessellation', driver='GPKG')
roads_tessellation.to_file('units.gpkg', layer='roads_tessellation', driver='GPKG')
mixed_tessellation.to_file('units.gpkg', layer='mixed_tessellation', driver='GPKG')