ENH: get_m and point.m property

docs/manual.rst

@@ -302,7 +302,7 @@ Its `x-y` bounding box is a ``(minx, miny, maxx, maxy)`` tuple.
   >>> point.bounds
   (0.0, 0.0, 0.0, 0.0)
 
-Coordinate values are accessed via `coords`, `x`, `y`, and `z` properties.
+Coordinate values are accessed via `coords`, `x`, `y`, `z`, and `m` properties.
 
 .. code-block:: pycon
 

shapely/_geometry.py

@@ -5,7 +5,7 @@
 
 from shapely import _geometry_helpers, geos_version, lib
 from shapely._enum import ParamEnum
-from shapely.decorators import multithreading_enabled
+from shapely.decorators import multithreading_enabled, requires_geos
 
 __all__ = [
     "GeometryType",
@@ -18,6 +18,7 @@
     "get_x",
     "get_y",
     "get_z",
+    "get_m",
     "get_exterior_ring",
     "get_num_points",
     "get_num_interior_rings",
@@ -257,7 +258,7 @@ def get_x(point, **kwargs):
 
     See also
     --------
-    get_y, get_z
+    get_y, get_z, get_m
 
     Examples
     --------
@@ -283,7 +284,7 @@ def get_y(point, **kwargs):
 
     See also
     --------
-    get_x, get_z
+    get_x, get_z, get_m
 
     Examples
     --------
@@ -303,14 +304,14 @@ def get_z(point, **kwargs):
     Parameters
     ----------
     point : Geometry or array_like
-        Non-point geometries or geometries without 3rd dimension will result
+        Non-point geometries or geometries without Z dimension will result
         in NaN being returned.
     **kwargs
         See :ref:`NumPy ufunc docs <ufuncs.kwargs>` for other keyword arguments.
 
     See also
     --------
-    get_x, get_y
+    get_x, get_y, get_m
 
     Examples
     --------
@@ -325,6 +326,40 @@ def get_z(point, **kwargs):
     return lib.get_z(point, **kwargs)
 
 
+@multithreading_enabled
+@requires_geos("3.12.0")
+def get_m(point, **kwargs):
+    """Returns the m-coordinate of a point.
+
+    .. versionadded:: 2.1.0
+
+    Parameters
+    ----------
+    point : Geometry or array_like
+        Non-point geometries or geometries without M dimension will result
+        in NaN being returned.
+    **kwargs
+        See :ref:`NumPy ufunc docs <ufuncs.kwargs>` for other keyword arguments.
+
+    See also
+    --------
+    get_x, get_y, get_z
+
+    Examples
+    --------
+    >>> from shapely import Point, from_wkt
+    >>> get_m(from_wkt("POINT ZM (1 2 3 4)"))
+    4.0
+    >>> get_m(from_wkt("POINT M (1 2 4)"))
+    4.0
+    >>> get_m(Point(1, 2, 3))
+    nan
+    >>> get_m(from_wkt("MULTIPOINT M ((1 1 1), (2 2 2))"))
+    nan
+    """
+    return lib.get_m(point, **kwargs)
+
+
 # linestrings
 
 

shapely/geometry/point.py

@@ -3,6 +3,7 @@
 import numpy as np
 
 import shapely
+from shapely import lib
 from shapely.errors import DimensionError
 from shapely.geometry.base import BaseGeometry
 
@@ -99,6 +100,15 @@ def z(self):
             raise DimensionError("This point has no z coordinate.")
         return shapely.get_z(self)
 
+    if lib.geos_version >= (3, 12, 0):
+
+        @property
+        def m(self):
+            """Return m coordinate."""
+            if not shapely.has_m(self):
+                raise DimensionError("This point has no m coordinate.")
+            return shapely.get_m(self)
+
     @property
     def __geo_interface__(self):
         return {"type": "Point", "coordinates": self.coords[0]}

shapely/tests/geometry/test_point.py

@@ -1,7 +1,7 @@
 import numpy as np
 import pytest
 
-from shapely import Point
+from shapely import geos_version, Point
 from shapely.coords import CoordinateSequence
 from shapely.errors import DimensionError
 
@@ -98,7 +98,7 @@ def test_from_invalid():
 class TestPoint:
     def test_point(self):
 
-        # Test 2D points
+        # Test XY point
         p = Point(1.0, 2.0)
         assert p.x == 1.0
         assert p.y == 2.0
@@ -107,13 +107,25 @@ def test_point(self):
         assert p.has_z is False
         with pytest.raises(DimensionError):
             p.z
-
-        # Check Z-dim
+        if geos_version >= (3, 12, 0):
+            assert p.has_m is False
+            with pytest.raises(DimensionError):
+                p.m
+        else:
+            assert not hasattr(p, "m")
+
+        # Check XYZ point
         p = Point(1.0, 2.0, 3.0)
         assert p.coords[:] == [(1.0, 2.0, 3.0)]
         assert str(p) == p.wkt
         assert p.has_z is True
         assert p.z == 3.0
+        if geos_version >= (3, 12, 0):
+            assert p.has_m is False
+            with pytest.raises(DimensionError):
+                p.m
+
+            # TODO: Check XYM and XYZM points
 
         # Coordinate access
         p = Point((3.0, 4.0))

shapely/tests/test_geometry.py

@@ -30,7 +30,9 @@
     multi_polygon,
     multi_polygon_z,
     point,
+    point_m,
     point_z,
+    point_zm,
     polygon,
     polygon_with_hole,
     polygon_with_hole_z,
@@ -177,6 +179,12 @@ def test_get_set_srid():
         shapely.get_x,
         shapely.get_y,
         shapely.get_z,
+        pytest.param(
+            shapely.get_m,
+            marks=pytest.mark.skipif(
+                shapely.geos_version < (3, 12, 0), reason="GEOS < 3.12"
+            ),
+        ),
     ],
 )
 @pytest.mark.parametrize(
@@ -203,6 +211,16 @@ def test_get_z_2d():
     assert np.isnan(shapely.get_z(point))
 
 
+@pytest.mark.skipif(
+    shapely.geos_version < (3, 12, 0),
+    reason="M coordinates not supported with GEOS < 3.12",
+)
+def test_get_m():
+    assert shapely.get_m([point_m, point_zm]).tolist() == [5.0, 5.0]
+    assert np.isnan(shapely.get_m(point))
+    assert np.isnan(shapely.get_m(point_z))
+
+
 @pytest.mark.parametrize("geom", all_types)
 def test_new_from_wkt(geom):
     if geom.is_empty and shapely.get_num_geometries(geom) > 0:

src/ufuncs.c

@@ -1055,6 +1055,18 @@ static int GetZ(void* context, void* a, double* b) {
   }
 }
 static void* get_z_data[1] = {GetZ};
+#if GEOS_SINCE_3_12_0
+static int GetM(void* context, void* a, double* b) {
+  char typ = GEOSGeomTypeId_r(context, a);
+  if (typ != 0) {
+    *(double*)b = NPY_NAN;
+    return 1;
+  } else {
+    return GEOSGeomGetM_r(context, a, b);
+  }
+}
+static void* get_m_data[1] = {GetM};
+#endif
 static void* area_data[1] = {GEOSArea_r};
 static void* length_data[1] = {GEOSLength_r};
 
@@ -3837,6 +3849,7 @@ int init_ufuncs(PyObject* m, PyObject* d) {
 
 #if GEOS_SINCE_3_12_0
   DEFINE_Y_b(has_m);
+  DEFINE_Y_d(get_m);
 #endif
 
   Py_DECREF(ufunc);