Mesh generation in a multipolygon

[RuudHurkmans]

Problem

In the current version of D-Hydamo, when the extent is a multipolygon, a mesh is created for each polygons bounds and then clipped.

In trying to reproduce it, when I now clip the 2nd mesh, every mesh from previous iterations is clipped and only the last mesh remains:

For illustration: the mesh based on each polygon's bounds without clipping:

There is probably a better way of doing this I'm open to suggestions.

The used shapefile is here:
2D_roostergebied_v5.zip

Code
The code is the same as in issue #614 , but with a different extent which is now a multipolygon:

from shapely.geometry import Point, Polygon, MultiPolygon, MultiLineString, LineString, box
from matplotlib.collections import LineCollection
import numpy as np
import meshkernel as mk
import matplotlib.pyplot as plt
import geopandas as gpd
from hydrolib.core.dflowfm.mdu.models import FMModel
from hydrolib.core.dflowfm.net.models import Mesh2d
from matplotlib.path import Path
from matplotlib.patches import PathPatch
from matplotlib.collections import PatchCollection

def _as_geometry_list(geometry, singletype, multitype):   
    if isinstance(geometry, singletype):
        return [geometry]
    elif isinstance(geometry, multitype):
        return [p for p in geometry.geoms]
    elif isinstance(geometry, list):
        lst = []
        for item in geometry:
            lst.extend(_as_geometry_list(item, singletype, multitype))
        return lst
    else:
        raise TypeError(f'Expected {singletype} or {multitype}. Got "{type(geometry)}"')
        

def as_polygon_list(polygon):   
    return _as_geometry_list(polygon, Polygon, MultiPolygon)

def from_polygon(polygon): 
    inner_outer_separator = -998
    # Add exterior
    x_ext, y_ext = np.array(polygon.exterior.coords[:]).T
    x_crds = [x_ext]
    y_crds = [y_ext]

    # Add interiors, seperated by inner_outer_separator
    for interior in polygon.interiors:
        x_int, y_int = np.array(interior.coords[:]).T
        x_crds.append([inner_outer_separator])
        x_crds.append(x_int)
        y_crds.append([inner_outer_separator])
        y_crds.append(y_int)
    gl = mk.GeometryList(
        x_coordinates=np.concatenate(x_crds), y_coordinates=np.concatenate(y_crds), inner_outer_separator=-998
    )
    return gl

extent = gpd.read_file(r"hydrolib/notebooks/data/2D_roostergebied_v5.shp")
bbox = box(198000, 392000,202000,398000)
extent = extent.geometry.clip( bbox).iloc[0]

fmmodel = FMModel()
network = fmmodel.geometry.netfile.network

def mesh_in_multipolgyon(network, multipolygon, dx, dy, dmo):
    mesh2d_hc = Mesh2d(meshkernel=mk.MeshKernel())
    if isinstance(multipolygon, Polygon):
        multipolygon=MultiPolygon([multipolygon])
    
    plist = multipolygon.geoms
    if len(plist) > 1:    
        for part in plist:     
            try:       
                xmin, ymin, xmax, ymax = part.bounds
                rows = int((ymax - ymin) / dy)
                columns = int((xmax - xmin) / dx)
                if rows < 1 and columns < 1:
                    break           
                mesh2d_hc.create_rectilinear(extent=part.bounds, dx=dx, dy=dy)                       
                mesh2d_hc.clip(from_polygon(part),dmo, False)
            except:
                print(f'Could not create mesh for this sub-polygon')
        mesh2d_output = mesh2d_hc.get_mesh2d()
    else:
        part = plist[0]
        mesh2d_hc.create_rectilinear(extent=part.bounds, dx=dx, dy=dy)
        mesh2d_hc.clip(from_polygon(part), dmo, False)
        mesh2d_output = mesh2d_hc.get_mesh2d()
    network._mesh2d._set_mesh2d(mesh2d_output.node_x, mesh2d_output.node_y, mesh2d_output.edge_nodes)
    return network

def plot_polygon(ax, poly, **kwargs):
    path = Path.make_compound_path(
        Path(np.asarray(poly.exterior.coords)[:, :2]),
        *[Path(np.asarray(ring.coords)[:, :2]) for ring in poly.interiors])

    patch = PathPatch(path, **kwargs)
    collection = PatchCollection([patch], **kwargs)
    
    ax.add_collection(collection, autolim=True)
    ax.autoscale_view()
    return collection

dmo = mk.DeleteMeshOption.INSIDE_NOT_INTERSECTED
network = mesh_in_multipolgyon(network, extent, 100, 100, dmo)

fig, ax = plt.subplots()
ax.set_aspect(1.0)
nodes2d = np.stack(
    [network._mesh2d.mesh2d_node_x, network._mesh2d.mesh2d_node_y], axis=1
)
mesh2d_kwargs = {"color": "r", "lw": 0.5}
edge_nodes = network._mesh2d.mesh2d_edge_nodes
lc_mesh2d = LineCollection(nodes2d[edge_nodes], **mesh2d_kwargs)
ax.add_collection(lc_mesh2d)
if isinstance(extent, Polygon):
    plot_polygon(ax, extent, facecolor='lightblue', edgecolor='lightblue')
else:
    for part in extent.geoms:
        plot_polygon(ax, part, facecolor='lightblue', edgecolor='lightblue')
ax.set_xlim(extent.bounds[0], extent.bounds[2])
ax.set_ylim(extent.bounds[1], extent.bounds[3])
plt.show()

[veenstrajelmer]

@RuudHurkmans could you maybe add the reproducible code to the issue description to make this issue self-describing? Could you also add the shapefile?

Last but not least, what is the expected/desired behaviour? Would you expect mesh in all the blue area? Would you also require mesh on the edges of the polygon or not?

[RuudHurkmans]

I added the complete code and the shapefile .

I would indeed expect a mesh in alle the blue polygons, but whether there are cells on the edges depends on the DMO I use. The issue here is that in every polygon in the iteration, previous results are clipped. I'm mainly looking for a way to prevent that.

[veenstrajelmer]

@RuudHurkmans the below code sort of does what you need. However, I have to bypass hydrolib-core in order to achieve this and in my view it is quite complex to get the desired result. If I do call hydrolib-core directly with network._mesh2d.clip(clip_pol_all, dmo, False) it raises "NotImplementedError: Deleting outside more than a single exterior (MultiPolygon) is not implemented."

from shapely.geometry import Polygon, MultiPolygon, box
from matplotlib.collections import LineCollection
import numpy as np
import meshkernel as mk
import matplotlib.pyplot as plt
plt.close("all")
import geopandas as gpd
from hydrolib.core.dflowfm.mdu.models import FMModel
from hydrolib.core.dflowfm.net.models import Mesh2d
from matplotlib.path import Path
from matplotlib.patches import PathPatch
from matplotlib.collections import PatchCollection

def _as_geometry_list(geometry, singletype, multitype):   
    if isinstance(geometry, singletype):
        return [geometry]
    elif isinstance(geometry, multitype):
        return [p for p in geometry.geoms]
    elif isinstance(geometry, list):
        lst = []
        for item in geometry:
            lst.extend(_as_geometry_list(item, singletype, multitype))
        return lst
    else:
        raise TypeError(f'Expected {singletype} or {multitype}. Got "{type(geometry)}"')
        

def as_polygon_list(polygon):   
    return _as_geometry_list(polygon, Polygon, MultiPolygon)

def from_polygon(polygon): 
    inner_outer_separator = -998
    # Add exterior
    x_ext, y_ext = np.array(polygon.exterior.coords[:]).T
    x_crds = [x_ext]
    y_crds = [y_ext]

    # Add interiors, seperated by inner_outer_separator
    for interior in polygon.interiors:
        x_int, y_int = np.array(interior.coords[:]).T
        x_crds.append([inner_outer_separator])
        x_crds.append(x_int)
        y_crds.append([inner_outer_separator])
        y_crds.append(y_int)
    gl = mk.GeometryList(
        x_coordinates=np.concatenate(x_crds), y_coordinates=np.concatenate(y_crds), inner_outer_separator=-998
    )
    return gl

def from_polygon_multi(polygon): 
    poly_separator = -999
    # Add exterior
    x_crds = []
    y_crds = []

    # Add interiors, seperated by inner_outer_separator
    for ipol, polygeom in enumerate(polygon.geoms):
        x_int, y_int = np.array(polygeom.exterior.coords[:]).T
        if ipol>0:
            x_crds.append([poly_separator])
            y_crds.append([poly_separator])
        x_crds.append(x_int)
        y_crds.append(y_int)
    gl = mk.GeometryList(
        x_coordinates=np.concatenate(x_crds), y_coordinates=np.concatenate(y_crds), inner_outer_separator=-998
    )
    return gl

extent = gpd.read_file(r"c:\Users\veenstra\Downloads\2D_roostergebied_v5\2D_roostergebied_v5.shp")
bbox = box(198000, 392000,202000,398000)
extent = extent.geometry.clip(bbox).iloc[0]
# [2D_roostergebied_v5.zip](https://github.com/Deltares/HYDROLIB-core/files/14330025/2D_roostergebied_v5.zip)

fmmodel = FMModel()
network = fmmodel.geometry.netfile.network

def mesh_in_multipolgyon(network, multipolygon, dx, dy, dmo):
    mesh2d_hc = Mesh2d(meshkernel=mk.MeshKernel())
    if isinstance(multipolygon, Polygon):
        multipolygon=MultiPolygon([multipolygon])
    
    plist = multipolygon.geoms
    for part in plist:     
        try:       
            xmin, ymin, xmax, ymax = part.bounds
            rows = int((ymax - ymin) / dy)
            columns = int((xmax - xmin) / dx)
            if rows < 1 and columns < 1:
                break
            # clip_pol_all.append(from_polygon(part))
            mesh2d_hc.create_rectilinear(extent=part.bounds, dx=dx, dy=dy)                       
        except:
            print('Could not create mesh for this sub-polygon')
    
    mesh2d_output = mesh2d_hc.get_mesh2d()
    network._mesh2d._set_mesh2d(mesh2d_output.node_x, mesh2d_output.node_y, mesh2d_output.edge_nodes)
    return network


def plot_polygon(ax, poly, **kwargs):
    path = Path.make_compound_path(
        Path(np.asarray(poly.exterior.coords)[:, :2]),
        *[Path(np.asarray(ring.coords)[:, :2]) for ring in poly.interiors])

    patch = PathPatch(path, **kwargs)
    collection = PatchCollection([patch], **kwargs)
    
    ax.add_collection(collection, autolim=True)
    ax.autoscale_view()
    return collection

dmo = mk.DeleteMeshOption.INSIDE_NOT_INTERSECTED
network = mesh_in_multipolgyon(network, extent, 100, 100, dmo)
clip_pol_all = from_polygon_multi(extent)
# network._mesh2d.clip(clip_pol_all, dmo, False)
network._mesh2d.meshkernel.mesh2d_delete(
    geometry_list=clip_pol_all,
    delete_option=dmo,
    invert_deletion=True,
)


fig, ax = plt.subplots()
ax.set_aspect(1.0)
nodes2d = np.stack(
    [network._mesh2d.mesh2d_node_x, network._mesh2d.mesh2d_node_y], axis=1
)
mesh2d_kwargs = {"color": "r", "lw": 0.5}
edge_nodes = network._mesh2d.mesh2d_edge_nodes
lc_mesh2d = LineCollection(nodes2d[edge_nodes], **mesh2d_kwargs)
ax.add_collection(lc_mesh2d)
if isinstance(extent, Polygon):
    plot_polygon(ax, extent, facecolor='lightblue', edgecolor='lightblue')
else:
    for part in extent.geoms:
        plot_polygon(ax, part, facecolor='lightblue', edgecolor='lightblue')
ax.set_xlim(extent.bounds[0], extent.bounds[2])
ax.set_ylim(extent.bounds[1], extent.bounds[3])
plt.show()

Gives:

Is this indeed what you are looking for? Except for the horizontal mesh that remains, which might be related to #614 Would you also expect there would be no mesh present in the holes? Either way, I guess the NotImplementedError in HYDROLIB-core could be picked up for implementation in the future if that is desired.

[RuudHurkmans]

@veenstrajelmer Yes that's wat i'm looking for. Ideallly the code would be able to deal with holes ( I was getting the same NotImplementedError), but if that's something for a later release it's fine, it's not critical for now. I hink for practical use in D-Hydamo issue #614 is most critical now.

[guydup]

I tried fixing the issue in #614 and the holes in this issue. The holes can be implemented for now with the code below. The issue in #614 is still present, even after manually finding the wrong faces and deleting them from the nodes/edges lists. It seems like the weird faces are produced in network._mesh2d._set_mesh2d()?

import geopandas as gpd
import matplotlib.pyplot as plt
import meshkernel as mk
import numpy as np
from matplotlib.collections import LineCollection
from shapely.geometry import MultiPolygon, Polygon, box

from hydrolib.core.dflowfm.mdu.models import FMModel
from hydrolib.core.dflowfm.net.models import Mesh2d


def from_polygon_multi(polygon, mode):
    poly_separator = -999
    x_crds, y_crds = [], []
    for polygeom in polygon.geoms:
        x_geom, y_geom = [], []
        if mode == "exterior":
            x_ext, y_ext = np.array(polygeom.exterior.coords[:]).T
            x_geom.append(x_ext)
            y_geom.append(y_ext)
        elif mode == "interior":
            for int_geom in polygeom.interiors:
                x_int, y_int = np.array(int_geom.coords[:]).T
                x_geom.append(x_int)
                y_geom.append(y_int)

        for xc, yc in zip(x_geom, y_geom):
            if len(x_crds) > 0:
                x_crds.append([poly_separator])
                y_crds.append([poly_separator])
            x_crds.append(xc)
            y_crds.append(yc)
    gl = mk.GeometryList(
        x_coordinates=np.concatenate(x_crds),
        y_coordinates=np.concatenate(y_crds),
        inner_outer_separator=-998,
    )
    return gl


def mesh_in_multipolgyon(multipolygon, dx, dy, delete="INSIDE_NOT_INTERSECTED"):
    if delete == "INSIDE_NOT_INTERSECTED":
        dmo_ext = mk.DeleteMeshOption.INSIDE_NOT_INTERSECTED
        dmo_int = mk.DeleteMeshOption.INSIDE_AND_INTERSECTED
    elif delete == "INSIDE_AND_INTERSECTED":
        dmo_ext = mk.DeleteMeshOption.INSIDE_AND_INTERSECTED
        dmo_int = mk.DeleteMeshOption.INSIDE_NOT_INTERSECTED
    else:
        raise ValueError(f"Unknown delete method '{delete}'")

    fmmodel = FMModel()
    network = fmmodel.geometry.netfile.network
    mesh2d_hc = Mesh2d(meshkernel=mk.MeshKernel())

    # Enforce multipolygon
    if isinstance(multipolygon, Polygon):
        multipolygon = MultiPolygon([multipolygon])

    for part in multipolygon.geoms:
        try:
            xmin, ymin, xmax, ymax = part.bounds
            rows = int((ymax - ymin) / dy)
            columns = int((xmax - xmin) / dx)
            if rows < 1 and columns < 1:
                break
            mesh2d_hc.create_rectilinear(extent=part.bounds, dx=dx, dy=dy)
        except Exception:
            print("Could not create mesh for this sub-polygon")

    mesh2d_output = mesh2d_hc.get_mesh2d()
    network._mesh2d._set_mesh2d(
        mesh2d_output.node_x, mesh2d_output.node_y, mesh2d_output.edge_nodes
    )

    # Keep polygon area
    clip_pol_ext = from_polygon_multi(extent, "exterior")
    network._mesh2d.meshkernel.mesh2d_delete(
        geometry_list=clip_pol_ext,
        delete_option=dmo_ext,
        invert_deletion=True,
    )

    # Clip holes
    clip_pol_int = from_polygon_multi(extent, "interior")
    network._mesh2d.meshkernel.mesh2d_delete(
        geometry_list=clip_pol_int,
        delete_option=dmo_int,
        invert_deletion=False,
    )

    # Find mesh faces that are not slivers
    node_x_org = network._mesh2d.mesh2d_node_x
    node_y_org = network._mesh2d.mesh2d_node_y
    nodes2d_org = np.stack([node_x_org, node_y_org], axis=1)
    edge_nodes_org = network._mesh2d.mesh2d_edge_nodes
    face_nodes_org = network._mesh2d.mesh2d_face_nodes
    grouped_nodes = nodes2d_org[face_nodes_org]
    geoms = gpd.GeoSeries(
        [Polygon(grouped_nodes[i, :, :]) for i in range(grouped_nodes.shape[0])]
    )
    if delete == "INSIDE_NOT_INTERSECTED":
        idx = geoms.sindex.query(extent, predicate="contains")
    else:
        idx = geoms.sindex.query(extent, predicate="intersects")

    # remove nodes and edges associated with the slivers
    kn = np.unique(face_nodes_org[idx, :].flatten())
    node_x = node_x_org[kn]
    node_y = node_y_org[kn]
    ke = np.isin(edge_nodes_org, kn).all(axis=1)
    edge_nodes = edge_nodes_org[ke, :]
    new_ind = dict(zip(kn, range(len(kn))))
    for i in range(edge_nodes.shape[0]):
        edge_nodes[i, 0] = new_ind[edge_nodes[i, 0]]
        edge_nodes[i, 1] = new_ind[edge_nodes[i, 1]]
    network._mesh2d._set_mesh2d(node_x, node_y, edge_nodes)

    fig, axs = plt.subplots(figsize=(10, 5), ncols=2)
    for ax in axs:
        gpd.GeoSeries([multipolygon]).plot(ax=ax, color="gray")
        geoms.plot(ax=ax, color="red")
        geoms[idx].plot(ax=ax, color="blue")
    # plot geopandas result
    axs[0].set_title("Geopandas fix")
    nodes2d = np.stack([node_x, node_y], axis=1)
    lc_mesh2d = LineCollection(nodes2d[edge_nodes], color="black", lw=1)
    axs[0].add_collection(lc_mesh2d)
    # plot manualy adjusted meshkernel result
    axs[1].set_title("Meshkernel result")
    nodes2d = np.stack(
        [network._mesh2d.mesh2d_node_x, network._mesh2d.mesh2d_node_y], axis=1
    )
    edge_nodes = network._mesh2d.mesh2d_edge_nodes
    lc_mesh2d = LineCollection(nodes2d[edge_nodes], color="black", lw=1)
    axs[1].add_collection(lc_mesh2d)
    fig.show()

    return network


extent = gpd.read_file("hydrolib/notebooks/data/2D_roostergebied_v5.shp")
bbox = box(198000, 392000, 202000, 398000)
extent = extent.geometry.clip(bbox).iloc[0]

network = mesh_in_multipolgyon(extent, 100, 100)