Deltares / MeshKernel

Deltares C++library for creating and editing meshes. It supports 1D & 2D unstructured meshes as well as curvilinear meshes.
https://deltares.github.io/MeshKernel
GNU Lesser General Public License v2.1
28 stars 8 forks source link

mesh2d_delete doesn't delete mesh2d_node_z #326

Open DanielTollenaar opened 5 months ago

DanielTollenaar commented 5 months ago

Describe the bug When clipping a netcdf over a geometrylist the mesh2d_node_x, mesh2d_node_y are clipped, but mesh2d_node_z isn't

To Reproduce Steps to reproduce the behavior:

  1. Download the mesh and reproduce-script we've created: https://we.tl/t-XHuUUF83RC
  2. Run the clip_mesh.py executing the code below (we use hydrolib-core to read the nc-file):
    
    from pathlib import Path
    import numpy as np
    from meshkernel import GeometryList
    import meshkernel as mk
    from hydrolib.core.dflowfm.net.models import Network

nc_path = Path(file).parent / "ovd-j98_6-v1a_net.nc"

read mesh using hydrolib-core

network = Network.from_file(Path(nc_path)) mesh2d = network._mesh2d

show node x/y/z

print("before") print(f"x: {mesh2d.mesh2d_node_x.shape}") print(f"y: {mesh2d.mesh2d_node_y.shape}") print(f"z: {mesh2d.mesh2d_node_z.shape}")

define gometrylist

x_coordinates = np.array( [219032.36821211, 219415.05146748, 236135.84080332, 257422.06085974, 257935.46607711, 221483.69564399, 219032.36821211], dtype=np.double)

y_coordinates = np.array( [502669.83347292, 504635.40577414, 520870.43426807, 522115.73282761, 511205.87195854, 495803.7154375, 502669.83347292], dtype=np.double)

geometrylist = GeometryList(x_coordinates, y_coordinates)

clip mesh

mesh2d.meshkernel.mesh2d_delete( geometry_list=geometrylist, delete_option=mk.DeleteMeshOption.INSIDE_NOT_INTERSECTED, invert_deletion=True, )

show_result

print("after") print(f"x: {mesh2d.mesh2d_node_x.shape}") print(f"y: {mesh2d.mesh2d_node_y.shape}") print(f"z: {mesh2d.mesh2d_node_z.shape}")


3. Note the size of `mesh_node2d_x`, `mesh_node2d_y` and `mesh_node2d_z` before and after `mesh2d_delete`:
![image](https://github.com/Deltares/MeshKernel/assets/14051252/401bdc31-59a0-4225-af34-95a426b9b178)

**Expected behavior**
A length of `mesh2d_node_z` of 226161 as are  `mesh2d_node_x` and `mesh2d_node_y`. 

**Screenshots**
We do think deletion of x and y is as expected. Running `network.plot()` after executing the script above gives a plausible result:
![image](https://github.com/Deltares/MeshKernel/assets/14051252/e56cf39a-7fe9-423a-b29a-d663760acfca)

**Version info (please complete the following information):**
 - OS: Windows
 - Versions:
   - Python: 3.8.19
   - hydrolib.core: 0.7.0
   - meshkernel: 4.1.0
veenstrajelmer commented 3 months ago

It is a bit tricky where this issue should land. It is indeed about meshkernel(py), but the issue lays in the fact that meshkernel does not support z-values at the moment. Hydrolib-core does have z-values, but these are the only properties that have to be carried around next to the meshkernel instance. They cannot be clipped with the meshkernel algorithm at the moment. I think the best way to fix this properly is to add z-values in meshkernel(py), but I am not sure if that is within the scope. Alternatively, we could do some interpolation of z-values in hydrolib core, since the new nodes/faces will be on the exact same locations as the old ones, this will probably result in identical values, but it would also require refactorings in other functions which would be more work than adding it in meshkernel. I will discuss this with @lucacarniato