diff --git a/meshkernel/meshkernel.py b/meshkernel/meshkernel.py
index 439621f4..6098d4bb 100644
--- a/meshkernel/meshkernel.py
+++ b/meshkernel/meshkernel.py
@@ -1648,6 +1648,45 @@ def mesh2d_get_orthogonality(self) -> GeometryList:
 
         return geometry_list_out
 
+    def mesh2d_get_property(self, property: Mesh2d.Property) -> GeometryList:
+        """Gets the polygons matching the metric value within the minimum and maximum value.
+
+        Args:
+
+            property (Mesh2d.Property): The property to retrieve
+
+        Returns:
+            GeometryList: The resulting geometry list containing the value of the properties
+        """
+
+        c_geometry_list_dimension = c_int()
+
+        self._execute_function(
+            self.lib.mkernel_mesh2d_get_property_dimension,
+            self._meshkernelid,
+            c_int(property),
+            byref(c_geometry_list_dimension),
+        )
+
+        n_coordinates = c_geometry_list_dimension.value
+        x_coordinates = np.empty(n_coordinates, dtype=np.double)
+        y_coordinates = np.empty(n_coordinates, dtype=np.double)
+        values = np.empty(n_coordinates, dtype=np.double)
+
+        property_list = GeometryList(
+            x_coordinates=x_coordinates, y_coordinates=y_coordinates, values=values
+        )
+        c_property_list = CGeometryList.from_geometrylist(property_list)
+
+        self._execute_function(
+            self.lib.mkernel_mesh2d_get_property,
+            self._meshkernelid,
+            c_int(property),
+            byref(c_property_list),
+        )
+
+        return property_list
+
     def mesh2d_get_smoothness(self):
         """Gets the smoothness, expressed as the ratio between the values of two neighboring faces areas.
 
diff --git a/meshkernel/py_structures.py b/meshkernel/py_structures.py
index 6d240502..b935bcc2 100644
--- a/meshkernel/py_structures.py
+++ b/meshkernel/py_structures.py
@@ -123,6 +123,7 @@ class Property(IntEnum):
         """Different properties on a 2D mesh."""
 
         ORTHOGONALITY = 0
+        EDGE_LENGTHS = 1
 
     def __init__(
         self,
diff --git a/tests/test_mesh2d_basics.py b/tests/test_mesh2d_basics.py
index 95217795..c2262eec 100644
--- a/tests/test_mesh2d_basics.py
+++ b/tests/test_mesh2d_basics.py
@@ -446,16 +446,16 @@ def test_mesh2d_get_node_index_no_node_in_search_radius(
     (
         True,
         DeleteMeshOption.INSIDE_NOT_INTERSECTED,
-        4,
-        4,
-        1,
+        16,
+        24,
+        9,
     ),
     (
         True,
         DeleteMeshOption.INSIDE_AND_INTERSECTED,
-        16,
-        24,
-        9,
+        4,
+        4,
+        1,
     ),
     (
         False,
@@ -2308,6 +2308,93 @@ def test_mesh2d_deletion_and_get_orthogonality(
     assert len(values) == len(mesh2d.edge_x)
 
 
+cases_get_property = [
+    (
+        Mesh2d.Property.ORTHOGONALITY,
+        np.array(
+            [
+                -999.0,
+                0.0,
+                0.0,
+                -999.0,
+                -999.0,
+                0.0,
+                0.0,
+                -999.0,
+                -999.0,
+                0.0,
+                0.0,
+                -999.0,
+                -999.0,
+                -999.0,
+                -999.0,
+                0.0,
+                0.0,
+                0.0,
+                0.0,
+                0.0,
+                0.0,
+                -999.0,
+                -999.0,
+                -999.0,
+            ],
+            dtype=np.double,
+        ),
+    ),
+    (
+        Mesh2d.Property.EDGE_LENGTHS,
+        np.array(
+            [
+                100.0,
+                100.0,
+                100.0,
+                100.0,
+                100.0,
+                100.0,
+                100.0,
+                100.0,
+                100.0,
+                100.0,
+                100.0,
+                100.0,
+                50.0,
+                50.0,
+                50.0,
+                50.0,
+                50.0,
+                50.0,
+                50.0,
+                50.0,
+                50.0,
+                50.0,
+                50.0,
+                50.0,
+            ],
+            dtype=np.double,
+        ),
+    ),
+]
+
+
+@pytest.mark.parametrize(
+    "property, expected_values",
+    cases_get_property,
+)
+def test_mesh2d_get_property(
+    meshkernel_with_mesh2d: MeshKernel,
+    property: Mesh2d.Property,
+    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)
+
+    assert property_list.values == approx(expected_values, 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
diff --git a/tests/test_mesh2d_refinement.py b/tests/test_mesh2d_refinement.py
index f5df66ed..1409b762 100644
--- a/tests/test_mesh2d_refinement.py
+++ b/tests/test_mesh2d_refinement.py
@@ -1,14 +1,6 @@
 import numpy as np
-import pytest
-from numpy import ndarray
 
-from meshkernel import (
-    GeometryList,
-    MakeGridParameters,
-    Mesh2d,
-    MeshKernel,
-    MeshRefinementParameters,
-)
+from meshkernel import GeometryList, MakeGridParameters, MeshKernel
 
 
 def test_mesh2d_casulli_refinement():