Open RuudHurkmans opened 7 months ago
veenstrajelmer
commented
7 months ago @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
commented
7 months ago 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
commented
7 months ago @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 https://github.com/Deltares/HYDROLIB-core/issues/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
commented
7 months ago @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
commented
7 months ago 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)When using meshkernel 4.2.0, the horizontal line does not show up anymore:
The minimal meshkernel version will be bumped as part of https://github.com/Deltares/HYDROLIB-core/issues/692.
However, to solve this issue the holes in the polygon should also not be filled with mesh. Also, this still requires calling meshkernel, while ideally all should not be necessary as the hydrolib-core functions take care of all desired behaviour.
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: