Behaviour of clipping to a polygon is limited due to new options of DeleteMeshOption

[RuudHurkmans]

Problem
When clipping to a polygon the old behaviour of clipping to the exact boundar of a polygon is not possible anymore.

The code below plots a mesh with a polygon clipped out of it. Depending on the DMO the behaviour is different.

dmo = mk.DeleteMeshOption.INSIDE_NOT_INTERSECTED

dmo = mk.DeleteMeshOption.INSIDE_AND_INTERSECTED

Both are not the behavirour of the old version with dmo = ALL_FACE_CIRCUMCENTERS:

(in this image a multipolygon was used, but the effect is the same)

** Code used to reproduce these images**

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

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

fmmodel = FMModel()
network = fmmodel.geometry.netfile.network
   
rectangle = box(-10, -10, 10, 10)
clipgeo = box(-6, -1, -4, 2)

dmo = mk.DeleteMeshOption.INSIDE_NOT_INTERSECTED

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

network.mesh2d_create_rectilinear_within_extent(rectangle.bounds, dx=1, dy=1)

network.mesh2d_clip_mesh(from_polygon(clipgeo), dmo, 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
)

edge_nodes = network._mesh2d.mesh2d_edge_nodes
lc_mesh2d = LineCollection(nodes2d[edge_nodes], **mesh2d_kwargs)
ax.add_collection(lc_mesh2d)
ax.plot(*clipgeo.exterior.coords.xy, color="r", ls="--")
plt.show()

Describe the solution you'd like
Hard to say whether this is critical, but in complex geometries of regional water courses it would be good to be able to clip the boundary as accurately as possible (they can be quite narrow). So I suggest to either investigate whether it is possible to reproduce the old behaviour or re-introduce the old deletemeshoption in meshkernel.

[veenstrajelmer]

@RuudHurkmans: meshkernelpy 4.1.0 has an additional DeleteMeshOption called FACES_WITH_INCLUDED_CIRCUMCENTERS. This deletes all cells of which the circumcenters are inside the polygon. This was also available in meshkernel<3..0.0 under a slightly different name (ALL_FACE_CIRCUMCENTERS). The result is the following:

This version of meshkernel will be supported in the next release of HYDROLIB-core, which I expect this week. If you want to test it already, you can use the main branch since #610 has already been merged.