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:>
../_images/enclosed_tessellation_5_1.png

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

  • Workers: 7
  • Cores: 28
  • Memory: 84.28 GB

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:>
../_images/enclosed_tessellation_14_1.png

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)
../_images/enclosed_tessellation_16_1.png
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

  • Workers: 7
  • Cores: 28
  • Memory: 84.28 GB
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.