ENH: Added Cellular Automata-based Enclosed Tessellation #686

yu-ta-sato · 2025-02-25T03:09:21Z

Based on the discussion #684, I added Cellular Automata-based Voronoi diagram for enclosed_tessellation, remaining the API the same as before.

Here's a sample code:

from momepy import enclosures, enclosed_tessellation
import geopandas as gpd
import matplotlib.pyplot as plt

blg_polygons = [
    shapely.geometry.Polygon([(15, 32), (35, 32), (35, 38), (15, 38)]),
    shapely.geometry.Polygon([(15, 22), (35, 22), (35, 28), (15, 28)]),
    shapely.geometry.Polygon([(15, 92), (35, 92), (35, 98), (15, 98)]),
    shapely.geometry.Polygon([(15, 82), (35, 82), (35, 88), (15, 88)]),
    shapely.geometry.Polygon([(45, 62), (65, 62), (65, 68), (45, 68)]),
    shapely.geometry.Polygon([(45, 52), (65, 52), (65, 58), (45, 58)]),
]
buildings_gdf = GeoDataFrame(
    {
        "building_id": list(range(1, len(blg_polygons) + 1)),
        "geometry": blg_polygons,
    },
    crs="EPSG:3857",
)

# Define sample barrier geometries.
barrier_geoms = [
    shapely.geometry.LineString([(0, 0), (80, 0)]),
    shapely.geometry.LineString([(80, 0), (80, 120)]),
    shapely.geometry.LineString([(80, 120), (0, 120)]),
    shapely.geometry.LineString([(0, 120), (0, 0)]),
    shapely.geometry.LineString([(40, 0), (40, 110)]),
    shapely.geometry.LineString([(10, 30), (40, 30)]),
    shapely.geometry.LineString([(10, 90), (40, 90)]),
    shapely.geometry.LineString([(40, 60), (70, 60)]),
]

barrier_gdf = GeoDataFrame(
    {
        "name": [
            "Bottom Edge",
            "Right Edge",
            "Top Edge",
            "Left Edge",
            "Main Vertical",
            "Left Cul-de-Sac (Bottom)",
            "Left Cul-de-Sac (Top)",
            "Right Cul-de-Sac (Middle)",
        ],
        "geometry": barrier_geoms,
    },
    crs="EPSG:3857",
)

There're three patterns of enclosed tessellation with comparison:

# Generate enclosure
enclosure = enclosures(barrier_gdf)

# Original enclosed_tessellation function
enclosed_tess_original = enclosed_tessellation(
    buildings_gdf,
    enclosures=enclosure,
)

fig, axes = plt.subplots(1, 3, figsize=(15, 5))

enclosure.plot(ax=axes[0], alpha=0.1)
buildings_gdf.plot(ax=axes[0], color="red")
barrier_gdf.plot(ax=axes[0], color="blue")
enclosed_tess_original.plot(ax=axes[0], alpha=0.3, color="red", edgecolor="black", linewidth=0.5)
axes[0].set_title("Original")

# New enclosed_tessellation function (without inner barriers)
enclosed_tess_no_inner = enclosed_tessellation(
    buildings_gdf,
    enclosures=enclosure,
    use_ca=True,
    cell_size=1.0,
    neighbor_mode="moore",
    fulfill=True)

enclosure.plot(ax=axes[1], alpha=0.1)
buildings_gdf.plot(ax=axes[1], color="red")
barrier_gdf.plot(ax=axes[1], color="blue")
enclosed_tess_no_inner.plot(ax=axes[1], alpha=0.3, color="red", edgecolor="black", linewidth=0.5)
axes[1].set_title("CA, no inner barriers")

# New enclosed_tessellation function (with inner barriers)
enclosed_tess_with_inner = enclosed_tessellation(
    buildings_gdf,
    enclosures=enclosure,
    use_ca=True,
    cell_size=1.0,
    neighbor_mode="moore",
    fulfill=True,
    barriers_for_inner=barrier_gdf)

enclosure.plot(ax=axes[2], alpha=0.1)
buildings_gdf.plot(ax=axes[2], color="red")
barrier_gdf.plot(ax=axes[2], color="blue")
enclosed_tess_with_inner.plot(ax=axes[2], alpha=0.3, color="red", edgecolor="black", linewidth=0.5)
axes[2].set_title("CA, with inner barriers")

plt.tight_layout()
plt.show()

In the Liverpool City Region by the Overture Maps, the output of the new algorithm with inner barriers (cell size: 1m)are like this:

If it's good to go, let me update the docs.

codecov · 2025-02-25T06:58:28Z

Codecov Report

Attention: Patch coverage is 15.86207% with 122 lines in your changes missing coverage. Please review.

Project coverage is 96.5%. Comparing base (4037c70) to head (d754ce4).
Report is 103 commits behind head on main.

Files with missing lines	Patch %	Lines
momepy/functional/_elements.py	15.9%	122 Missing ⚠️

Additional details and impacted files

@@           Coverage Diff           @@
##            main    #686     +/-   ##
=======================================
- Coverage   97.4%   96.5%   -0.8%     
=======================================
  Files         26      40     +14     
  Lines       4328    7350   +3022     
=======================================
+ Hits        4214    7094   +2880     
- Misses       114     256    +142

Files with missing lines	Coverage Δ
momepy/functional/_elements.py	`53.4% <15.9%> (ø)`

... and 1 file with indirect coverage changes

🚀 New features to boost your workflow:

❄ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

jGaboardi · 2025-02-25T14:19:31Z

@yu-ta-sato This is really cool! Will you be able to add appropriate testing for this new functionality?

yu-ta-sato · 2025-02-25T15:23:46Z

@jGaboardi Yes of course! I wanted to be on the same page with you guys in advance, especially the direction of its implementation. If needed, I can separate the function from enclosed_tessellation. Will let @martinfleis have a look at the scripts as well.

martinfleis · 2025-02-25T16:03:59Z

Hold on with testing, I'll have some requests but need more time to properly look at it.

martinfleis · 2025-03-04T13:02:45Z

Hey, enclosed tessellation is supposed to be continuous coverage by definition. As illustrated below, that is not the case at the moment. Also, I would very much prefer to ensure that the lines are as straight as possible, without this artifact of the grid we see here. We can either do that in conversion of the grid to polygons or afterwards using coverage simplification. We will need to touch the API as well but let's focus on ensuring this behaves the way we want first.

yu-ta-sato · 2025-03-04T19:34:29Z

We can either do that in conversion of the grid to polygons or afterwards using coverage simplification.

The cellular-automata grids are already polygonised in the returning outputs (by _get_cell_polygon in _voronoi_by_ca). For the latter, Ramer–Douglas–Peucker algorithm might be applied for the smoothing?

martinfleis · 2025-03-04T19:49:35Z

For the latter, Ramer–Douglas–Peucker algorithm might be applied for the smoothing?

No. We need to ensure that whatever comes from the grid is continuous coverage and then use upcoming coverage_simplify which ensures that the simplification does not break topological relations between neighbouring cells.

yu-ta-sato · 2025-03-04T20:07:10Z

Did you enable fulfill=True? I prepared for the parameter so that the cells of boundaries are assigned to the closest tessellations to ensure the coverage (_assign_adjacent_seed_cells). Hope this is what you mean by the continuous coverage

martinfleis · 2025-03-04T21:51:41Z

Ah, no. I did not. However, this should not be an option, it should be set by default.

yu-ta-sato · 2025-03-12T16:40:16Z

Need to chase the on-going discussion on shapely a bit more, but can we implement coverage_simplify on the polygons once it gets stable in its v2.1? Meanwhile I can work on the testing and API design.

martinfleis · 2025-03-12T17:08:41Z

We can already do that. It will just be tested in the dev CI environment only.