Mesh generation in a polygon with holes

[RuudHurkmans]

Generate a mesh in a single complex polygon

Based on earlier test with simple polygons, I rewrote code to generate a mesh within a polygon by generating a mesh within the polygon bounds and then clipping everything outside.

Furthermore, I think an third option for the deletemeshoption would be useful for this case as well.

It's almost correct but I don't understand the vertical lines.

Code (Sorry, it's a bit much but i needed to include some helper functions)

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")#.iloc[0]
#bbox = box(198000, 392000,202000,398000)
bbox = box(180000, 380000,190000,400000)

#extent = extent.geometry.clip( bbox).iloc[0]
extent = extent.geometry.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 it be that the polygon is not contiguous? For instance at the vertical lines the polygon is actually disconnected and therefore consists of multiple polygons? Also, could you attach the required data to the issue?

[RuudHurkmans]

@veenstrajelmer not that I can see. It's just one single polygon, not a multigeometry feature. I forgot to attach the used geometry-file. Here it is:

2D_roostergebied_v5.zip

[veenstrajelmer]

When using meshkernel 4.2.0, the result seems to be as expected:

The minimal meshkernel version will be bumped as part of #692.