ENH: Added Cellular Automata-based Enclosed Tessellation

[yu-ta-sato]

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]

Codecov Report

❌ Patch coverage is 98.91892% with 2 lines in your changes missing coverage. Please review.
✅ Project coverage is 98.3%. Comparing base (4037c70) to head (7d534c7).
⚠️ Report is 122 commits behind head on main.

Files with missing lines	Patch %	Lines
momepy/elements.py	98.6%	2 Missing ⚠️

Additional details and impacted files

@@           Coverage Diff           @@
##            main    #686     +/-   ##
=======================================
+ Coverage   97.4%   98.3%   +1.0%     
=======================================
  Files         26      26             
  Lines       4328    4397     +69     
=======================================
+ Hits        4214    4323    +109     
+ Misses       114      74     -40     

Files with missing lines	Coverage Δ
momepy/tests/test_elements.py	100.0% <100.0%> (ø)
momepy/elements.py	98.8% <98.6%> (+1.9%)	⬆️

... and 4 files with indirect coverage changes

🚀 New features to boost your workflow:

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

[jGaboardi]

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

[yu-ta-sato]

@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]

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

[martinfleis]

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]

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]

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]

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]

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

[yu-ta-sato]

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]

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

[martinfleis]

What is the rationale behind this? Why are you ignoring holes when a MultiPolygon is provided? Feels inconsistent.

[yu-ta-sato]

The initial intention was to let the poly be accommodated with inner barriers, but now they're separated. Somehow this part remained, so will remove it.

[martinfleis]

Where does 10 come from here? Can't use hardcoded values here.

[yu-ta-sato]

This amount of buffer was needed to ensure that all the extent are assigned to either of the cell state, which could be arbitrary number (This kind of vacant cells occurred between three tessellations). Will set the number based on the cell size automatically.

[martinfleis]

We should raise this before we start doing anything else. Fail fast.

[martinfleis]

Also, the error should say something like Shapely 2.1.0 or higher is required for tessellation with inner barriers provided.

[martinfleis]

Isn't inner_barriers already a GeoSeries?

[yu-ta-sato]

I checked it, you're right.

[martinfleis]

Hey, I pushed some changes and got this up to speed with main.

The main issue I see now is performance. On the tiny bubenec dataset, this algorithm takes 13.4 s ± 1.23 s per loop (mean ± std. dev. of 7 runs, 1 loop each) while the standard one 172 ms ± 57.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each). That is two orders of magnitude slower. @jGaboardi how do you feel about that?

[jGaboardi]

That is two orders of magnitude slower. @jGaboardi how do you feel about that?

I am not excited about it, but since it is optional functionality maybe OK?

[yu-ta-sato]

@martinfleis @jGaboardi Thanks for the revisit and reminder!

I agree with the concern about the notoriously huge computational order. Although the algorithm of Cell Automata is inebitably complex, I fixed some part of inefficiency as follows:

Replaced per-cell polygon creation (_get_cell_polygon) loops with shapely.box creation and intersection checks.

Before: used np.meshgrid → np.flatten → loop calling _get_cell_polygon for each (x,y) to build a GeoSeries, then used GeoSeries.intersects to test intersections

Now: uses np.meshgrid → np.ravel to get coordinate arrays, creates all cell boxes at once with shapely.box (vectorised), and runs shapely.intersects on the whole array for vectorised intersection tests

In my local environment for bubenec dataset, now it's reduced like this:

Before: 5.33 seconds
Now : 2.52 seconds

Hope it works!

[martinfleis]

Good job!

@yu-ta-sato can you also try to craft some tests for this and expand the user guide on enclosed tessellation? You can add the test to the existing block.

[yu-ta-sato]

@martinfleis yes, I was about to work on them!
Let me come back to you once the implementation is done.

[yu-ta-sato]

Somehow the current code penetrates through the inner barriers from the connected point with the outers.
I'm still investigating the root cause, but it seems it was already happening even after reverting back to before I vectorised the polygon generation.

[martinfleis]

That is probably something I caused in my commits but it is not immediately clear to me where.

[yu-ta-sato]

I identified the root cause! The intersection operation in the generation process of inner barrier was removed throughout the refactoring. Now it's working properly. Let me finalise the test codes based on this cul-de-sac sample.

[yu-ta-sato]

Tests and docstrings are done except linting. Tests are not covered for line 606 and 619 in elements.py from the test cases prepared. Also, shapely.union_all is failed in the test of 311-oldest. I'd like to follow your direction of how to deal with shapely >=2 so it would be helpful if you share whether to modify the test env or test codes.

[jGaboardi]

@yu-ta-sato -- this is awesome!

Also, shapely.union_all is failed in the test of 311-oldest. I'd like to follow your direction of how to deal with shapely >=2 so it would be helpful if you share whether to modify the test env or test codes.

I think these failures are actually happening due to geopandas<1 in the oldest environment. However, considering SPEC000 we'll probably want drop support for geopandas<1. I open another issue explicitly for this and we'll get @martinfleis's thoughts on how to handle this specific case.

gantt
dateFormat YYYY-MM-DD
axisFormat %m / %Y
title Support Window

section geopandas
0.13.0 : 2023-05-06,2025-05-05
0.14.0 : 2023-09-15,2025-09-14
1.0.0 : 2024-06-24,2026-06-24
1.1.0 : 2025-06-01,2027-06-01

Loading