In [1]:
import os

os.environ.get('GDS_ENV_VERSION')
Out[1]:
'6.0alpha'

Generate illustrations of tessellation

This notebook contains one function pipeline, which for a given point (lat, lon) generates a sequence of seven images illustrating the process of creation of morphologicla tessellation within 250m buffer. The function is used to generate animations and figures in the blogpost.

In [2]:
import geopandas as gpd
import momepy as mm
import pygeos
import numpy as np
from scipy.spatial import Voronoi
import pandas as pd
from mapclassify import greedy
import contextily as ctx
import matplotlib.pyplot as plt
from palettable.wesanderson import FantasticFox2_5
Matplotlib created a temporary config/cache directory at /tmp/matplotlib-dq2ou8_z because the default path (/home/jovyan/.config/matplotlib) is not a writable directory; it is highly recommended to set the MPLCONFIGDIR environment variable to a writable directory, in particular to speed up the import of Matplotlib and to better support multiprocessing.
In [187]:
def pipeline(lat, lon, path, prefix, dist=250, figsize=(12, 12)):
    point = (lat, lon)
    gdf = ox.geometries.geometries_from_point(point, dist=dist, tags={'building':True})


    gdf_projected = ox.projection.project_gdf(gdf)

    bounds = gdf_projected.total_bounds
    limit = Point(np.mean([bounds[0], bounds[2]]), np.mean([bounds[1], bounds[3]])).buffer(250)
    blg = gpd.clip(gdf_projected, limit).explode()
    bounds = limit.bounds

    # figure 1 - aerial
    fig, ax = plt.subplots(figsize=figsize)
    ax.axis([bounds[0], bounds[2], bounds[1], bounds[3]])
    gpd.GeoSeries([limit.buffer(150).difference(limit)]).plot(ax=ax, color='white')
    ctx.add_basemap(ax, crs=blg.crs, source=ctx.providers.Esri.WorldImagery)
    ax.set_axis_off()
    plt.savefig(path + prefix + "01.png", bbox_inches='tight')
    plt.close()
    print("Figure 1 saved to " + path + prefix + "01.png")

    # figure 2 - overlay
    fig, ax = plt.subplots(figsize=figsize)
    ax.axis([bounds[0], bounds[2], bounds[1], bounds[3]])
    gpd.GeoSeries([limit.buffer(150).difference(limit)]).plot(ax=ax, color='white')
    ctx.add_basemap(ax, crs=blg.crs, source=ctx.providers.Esri.WorldImagery)
    blg.plot(ax=ax, color='#0ea48f', edgecolor='k', alpha=.6)
    ax.set_axis_off()
    plt.savefig(path + prefix + "02.png", bbox_inches='tight')
    plt.close()
    print("Figure 2 saved to " + path + prefix + "02.png")


    # figure 3 - footprints
    fig, ax = plt.subplots(figsize=figsize)
    ax.axis([bounds[0], bounds[2], bounds[1], bounds[3]])
    blg.plot(ax=ax, color='#0ea48f', edgecolor='k').set_axis_off()
    plt.savefig(path + prefix + "03.png", bbox_inches='tight')
    plt.close()
    print("Figure 3 saved to " + path + prefix + "03.png")

    shrinked = blg.buffer(-2)
    shrinked = shrinked[~shrinked.is_empty]

    # figure 4 - shrinked
    fig, ax = plt.subplots(figsize=figsize)
    ax.axis([bounds[0], bounds[2], bounds[1], bounds[3]])
    blg.plot(ax=ax, facecolor='none', linewidth=.5, edgecolor='k')
    shrinked.plot(ax=ax, color='#0ea48f')
    ax.set_axis_off()
    plt.savefig(path + prefix + "04.png", bbox_inches='tight')
    plt.close()
    print("Figure 4 saved to " + path + prefix + "04.png")

    distance = 4
    points = np.empty((0, 2))
    ids = []
    lines = shrinked.boundary.values.data
    lengths = shrinked.length
    for ix, line, length in zip(shrinked.index, lines, lengths):
        if length > distance:
            pts = pygeos.line_interpolate_point(
                line,
                np.linspace(0.1, length - 0.1, num=int((length - 0.1) // distance)),
            )  # .1 offset to keep a gap between two segments
            if len(pts) > 0:
                points = np.append(points, pygeos.get_coordinates(pts), axis=0)
                ids += [ix] * len(pts)

    # figure 5 - points
    fig, ax = plt.subplots(figsize=figsize)
    ax.axis([bounds[0], bounds[2], bounds[1], bounds[3]])
    blg.plot(ax=ax, facecolor='none', linewidth=.5, edgecolor='k')
    gpd.GeoSeries(pygeos.points(points)).plot(ax=ax, markersize=1, color='#0ea48f')
    ax.set_axis_off()
    plt.savefig(path + prefix + "05.png", bbox_inches='tight')
    plt.close()
    print("Figure 5 saved to " + path + prefix + "05.png")

    # add hull to resolve issues with infinity
    # this is just a correction step ensuring the algorithm will work correctly
    stop = points.shape[0]
    series = gpd.GeoSeries(limit)
    hull = series.geometry[[0]].buffer(500)
    line = hull.boundary.values.data[0]
    length = hull.length[0]
    pts = pygeos.line_interpolate_point(
        line,
        np.linspace(0.1, length - 0.1, num=int((length - 0.1) // distance)),
    )  # .1 offset to keep a gap between two segments
    points = np.append(points, pygeos.get_coordinates(pts), axis=0)
    ids += [-1] * len(pts)

    voronoi_diagram = Voronoi(np.array(points))

    vertices = pd.Series(voronoi_diagram.regions).take(voronoi_diagram.point_region)
    polygons = []
    for region in vertices:
        if -1 not in region:
            polygons.append(pygeos.polygons(voronoi_diagram.vertices[region]))
        else:
            polygons.append(None)

    regions_gdf = gpd.GeoDataFrame(
        {'unique_id': ids}, geometry=polygons
    ).dropna()
    regions_gdf = regions_gdf.loc[
        regions_gdf['unique_id'] != -1
    ]  # delete hull-based cells
    voronoi_tessellation = gpd.clip(regions_gdf, limit)


    # figure 6 - voronoi
    fig, ax = plt.subplots(figsize=figsize)
    ax.axis([bounds[0], bounds[2], bounds[1], bounds[3]])
    gpd.GeoSeries(pygeos.points(points[:stop])).plot(ax=ax, markersize=1, zorder=3, color='#0ea48f')
    voronoi_tessellation.plot(ax=ax, facecolor='none', linewidth=.2, edgecolor='gray')
    ax.set_axis_off()
    plt.savefig(path + prefix + "06.png", bbox_inches='tight')
    plt.close()
    print("Figure 6 saved to " + path + prefix + "06.png")

    # figure 7 - tessellation
    fig, ax = plt.subplots(figsize=figsize)
    ax.axis([bounds[0], bounds[2], bounds[1], bounds[3]])
    blg = blg[blg.geom_type == 'Polygon']
    blg = blg.reset_index(drop=True)
    blg['uid'] = range(len(blg))
    tessellation = mm.Tessellation(blg, 'uid', limit, verbose=False).tessellation
    tessellation.plot(greedy(tessellation, strategy='smallest_last'), ax=ax, categorical=True, edgecolor='w', alpha=.6, cmap=FantasticFox2_5.mpl_colormap)
    ax.set_axis_off()
    plt.savefig(path + prefix + "07.png", bbox_inches='tight')
    plt.close()
    print("Figure 7 saved to " + path + prefix + "07.png")
/opt/conda/lib/python3.8/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning: `should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.
  and should_run_async(code)
In [186]:
pipeline(33.9488360, -118.2372975, path='./', prefix='la_', figsize=(15, 15))
/opt/conda/lib/python3.8/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning: `should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.
  and should_run_async(code)
/opt/conda/lib/python3.8/site-packages/osmnx/footprints.py:44: UserWarning: The `footprints` module has been deprecated and will be removed in a future release. Instead, use the `geometries` module's `geometries_from_point` function, passing `tags={'building':True}`.
  warnings.warn(msg)
Figure 1 saved to ./la_01.png
Figure 2 saved to ./la_02.png
Figure 3 saved to ./la_03.png
Figure 4 saved to ./la_04.png
Figure 5 saved to ./la_05.png
Figure 6 saved to ./la_06.png
/opt/conda/lib/python3.8/site-packages/momepy/elements.py:474: UserWarning: Tessellation does not fully match buildings. 5 element(s) collapsed during generation - unique_id: {37, 133, 363, 375, 376}
  warnings.warn(
Figure 7 saved to ./la_07.png
In [188]:
pipeline(41.3907594, 2.1573404, path='./', prefix='bcn_', figsize=(15, 15))
Figure 1 saved to ./bcn_01.png
Figure 2 saved to ./bcn_02.png
Figure 3 saved to ./bcn_03.png
Figure 4 saved to ./bcn_04.png
Figure 5 saved to ./bcn_05.png
Figure 6 saved to ./bcn_06.png
/opt/conda/lib/python3.8/site-packages/momepy/elements.py:474: UserWarning: Tessellation does not fully match buildings. 5 element(s) collapsed during generation - unique_id: {206, 277, 279, 345, 252}
  warnings.warn(
Figure 7 saved to ./bcn_07.png
In [189]:
pipeline(38.995888, -77.135073, path='./', prefix='atl_', figsize=(15, 15))
/opt/conda/lib/python3.8/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning: `should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.
  and should_run_async(code)
Figure 1 saved to ./atl_01.png
Figure 2 saved to ./atl_02.png
Figure 3 saved to ./atl_03.png
Figure 4 saved to ./atl_04.png
Figure 5 saved to ./atl_05.png
Figure 6 saved to ./atl_06.png
/opt/conda/lib/python3.8/site-packages/momepy/elements.py:474: UserWarning: Tessellation does not fully match buildings. 1 element(s) collapsed during generation - unique_id: {21}
  warnings.warn(
Figure 7 saved to ./atl_07.png
In [190]:
pipeline(44.4942640, 11.3473233, path='./', prefix='bol_', figsize=(15, 15))
/opt/conda/lib/python3.8/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning: `should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.
  and should_run_async(code)
Figure 1 saved to ./bol_01.png
Figure 2 saved to ./bol_02.png
Figure 3 saved to ./bol_03.png
Figure 4 saved to ./bol_04.png
Figure 5 saved to ./bol_05.png
Figure 6 saved to ./bol_06.png
/opt/conda/lib/python3.8/site-packages/momepy/elements.py:474: UserWarning: Tessellation does not fully match buildings. 6 element(s) collapsed during generation - unique_id: {352, 808, 684, 685, 716, 175}
  warnings.warn(
Figure 7 saved to ./bol_07.png
In [193]:
pipeline(-15.8038355, -47.8918796, path='./', prefix='bra_', figsize=(15, 15))
/opt/conda/lib/python3.8/site-packages/ipykernel/ipkernel.py:287: DeprecationWarning: `should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.
  and should_run_async(code)
Figure 1 saved to ./bra_01.png
Figure 2 saved to ./bra_02.png
Figure 3 saved to ./bra_03.png
Figure 4 saved to ./bra_04.png
Figure 5 saved to ./bra_05.png
Figure 6 saved to ./bra_06.png
Figure 7 saved to ./bra_07.png
In [ ]: