Refs #35058 -- Added support for measured geometries to GDAL Point.

Co-authored-by: Nick Pope <nick@nickpope.me.uk>

6 files changed, 157 insertions, 7 deletions

diff --git a/django/contrib/gis/gdal/geometries.py b/django/contrib/gis/gdal/geometries.py
index 57f81ca476..f4ba5816a2 100644
--- a/django/contrib/gis/gdal/geometries.py
+++ b/django/contrib/gis/gdal/geometries.py
@@ -287,6 +287,22 @@ class OGRGeometry(GDALBase):
         else:
             raise ValueError(f"Input to 'set_3d' must be a boolean, got '{value!r}'.")
 
+    @property
+    def is_measured(self):
+        """Return True if the geometry has M coordinates."""
+        return capi.is_measured(self.ptr)
+
+    def set_measured(self, value):
+        """Set if this geometry has M coordinates."""
+        if value is True:
+            capi.set_measured(self.ptr, 1)
+        elif value is False:
+            capi.set_measured(self.ptr, 0)
+        else:
+            raise ValueError(
+                f"Input to 'set_measured' must be a boolean, got '{value!r}'."
+            )
+
     # #### SpatialReference-related Properties ####
 
     # The SRS property
@@ -386,14 +402,22 @@ class OGRGeometry(GDALBase):
         sz = self.wkb_size
         # Creating the unsigned character buffer, and passing it in by reference.
         buf = (c_ubyte * sz)()
-        capi.to_wkb(self.ptr, byteorder, byref(buf))
+        # For backward compatibility, export old-style 99-402 extended
+        # dimension types when geometry does not have an M dimension.
+        # https://gdal.org/api/vector_c_api.html#_CPPv417OGR_G_ExportToWkb12OGRGeometryH15OGRwkbByteOrderPh
+        to_wkb = capi.to_iso_wkb if self.is_measured else capi.to_wkb
+        to_wkb(self.ptr, byteorder, byref(buf))
         # Returning a buffer of the string at the pointer.
         return memoryview(string_at(buf, sz))
 
     @property
     def wkt(self):
         "Return the WKT representation of the Geometry."
-        return capi.to_wkt(self.ptr, byref(c_char_p()))
+        # For backward compatibility, export old-style 99-402 extended
+        # dimension types when geometry does not have an M dimension.
+        # https://gdal.org/api/vector_c_api.html#_CPPv417OGR_G_ExportToWkt12OGRGeometryHPPc
+        to_wkt = capi.to_iso_wkt if self.is_measured else capi.to_wkt
+        return to_wkt(self.ptr, byref(c_char_p()))
 
     @property
     def ewkt(self):
@@ -569,10 +593,20 @@ class Point(OGRGeometry):
             return capi.getz(self.ptr, 0)
 
     @property
+    def m(self):
+        """Return the M coordinate for this Point."""
+        if self.is_measured:
+            return capi.getm(self.ptr, 0)
+
+    @property
     def tuple(self):
         "Return the tuple of this point."
+        if self.is_3d and self.is_measured:
+            return self.x, self.y, self.z, self.m
         if self.is_3d:
-            return (self.x, self.y, self.z)
+            return self.x, self.y, self.z
+        if self.is_measured:
+            return self.x, self.y, self.m
         return self.x, self.y
 
     coords = tuple
@@ -753,7 +787,9 @@ GEO_CLASSES = {
     6: MultiPolygon,
     7: GeometryCollection,
     101: LinearRing,
-    1 + OGRGeomType.wkb25bit: Point,
+    2001: Point,  # POINT M
+    3001: Point,  # POINT ZM
+    1 + OGRGeomType.wkb25bit: Point,  # POINT Z
     2 + OGRGeomType.wkb25bit: LineString,
     3 + OGRGeomType.wkb25bit: Polygon,
     4 + OGRGeomType.wkb25bit: MultiPoint,
diff --git a/django/contrib/gis/gdal/prototypes/geom.py b/django/contrib/gis/gdal/prototypes/geom.py
index aaada31d5b..6dc6967ecf 100644
--- a/django/contrib/gis/gdal/prototypes/geom.py
+++ b/django/contrib/gis/gdal/prototypes/geom.py
@@ -51,6 +51,7 @@ to_kml = string_output(
 getx = pnt_func(lgdal.OGR_G_GetX)
 gety = pnt_func(lgdal.OGR_G_GetY)
 getz = pnt_func(lgdal.OGR_G_GetZ)
+getm = pnt_func(lgdal.OGR_G_GetM)
 
 # Geometry creation routines.
 if GDAL_VERSION >= (3, 3):
@@ -82,6 +83,8 @@ geom_sym_diff = geom_output(lgdal.OGR_G_SymmetricDifference, [c_void_p, c_void_p
 geom_union = geom_output(lgdal.OGR_G_Union, [c_void_p, c_void_p])
 is_3d = bool_output(lgdal.OGR_G_Is3D, [c_void_p])
 set_3d = void_output(lgdal.OGR_G_Set3D, [c_void_p, c_int], errcheck=False)
+is_measured = bool_output(lgdal.OGR_G_IsMeasured, [c_void_p])
+set_measured = void_output(lgdal.OGR_G_SetMeasured, [c_void_p, c_int], errcheck=False)
 
 # Geometry modification routines.
 add_geom = void_output(lgdal.OGR_G_AddGeometry, [c_void_p, c_void_p])
@@ -94,9 +97,13 @@ destroy_geom = void_output(lgdal.OGR_G_DestroyGeometry, [c_void_p], errcheck=Fal
 to_wkb = void_output(
     lgdal.OGR_G_ExportToWkb, None, errcheck=True
 )  # special handling for WKB.
+to_iso_wkb = void_output(lgdal.OGR_G_ExportToIsoWkb, None, errcheck=True)
 to_wkt = string_output(
     lgdal.OGR_G_ExportToWkt, [c_void_p, POINTER(c_char_p)], decoding="ascii"
 )
+to_iso_wkt = string_output(
+    lgdal.OGR_G_ExportToIsoWkt, [c_void_p, POINTER(c_char_p)], decoding="ascii"
+)
 to_gml = string_output(
     lgdal.OGR_G_ExportToGML, [c_void_p], str_result=True, decoding="ascii"
 )
@@ -115,7 +122,7 @@ get_geom_srs = srs_output(lgdal.OGR_G_GetSpatialReference, [c_void_p])
 get_area = double_output(lgdal.OGR_G_GetArea, [c_void_p])
 get_centroid = void_output(lgdal.OGR_G_Centroid, [c_void_p, c_void_p])
 get_dims = int_output(lgdal.OGR_G_GetDimension, [c_void_p])
-get_coord_dim = int_output(lgdal.OGR_G_GetCoordinateDimension, [c_void_p])
+get_coord_dim = int_output(lgdal.OGR_G_CoordinateDimension, [c_void_p])
 set_coord_dim = void_output(
     lgdal.OGR_G_SetCoordinateDimension, [c_void_p, c_int], errcheck=False
 )
diff --git a/django/contrib/gis/utils/layermapping.py b/django/contrib/gis/utils/layermapping.py
index 3e33c27772..d1fe7c5721 100644
--- a/django/contrib/gis/utils/layermapping.py
+++ b/django/contrib/gis/utils/layermapping.py
@@ -488,6 +488,10 @@ class LayerMapping:
         if necessary (for example if the model field is MultiPolygonField while
         the mapped shapefile only contains Polygons).
         """
+        # Measured geometries are not yet supported by GeoDjango models.
+        if geom.is_measured:
+            geom.set_measured(False)
+
         # Downgrade a 3D geom to a 2D one, if necessary.
         if self.coord_dim == 2 and geom.is_3d:
             geom.set_3d(False)
diff --git a/docs/ref/contrib/gis/gdal.txt b/docs/ref/contrib/gis/gdal.txt
index 4b7a706fc4..f6e9ef0964 100644
--- a/docs/ref/contrib/gis/gdal.txt
+++ b/docs/ref/contrib/gis/gdal.txt
@@ -581,6 +581,27 @@ coordinate transformation:
         >>> p.wkt
         "POINT (1 2)"
 
+    .. attribute:: is_measured
+
+    .. versionadded:: 5.1
+
+    A boolean indicating if this geometry has M coordinates.
+
+    .. method:: set_measured(value)
+
+    .. versionadded:: 5.1
+
+    A method to add or remove the M coordinate dimension.
+
+    .. code-block:: pycon
+
+        >>> p = OGRGeometry("POINT (1 2)")
+        >>> p.is_measured
+        False
+        >>> p.set_measured(True)
+        >>> p.wkt
+        "POINT M (1 2 0)"
+
     .. attribute:: geom_count
 
     Returns the number of elements in this geometry:
@@ -864,6 +885,18 @@ coordinate transformation:
         >>> OGRGeometry("POINT (1 2 3)").z
         3.0
 
+    .. attribute:: m
+
+    .. versionadded:: 5.1
+
+    Returns the M coordinate of this point, or ``None`` if the Point does not
+    have an M coordinate:
+
+    .. code-block:: pycon
+
+        >>> OGRGeometry("POINT ZM (1 2 3 4)").m
+        4.0
+
 .. class:: LineString
 
     .. attribute:: x
diff --git a/docs/releases/5.1.txt b/docs/releases/5.1.txt
index 9e8ff8c049..463551e85c 100644
--- a/docs/releases/5.1.txt
+++ b/docs/releases/5.1.txt
@@ -78,6 +78,11 @@ Minor features
 * The new :meth:`.OGRGeometry.set_3d` method allows addition and removal of the
   ``Z`` coordinate dimension.
 
+* :class:`~django.contrib.gis.gdal.OGRGeometry` and
+  :class:`~django.contrib.gis.gdal.Point` now support measured geometries
+  via the new :attr:`.OGRGeometry.is_measured` and :attr:`.Point.m` properties,
+  and the :meth:`.OGRGeometry.set_measured` method.
+
 :mod:`django.contrib.messages`
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
diff --git a/tests/gis_tests/gdal_tests/test_geom.py b/tests/gis_tests/gdal_tests/test_geom.py
index 13e7d3e70d..a2aa9d8e36 100644
--- a/tests/gis_tests/gdal_tests/test_geom.py
+++ b/tests/gis_tests/gdal_tests/test_geom.py
@@ -672,7 +672,7 @@ class OGRGeomTest(SimpleTestCase, TestDataMixin):
             ("PolyhedralSurface Z", 1015, False),
             ("TIN Z", 1016, False),
             ("Triangle Z", 1017, False),
-            ("Point M", 2001, False),
+            ("Point M", 2001, True),
             ("LineString M", 2002, False),
             ("Polygon M", 2003, False),
             ("MultiPoint M", 2004, False),
@@ -687,7 +687,7 @@ class OGRGeomTest(SimpleTestCase, TestDataMixin):
             ("PolyhedralSurface M", 2015, False),
             ("TIN M", 2016, False),
             ("Triangle M", 2017, False),
-            ("Point ZM", 3001, False),
+            ("Point ZM", 3001, True),
             ("LineString ZM", 3002, False),
             ("Polygon ZM", 3003, False),
             ("MultiPoint ZM", 3004, False),
@@ -812,6 +812,71 @@ class OGRGeomTest(SimpleTestCase, TestDataMixin):
                 self.assertEqual(g.wkt, geom)
                 self.assertEqual(g.wkb.hex(), wkb)
 
+    def test_measure_is_measure_and_set_measure(self):
+        geom = OGRGeometry("POINT (1 2 3)")
+        self.assertIs(geom.is_measured, False)
+        geom.set_measured(True)
+        self.assertIs(geom.is_measured, True)
+        self.assertEqual(geom.wkt, "POINT ZM (1 2 3 0)")
+        geom.set_measured(False)
+        self.assertIs(geom.is_measured, False)
+        self.assertEqual(geom.wkt, "POINT (1 2 3)")
+        msg = "Input to 'set_measured' must be a boolean, got 'None'"
+        with self.assertRaisesMessage(ValueError, msg):
+            geom.set_measured(None)
+
+    def test_point_m_coordinate(self):
+        geom = OGRGeometry("POINT ZM (1 2 3 4)")
+        self.assertEqual(geom.m, 4)
+        geom = OGRGeometry("POINT (1 2 3 4)")
+        self.assertEqual(geom.m, 4)
+        geom = OGRGeometry("POINT M (1 2 3)")
+        self.assertEqual(geom.m, 3)
+        geom = OGRGeometry("POINT Z (1 2 3)")
+        self.assertEqual(geom.m, None)
+
+    def test_point_m_tuple(self):
+        geom = OGRGeometry("POINT ZM (1 2 3 4)")
+        self.assertEqual(geom.tuple, (geom.x, geom.y, geom.z, geom.m))
+        geom = OGRGeometry("POINT M (1 2 3)")
+        self.assertEqual(geom.tuple, (geom.x, geom.y, geom.m))
+        geom = OGRGeometry("POINT Z (1 2 3)")
+        self.assertEqual(geom.tuple, (geom.x, geom.y, geom.z))
+        geom = OGRGeometry("POINT (1 2 3)")
+        self.assertEqual(geom.tuple, (geom.x, geom.y, geom.z))
+
+    def test_point_m_wkt_wkb(self):
+        wkt = "POINT ZM (1 2 3 4)"
+        geom = OGRGeometry(wkt)
+        self.assertEqual(geom.wkt, wkt)
+        self.assertEqual(
+            geom.wkb.hex(),
+            "01b90b0000000000000000f03f00000000000000"
+            "4000000000000008400000000000001040",
+        )
+        wkt = "POINT M (1 2 3)"
+        geom = OGRGeometry(wkt)
+        self.assertEqual(geom.wkt, wkt)
+        self.assertEqual(
+            geom.wkb.hex(),
+            "01d1070000000000000000f03f00000000000000400000000000000840",
+        )
+
+    def test_point_m_dimension_types(self):
+        geom = OGRGeometry("POINT ZM (1 2 3 4)")
+        self.assertEqual(geom.geom_type.name, "PointZM")
+        self.assertEqual(geom.geom_type.num, 3001)
+        geom = OGRGeometry("POINT M (1 2 3)")
+        self.assertEqual(geom.geom_type.name, "PointM")
+        self.assertEqual(geom.geom_type.num, 2001)
+
+    def test_point_m_dimension_geos(self):
+        """GEOSGeometry does not yet support the M dimension."""
+        geom = OGRGeometry("POINT ZM (1 2 3 4)")
+        self.assertEqual(geom.geos.wkt, "POINT Z (1 2 3)")
+        geom = OGRGeometry("POINT M (1 2 3)")
+        self.assertEqual(geom.geos.wkt, "POINT (1 2)")
+
 
 class DeprecationTests(SimpleTestCase):
     def test_coord_setter_deprecation(self):