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:>
../_images/unit_5_1.png
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:>
../_images/unit_6_1.png
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:>
../_images/unit_15_1.png
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:>
../_images/unit_19_1.png
union = buildings.buffer(100).unary_union
ind = roads.sindex.query(union, predicate='contains')
roads[~roads.index.isin(ind)].plot()
<AxesSubplot:>
../_images/unit_22_1.png
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:>
../_images/unit_27_1.png
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:>
../_images/unit_33_1.png
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:>
../_images/unit_40_1.png
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')