GRIDEDIT-1548: Add Casulli refinement based on depths

meshkernel/__init__.py

@@ -11,6 +11,7 @@
     GeometryList,
     GriddedSamples,
     InterpolationValues,
+    InterpolationType,
     MakeGridParameters,
     Mesh1d,
     Mesh2d,

meshkernel/meshkernel.py

@@ -41,6 +41,7 @@
     DeleteMeshOption,
     GeometryList,
     GriddedSamples,
+    InterpolationType,
     MakeGridParameters,
     Mesh1d,
     Mesh2d,
@@ -1585,6 +1586,75 @@ def mesh2d_casulli_refinement_on_polygon(
             byref(c_polygon),
         )
 
+    def mkernel_delete_property(self, propertyId: int):
+        """
+        Deletes a property and its calculator.
+
+        Args:
+            propertyId (int): The property id.
+        """
+
+        self._execute_function(
+            self.lib.mkernel_deallocate_property,
+            c_int(propertyId),
+        )
+
+    def mkernel_set_property(self,
+                             projection: ProjectionType,
+                             interpolationType: InterpolationType,
+                             sampleData: GeometryList) -> int:
+
+        """
+        Sets the property data for the mesh, the sample data points do not have to match the mesh2d nodes.
+
+        Args:
+            projection (ProjectionType): The projection type used by the sample data.
+            interpolationType (InterpolationType): The type of interpolation required, triangulation or averaging.
+            sampleData (GeometryList): The sample data and associated sample data points
+
+        Returns:
+            int: The id of the property.
+        """
+
+        c_sampleData = CGeometryList.from_geometrylist(sampleData)
+        propertyId = c_int()
+
+        self._execute_function(
+            self.lib.mkernel_mesh2d_set_property,
+            c_int(projection.value),
+            interpolationType,
+            byref(c_sampleData),
+            byref(propertyId)
+        )
+
+        return propertyId.value
+
+    def mkernel_mesh2d_casulli_refinement_wrt_depths(self,
+                                                     polygons: GeometryList,
+                                                     propertyId: int,
+                                                     meshRefinementParameters: MeshRefinementParameters,
+                                                     minimumRefinementDepth: float ) -> None:
+
+        """
+        Refine mesh using the Casulli refinement algorithm based on the depth values.
+
+        Args:
+            polygons (GeometryList): The polygon within which the refinement is computed.
+            propertyId (int): The identifier of the property and its associated interpolator.
+            meshRefinementParameters (MeshRefinementParameters): Parameters indicating how the mesh is to be refined.
+            minimumRefinementDepth (float): Nodes with depth value less than this value will not be marked for refinement.
+        """
+
+        c_polygons = CGeometryList.from_geometrylist(polygons)
+        self._execute_function(
+            self.lib.mkernel_mesh2d_casulli_refinement_wrt_depths,
+            self._meshkernelid,
+            byref(c_polygons),
+            c_int(propertyId),
+            byref(meshRefinementParameters),
+            float(minimumRefinementDepth))
+
+
     def mesh2d_compute_orthogonalization(
         self,
         project_to_land_boundary_option: ProjectToLandBoundaryOption,
@@ -1648,7 +1718,7 @@ def mesh2d_get_orthogonality(self) -> GeometryList:
 
         return geometry_list_out
 
-    def mesh2d_get_property(self, property: Mesh2d.Property) -> GeometryList:
+    def mesh2d_get_property(self, mesh2dLocation: Mesh2dLocation, property: Mesh2d.Property) -> GeometryList:
         """Gets the polygons matching the metric value within the minimum and maximum value.
 
         Args:
@@ -1682,6 +1752,7 @@ def mesh2d_get_property(self, property: Mesh2d.Property) -> GeometryList:
             self.lib.mkernel_mesh2d_get_property,
             self._meshkernelid,
             c_int(property),
+            c_int(mesh2dLocation.value),
             byref(c_property_list),
         )
 

meshkernel/py_structures.py

@@ -57,7 +57,6 @@ class Mesh2dLocation(IntEnum):
     NODES = 1
     EDGES = 2
 
-
 @unique
 class AveragingMethod(IntEnum):
     """The averaging methods."""
@@ -99,6 +98,12 @@ class InterpolationValues(IntEnum):
     INT = 2
     DOUBLE = 3
 
+@unique
+class InterpolationType(IntEnum):
+    """The Mesh2d location types."""
+
+    TRIANGULATION = 0
+    AVERAGING = 1
 
 class Mesh2d:
     """This class is used for getting and setting two-dimensional mesh data.

tests/test_mesh2d_basics.py

@@ -2311,6 +2311,7 @@ def test_mesh2d_deletion_and_get_orthogonality(
 cases_get_property = [
     (
         Mesh2d.Property.ORTHOGONALITY,
+        Mesh2dLocation.EDGES,
         np.array(
             [
                 -999.0,
@@ -2343,6 +2344,7 @@ def test_mesh2d_deletion_and_get_orthogonality(
     ),
     (
         Mesh2d.Property.EDGE_LENGTHS,
+        Mesh2dLocation.EDGES,
         np.array(
             [
                 100.0,
@@ -2383,14 +2385,15 @@ def test_mesh2d_deletion_and_get_orthogonality(
 def test_mesh2d_get_property(
     meshkernel_with_mesh2d: MeshKernel,
     property: Mesh2d.Property,
+    mesh2dLocation: Mesh2dLocation,
     expected_values: ndarray,
 ):
     """Test mesh2d_get_property,
     getting the mesh2d property values
     """
     mk = meshkernel_with_mesh2d(rows=3, columns=3, spacing_x=50.0, spacing_y=100.0)
 
-    property_list = mk.mesh2d_get_property(property)
+    property_list = mk.mesh2d_get_property(mesh2dLocation, property)
 
     assert property_list.values == approx(expected_values, abs=1e-6)
 
@@ -2421,3 +2424,31 @@ def test_mesh2d_get_filtered_face_polygons():
 
     assert face_polygons.x_coordinates == approx(expected_coordinates_x, abs=1e-6)
     assert face_polygons.y_coordinates == approx(expected_coordinates_y, abs=1e-6)
+
+
+def test_mesh2d_get_filtered_face_polygons():
+    """Test mesh2d_get_filtered_face_polygons,
+    getting the polygons of faces with all edges having bad orthogonality values
+    """
+    mk = MeshKernel()
+
+    edge_nodes = np.array(
+        [0, 1, 1, 2, 2, 3, 0, 3, 1, 4, 0, 4, 0, 5, 3, 5, 3, 6, 2, 6, 2, 7, 1, 7],
+        dtype=np.int32,
+    )
+
+    node_x = np.array([57.0, 49.1, 58.9, 66.7, 48.8, 65.9, 67.0, 49.1], dtype=np.double)
+    node_y = np.array([23.6, 14.0, 6.9, 16.2, 23.4, 24.0, 7.2, 6.7], dtype=np.double)
+
+    input_mesh2d = Mesh2d(node_x, node_y, edge_nodes)
+    mk.mesh2d_set(input_mesh2d)
+
+    face_polygons = mk.mesh2d_get_filtered_face_polygons(
+        Mesh2d.Property.ORTHOGONALITY, 0.04, 1.0
+    )
+
+    expected_coordinates_x = np.array([57.0, 49.1, 58.9, 66.7, 57.0], dtype=np.double)
+    expected_coordinates_y = np.array([23.6, 14.0, 6.9, 16.2, 23.6], dtype=np.double)
+
+    assert face_polygons.x_coordinates == approx(expected_coordinates_x, abs=1e-6)
+    assert face_polygons.y_coordinates == approx(expected_coordinates_y, abs=1e-6)