Enclosed tessellation (proof of a concept)¶
This notebooks is a proof of a concept of enclosed tessellation, i.e. two-step partitioning of space based on building footprints and boundaries (e.g. street network, railway).
Note: An algorithm to generate enclosed tessellation has been implemented in momepy 0.4.0 within momepy.Tessellation.
Load data¶
import os
import geopandas as gpd
import pandas as pd
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)
Enclosed tessellation toy example¶
For the POC, we’ll use road network, rivers and railwat as barriers and building footprints. All now stored in PostGIS database.
x, y = 352125.32, 492802.86 # coordinates in epsg 27700
buffer = 5000 # radius in [m]
sql = f'SELECT * FROM openroads_200803_topological WHERE ST_Intersects(geometry, ST_Buffer(ST_SetSRID(ST_Point({x}, {y}), 27700), {buffer}))'
roads = gpd.read_postgis(sql, engine, geom_col='geometry')
sql = f'SELECT * FROM openrivers_200909 WHERE ST_Intersects(geometry, ST_Buffer(ST_SetSRID(ST_Point({x}, {y}), 27700), {buffer}))'
rivers = gpd.read_postgis(sql, engine, geom_col='geometry')
sql = f'SELECT * FROM openmap_railwaytrack_200824 WHERE ST_Intersects(geometry, ST_Buffer(ST_SetSRID(ST_Point({x}, {y}), 27700), {buffer}))'
railway = gpd.read_postgis(sql, engine, geom_col='geometry')
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')
ax = roads.plot(figsize=(16, 16))
rivers.plot(ax=ax, color='r')
railway.plot(ax=ax, color='k')
<AxesSubplot:>
Enclosed tessellation is a dask-based parallelised algorithm, so let’s start a client first.
from dask.distributed import Client
client = Client()
client
Client
|
Cluster
|
enclosed_tessellation is loaded from tessellation.py in the same directory.
from tessellation import enclosed_tessellation
from shapely.geometry import Point
To specify an external limit of tesellation, we use the same limit we used for querying the data above.
# assign unique IDs
buildings['uID'] = range(len(buildings))
# get road-based polygons
limit = Point(x, y).buffer(buffer)
# merge barriers
barriers = pd.concat([roads.geometry, rivers.geometry, railway.geometry])
res = enclosed_tessellation(buildings, roads, limit, unique_id="uID")
res.plot(figsize=(12, 12))
<AxesSubplot:>
Let’s zoom in to explore the generated geometry.
ax = res.plot(figsize=(12, 12), edgecolor='w')
buildings.plot(ax=ax, color='w', alpha=.6)
ax.set_xlim(351000, 352000)
ax.set_ylim(491000, 492000)
(491000.0, 492000.0)
client.close()
Real life example - Edinburgh MasterMap and OpenMap¶
mastermap = gpd.read_parquet('../../urbangrammar_samba/OS_MasterMap_Buildings_sample/edinburgh.pq')
from shapely.geometry import box
from shapely.wkb import dumps
bounds = mastermap.total_bounds
limit = box(*bounds)
sql = f"SELECT * FROM openroads_200803_topological WHERE ST_Intersects(geometry, ST_GeomFromText('{limit.wkt}',27700))"
roads = gpd.read_postgis(sql, engine, geom_col='geometry')
sql = f"SELECT * FROM openmap_buildings_200814 WHERE ST_DWithin(geometry, ST_GeomFromText('{limit.wkt}',27700), 0)"
openmap = gpd.read_postgis(sql, engine, geom_col='geometry')
/opt/conda/lib/python3.7/asyncio/base_events.py:626: ResourceWarning: unclosed event loop <_UnixSelectorEventLoop running=False closed=False debug=False>
source=self)
from shapely.ops import polygonize
polygons = polygonize(roads.geometry.append(gpd.GeoSeries([limit.boundary])).unary_union)
enclosures = gpd.array.from_shapely(list(polygons), crs=roads.crs)
OpenMap¶
client = Client()
openmap['uID'] = range(len(openmap))
openmap_tess = enclosed_tessellation(roads, openmap, limit, unique_id="uID", enclosures=enclosures)
openmap_tess.to_parquet('../../urbangrammar_samba/OS_MasterMap_Buildings_sample/edinburgh_om_tess.pq')
/opt/conda/lib/python3.7/site-packages/ipykernel_launcher.py:1: UserWarning: this is an initial implementation of Parquet/Feather file support and associated metadata. This is tracking version 0.1.0 of the metadata specification at https://github.com/geopandas/geo-arrow-spec
This metadata specification does not yet make stability promises. We do not yet recommend using this in a production setting unless you are able to rewrite your Parquet/Feather files.
To further ignore this warning, you can do:
import warnings; warnings.filterwarnings('ignore', message='.*initial implementation of Parquet.*')
"""Entry point for launching an IPython kernel.
MasterMap¶
Sometimes it is better to restart dask client to start with clean state.
client.restart()
distributed.nanny - WARNING - Restarting worker
distributed.nanny - WARNING - Restarting worker
distributed.nanny - WARNING - Restarting worker
distributed.nanny - WARNING - Restarting worker
distributed.nanny - WARNING - Restarting worker
Client
|
Cluster
|
mastermap['uID'] = range(len(mastermap))
mastermap_tess = enclosed_tessellation(roads, mastermap, limit, unique_id="uID", enclosures=enclosures)
mastermap_tess.to_parquet('../../urbangrammar_samba/OS_MasterMap_Buildings_sample/edinburgh_mm_tess.pq')
/opt/conda/lib/python3.7/site-packages/ipykernel_launcher.py:1: UserWarning: this is an initial implementation of Parquet/Feather file support and associated metadata. This is tracking version 0.1.0 of the metadata specification at https://github.com/geopandas/geo-arrow-spec
This metadata specification does not yet make stability promises. We do not yet recommend using this in a production setting unless you are able to rewrite your Parquet/Feather files.
To further ignore this warning, you can do:
import warnings; warnings.filterwarnings('ignore', message='.*initial implementation of Parquet.*')
"""Entry point for launching an IPython kernel.
client.close()
mastermap.to_parquet('../../urbangrammar_samba/OS_MasterMap_Buildings_sample/edinburgh.pq')
/opt/conda/lib/python3.7/site-packages/ipykernel_launcher.py:1: UserWarning: this is an initial implementation of Parquet/Feather file support and associated metadata. This is tracking version 0.1.0 of the metadata specification at https://github.com/geopandas/geo-arrow-spec
This metadata specification does not yet make stability promises. We do not yet recommend using this in a production setting unless you are able to rewrite your Parquet/Feather files.
To further ignore this warning, you can do:
import warnings; warnings.filterwarnings('ignore', message='.*initial implementation of Parquet.*')
"""Entry point for launching an IPython kernel.
openmap.to_parquet('../../urbangrammar_samba/OS_MasterMap_Buildings_sample/edinburgh_om.pq')
/opt/conda/lib/python3.7/site-packages/ipykernel_launcher.py:1: UserWarning: this is an initial implementation of Parquet/Feather file support and associated metadata. This is tracking version 0.1.0 of the metadata specification at https://github.com/geopandas/geo-arrow-spec
This metadata specification does not yet make stability promises. We do not yet recommend using this in a production setting unless you are able to rewrite your Parquet/Feather files.
To further ignore this warning, you can do:
import warnings; warnings.filterwarnings('ignore', message='.*initial implementation of Parquet.*')
"""Entry point for launching an IPython kernel.