Repository: shapely/shapely
Pull request #2235 (2025-02-23 23:38:03)

ENH: Add `include_m` to `to_ragged_array`
This enhancement allows [`to_ragged_array()`](https://shapely.readthedocs.io/en/2.0.7/reference/shapely.to_ragged_array.html) to get M coordinate data with `include_m`, which behaves similar to `include_z`, where it can be None (default), False or True. Setting `include_m=True` with GEOS<3.12.0 will get all NaN coordinate data, as these GEOS versions don't support geometries with M dimensions.

The tests `test_include_z` and `test_include_z_false` are replaced with `test_to_ragged_array(geom, include_z, include_m)` with all combinations of parameters, including XYM and XYZM geometries for GEOS >= 3.12.0.

A minor unrelated fix is to correct a typo: `stringt` -> `string`.

Next steps:

- [ ] Update `from_ragged_array()` to support XYM and XYZM geometries

coveralls: ## Pull Request Test Coverage Report for [Build 13488006998](https://coveralls.io/builds/72384707)


### Details

* **31** of **31**   **(100.0%)**  changed or added relevant lines in **1** file are covered.
* No unchanged relevant lines lost coverage.
* Overall coverage increased (+**0.3%**) to **87.844%**

---



|  Totals | [![Coverage Status](https://coveralls.io/builds/72384707/badge)](https://coveralls.io/builds/72384707) |
| :-- | --: |
| Change from base [Build 13479085130](https://coveralls.io/builds/72379212): |  0.3% |
| Covered Lines: | 2652 |
| Relevant Lines: | 3019 |

---
##### 💛  - [Coveralls](https://coveralls.io)


Patch:
diff --git a/docs/manual.rst b/docs/manual.rst
index 09a4353..02e6f3e 100644
--- a/docs/manual.rst
+++ b/docs/manual.rst
@@ -350,7 +350,8 @@ Figure 1. A simple `LineString` on the left, a complex `LineString` on the
 right. The (`MultiPoint`) boundary of each is shown in black, the other points
 that describe the lines are shown in grey.
 
-A `LineString` has zero area and non-zero length.
+A `LineString` has zero area and non-zero length if it isn't
+:ref:`empty <empties>`.
 
 .. code-block:: pycon
 
@@ -430,7 +431,8 @@ grey. A ring's boundary is `empty`.
    Shapely will not prevent the creation of such rings, but exceptions will be
    raised when they are operated on.
 
-A `LinearRing` has zero area and non-zero length.
+A `LinearRing` has zero area and non-zero length if it isn't
+:ref:`empty <empties>`.
 
 .. code-block:: pycon
 
@@ -466,7 +468,7 @@ instance, respectively.
   >>> LinearRing(ring)
   <LINEARRING (0 0, 1 1, 1 0, 0 0)>
 
-As with `LineString`, a sequence of `Point` instances is not a valid
+As with `LineString`, a sequence of `Point` instances is a valid
 constructor parameter.
 
 .. _polygons:
@@ -484,7 +486,8 @@ Polygons
 
 Rings of a `valid` `Polygon` may not cross each other, but may touch at a
 single point only.  Again, Shapely will not prevent the creation of invalid
-features, but exceptions will be raised when they are operated on.
+features, but when they are operated on the results might be wrong or
+exceptions might be raised.
 
 .. plot:: code/polygon.py
 
@@ -499,7 +502,8 @@ Figure 4. On the left, a `Polygon` that is `invalid` because its exterior and
 interior rings touch along a line, and on the right, a `Polygon` that is
 `invalid` because its interior rings touch along a line.
 
-A `Polygon` has non-zero area and non-zero length.
+A `Polygon` has non-zero area and non-zero length if it isn't
+:ref:`empty <empties>`.
 
 .. code-block:: pycon
 
@@ -676,7 +680,8 @@ Figure 6. On the left, a `simple`, disconnected `MultiLineString`, and on the
 right, a non-simple `MultiLineString`. The points defining the objects are
 shown in gray, the boundaries of the objects in black.
 
-A `MultiLineString` has zero area and non-zero length.
+A `MultiLineString` has zero area and non-zero length if it isn't
+:ref:`empty <empties>`.
 
 .. code-block:: pycon
 
@@ -1028,7 +1033,8 @@ Operations on non-simple `LineStrings` are fully supported by Shapely.
 
 A valid `Polygon` may not possess any overlapping exterior or interior rings. A
 valid `MultiPolygon` may not collect any overlapping polygons. A valid `LineString`
-must have non-zero length. Operations on invalid features may fail.
+must have non-zero length if it isn't :ref:`empty <empties>`. Operations on invalid
+features may fail.
 
 .. code-block:: pycon
 
diff --git a/shapely/_ragged_array.py b/shapely/_ragged_array.py
index c813120..9438f4b 100644
--- a/shapely/_ragged_array.py
+++ b/shapely/_ragged_array.py
@@ -26,10 +26,9 @@ https://github.com/geoarrow/geoarrow
 
 import numpy as np
 
-from shapely import creation
+from shapely import creation, geos_version
 from shapely._geometry import (
     GeometryType,
-    get_coordinate_dimension,
     get_parts,
     get_rings,
     get_type_id,
@@ -43,13 +42,15 @@ from shapely.predicates import is_empty, is_missing
 
 __all__ = ["from_ragged_array", "to_ragged_array"]
 
+_geos_ge_312 = geos_version >= (3, 12, 0)
+
 
 # # GEOS -> coords/offset arrays (to_ragged_array)
 
 
-def _get_arrays_point(arr, include_z):
+def _get_arrays_point(arr, include_z, include_m):
     # only one array of coordinates
-    coords = get_coordinates(arr, include_z=include_z)
+    coords = get_coordinates(arr, include_z=include_z, include_m=include_m)
 
     # empty points are represented by NaNs
     # + missing geometries should also be present with some value
@@ -83,53 +84,59 @@ def _indices_to_offsets(indices, n):
     return offsets
 
 
-def _get_arrays_multipoint(arr, include_z):
+def _get_arrays_multipoint(arr, include_z, include_m):
     # explode/flatten the MultiPoints
     _, part_indices = get_parts(arr, return_index=True)
     # the offsets into the multipoint parts
     offsets = _indices_to_offsets(part_indices, len(arr))
 
     # only one array of coordinates
-    coords = get_coordinates(arr, include_z=include_z)
+    coords = get_coordinates(arr, include_z=include_z, include_m=include_m)
 
     return coords, (offsets,)
 
 
-def _get_arrays_linestring(arr, include_z):
+def _get_arrays_linestring(arr, include_z, include_m):
     # the coords and offsets into the coordinates of the linestrings
-    coords, indices = get_coordinates(arr, return_index=True, include_z=include_z)
+    coords, indices = get_coordinates(
+        arr, return_index=True, include_z=include_z, include_m=include_m
+    )
     offsets = _indices_to_offsets(indices, len(arr))
 
     return coords, (offsets,)
 
 
-def _get_arrays_multilinestring(arr, include_z):
+def _get_arrays_multilinestring(arr, include_z, include_m):
     # explode/flatten the MultiLineStrings
     arr_flat, part_indices = get_parts(arr, return_index=True)
     # the offsets into the multilinestring parts
     offsets2 = _indices_to_offsets(part_indices, len(arr))
 
     # the coords and offsets into the coordinates of the linestrings
-    coords, indices = get_coordinates(arr_flat, return_index=True, include_z=include_z)
+    coords, indices = get_coordinates(
+        arr_flat, return_index=True, include_z=include_z, include_m=include_m
+    )
     offsets1 = _indices_to_offsets(indices, len(arr_flat))
 
     return coords, (offsets1, offsets2)
 
 
-def _get_arrays_polygon(arr, include_z):
+def _get_arrays_polygon(arr, include_z, include_m):
     # explode/flatten the Polygons into Rings
     arr_flat, ring_indices = get_rings(arr, return_index=True)
     # the offsets into the exterior/interior rings of the multipolygon parts
     offsets2 = _indices_to_offsets(ring_indices, len(arr))
 
     # the coords and offsets into the coordinates of the rings
-    coords, indices = get_coordinates(arr_flat, return_index=True, include_z=include_z)
+    coords, indices = get_coordinates(
+        arr_flat, return_index=True, include_z=include_z, include_m=include_m
+    )
     offsets1 = _indices_to_offsets(indices, len(arr_flat))
 
     return coords, (offsets1, offsets2)
 
 
-def _get_arrays_multipolygon(arr, include_z):
+def _get_arrays_multipolygon(arr, include_z, include_m):
     # explode/flatten the MultiPolygons
     arr_flat, part_indices = get_parts(arr, return_index=True)
     # the offsets into the multipolygon parts
@@ -141,13 +148,15 @@ def _get_arrays_multipolygon(arr, include_z):
     offsets2 = _indices_to_offsets(ring_indices, len(arr_flat))
 
     # the coords and offsets into the coordinates of the rings
-    coords, indices = get_coordinates(arr_flat2, return_index=True, include_z=include_z)
+    coords, indices = get_coordinates(
+        arr_flat2, return_index=True, include_z=include_z, include_m=include_m
+    )
     offsets1 = _indices_to_offsets(indices, len(arr_flat2))
 
     return coords, (offsets1, offsets2, offsets3)
 
 
-def to_ragged_array(geometries, include_z=None):
+def to_ragged_array(geometries, include_z=None, include_m=None):
     """Convert geometries to a ragged array representation.
 
     This function converts an array of geometries to a ragged array
@@ -164,14 +173,21 @@ def to_ragged_array(geometries, include_z=None):
     ----------
     geometries : array_like
         Array of geometries (1-dimensional).
-    include_z : bool, default None
-        If False, return 2D geometries. If True, include the third dimension
-        in the output (if a geometry has no third dimension, the z-coordinates
-        will be NaN). By default, will infer the dimensionality from the
+    include_z, include_m : bool, default None
+        If both are False, return XY (2D) geometries.
+        If both are True, return XYZM (4D) geometries.
+        If either is True, return either XYZ or XYM (3D) geometries.
+        If a geometry has no Z or M dimension, extra coordinate data will be NaN.
+        By default, will infer the dimensionality from the
         input geometries. Note that this inference can be unreliable with
         empty geometries (for a guaranteed result, it is recommended to
         specify the keyword).
 
+        .. versionadded:: 2.1.0
+            The ``include_m`` parameter was added to support XYM (3D) and
+            XYZM (4D) geometries available with GEOS 3.12.0 or later.
+            With older GEOS versions, M dimension coordinates will be NaN.
+
     Returns
     -------
     tuple of (geometry_type, coords, offsets)
@@ -179,8 +195,8 @@ def to_ragged_array(geometries, include_z=None):
             The type of the input geometries (required information for
             roundtrip).
         coords : np.ndarray
-            Contiguous array of shape (n, 2) or (n, 3) of all coordinates
-            of all input geometries.
+            Contiguous array of shape (n, 2), (n, 3), or (n, 4) of all
+            coordinates of all input geometries.
         offsets: tuple of np.ndarray
             Offset arrays that make it possible to reconstruct the
             geometries from the flat coordinates array. The number of
@@ -257,43 +273,49 @@ def to_ragged_array(geometries, include_z=None):
            [ 0.,  0.]])
 
     """
+    from shapely import has_m, has_z  # avoid circular import
+
     geometries = np.asarray(geometries)
     if include_z is None:
-        include_z = np.any(
-            get_coordinate_dimension(geometries[~is_empty(geometries)]) == 3
-        )
+        include_z = np.any(has_z(geometries[~is_empty(geometries)]))
+    if include_m is None:
+        if _geos_ge_312:
+            include_m = np.any(has_m(geometries[~is_empty(geometries)]))
+        else:
+            include_m = False
 
     geom_types = np.unique(get_type_id(geometries))
     # ignore missing values (type of -1)
     geom_types = geom_types[geom_types >= 0]
 
+    get_arrays_args = geometries, include_z, include_m
     if len(geom_types) == 1:
         typ = GeometryType(geom_types[0])
         if typ == GeometryType.POINT:
-            coords, offsets = _get_arrays_point(geometries, include_z)
+            coords, offsets = _get_arrays_point(*get_arrays_args)
         elif typ == GeometryType.LINESTRING:
-            coords, offsets = _get_arrays_linestring(geometries, include_z)
+            coords, offsets = _get_arrays_linestring(*get_arrays_args)
         elif typ == GeometryType.POLYGON:
-            coords, offsets = _get_arrays_polygon(geometries, include_z)
+            coords, offsets = _get_arrays_polygon(*get_arrays_args)
         elif typ == GeometryType.MULTIPOINT:
-            coords, offsets = _get_arrays_multipoint(geometries, include_z)
+            coords, offsets = _get_arrays_multipoint(*get_arrays_args)
         elif typ == GeometryType.MULTILINESTRING:
-            coords, offsets = _get_arrays_multilinestring(geometries, include_z)
+            coords, offsets = _get_arrays_multilinestring(*get_arrays_args)
         elif typ == GeometryType.MULTIPOLYGON:
-            coords, offsets = _get_arrays_multipolygon(geometries, include_z)
+            coords, offsets = _get_arrays_multipolygon(*get_arrays_args)
         else:
             raise ValueError(f"Geometry type {typ.name} is not supported")
 
     elif len(geom_types) == 2:
         if set(geom_types) == {GeometryType.POINT, GeometryType.MULTIPOINT}:
             typ = GeometryType.MULTIPOINT
-            coords, offsets = _get_arrays_multipoint(geometries, include_z)
+            coords, offsets = _get_arrays_multipoint(*get_arrays_args)
         elif set(geom_types) == {GeometryType.LINESTRING, GeometryType.MULTILINESTRING}:
             typ = GeometryType.MULTILINESTRING
-            coords, offsets = _get_arrays_multilinestring(geometries, include_z)
+            coords, offsets = _get_arrays_multilinestring(*get_arrays_args)
         elif set(geom_types) == {GeometryType.POLYGON, GeometryType.MULTIPOLYGON}:
             typ = GeometryType.MULTIPOLYGON
-            coords, offsets = _get_arrays_multipolygon(geometries, include_z)
+            coords, offsets = _get_arrays_multipolygon(*get_arrays_args)
         else:
             raise ValueError(
                 "Geometry type combination is not supported "


Test patch:
diff --git a/shapely/tests/common.py b/shapely/tests/common.py
index fd06191..ca24b55 100644
--- a/shapely/tests/common.py
+++ b/shapely/tests/common.py
@@ -37,7 +37,7 @@ empty_multi_point = shapely.from_wkt("MULTIPOINT EMPTY")
 empty_multi_line_string = shapely.from_wkt("MULTILINESTRING EMPTY")
 empty_multi_polygon = shapely.from_wkt("MULTIPOLYGON EMPTY")
 multi_point_empty = shapely.multipoints([empty_point])
-multi_line_stringt_empty = shapely.multilinestrings([empty_line_string])
+multi_line_string_empty = shapely.multilinestrings([empty_line_string])
 multi_polygon_empty = shapely.multipolygons([empty_polygon])
 geometry_collection_empty = shapely.geometrycollections([empty_line_string])
 # XYZ
@@ -68,7 +68,7 @@ empty_multi_point_z = shapely.from_wkt("MULTIPOINT Z EMPTY")
 empty_multi_line_string_z = shapely.from_wkt("MULTILINESTRING Z EMPTY")
 empty_multi_polygon_z = shapely.from_wkt("MULTIPOLYGON Z EMPTY")
 multi_point_empty_z = shapely.multipoints([empty_point_z])
-multi_line_stringt_empty_z = shapely.multilinestrings([empty_line_string_z])
+multi_line_string_empty_z = shapely.multilinestrings([empty_line_string_z])
 multi_polygon_empty_z = shapely.multipolygons([empty_polygon_z])
 geometry_collection_empty_z = shapely.geometrycollections([empty_line_string_z])
 # XYM
@@ -95,7 +95,7 @@ empty_multi_point_m = shapely.from_wkt("MULTIPOINT M EMPTY")
 empty_multi_line_string_m = shapely.from_wkt("MULTILINESTRING M EMPTY")
 empty_multi_polygon_m = shapely.from_wkt("MULTIPOLYGON M EMPTY")
 multi_point_empty_m = shapely.multipoints([empty_point_m])
-multi_line_stringt_empty_m = shapely.multilinestrings([empty_line_string_m])
+multi_line_string_empty_m = shapely.multilinestrings([empty_line_string_m])
 multi_polygon_empty_m = shapely.multipolygons([empty_polygon_m])
 geometry_collection_empty_m = shapely.geometrycollections([empty_line_string_m])
 # XYZM
@@ -126,7 +126,7 @@ empty_multi_point_zm = shapely.from_wkt("MULTIPOINT ZM EMPTY")
 empty_multi_line_string_zm = shapely.from_wkt("MULTILINESTRING ZM EMPTY")
 empty_multi_polygon_zm = shapely.from_wkt("MULTIPOLYGON ZM EMPTY")
 multi_point_empty_zm = shapely.multipoints([empty_point_zm])
-multi_line_stringt_empty_zm = shapely.multilinestrings([empty_line_string_zm])
+multi_line_string_empty_zm = shapely.multilinestrings([empty_line_string_zm])
 multi_polygon_empty_zm = shapely.multipolygons([empty_polygon_zm])
 geometry_collection_empty_zm = shapely.geometrycollections([empty_line_string_zm])
 
@@ -148,7 +148,7 @@ all_types = (
     empty_multi_line_string,
     empty_multi_polygon,
     multi_point_empty,
-    multi_line_stringt_empty,
+    multi_line_string_empty,
     multi_polygon_empty,
     geometry_collection_empty,
 )
@@ -171,7 +171,7 @@ all_types_z = (
     empty_multi_line_string_z,
     empty_multi_polygon_z,
     multi_point_empty_z,
-    multi_line_stringt_empty_z,
+    multi_line_string_empty_z,
     multi_polygon_empty_z,
     geometry_collection_empty_z,
 )
@@ -193,7 +193,7 @@ all_types_m = (
     empty_multi_line_string_m,
     empty_multi_polygon_m,
     multi_point_empty_m,
-    multi_line_stringt_empty_m,
+    multi_line_string_empty_m,
     multi_polygon_empty_m,
     geometry_collection_empty_m,
 )
@@ -215,7 +215,7 @@ all_types_zm = (
     empty_multi_line_string_zm,
     empty_multi_polygon_zm,
     multi_point_empty_zm,
-    multi_line_stringt_empty_zm,
+    multi_line_string_empty_zm,
     multi_polygon_empty_zm,
     geometry_collection_empty_zm,
 )
diff --git a/shapely/tests/test_ragged_array.py b/shapely/tests/test_ragged_array.py
index 023b16c..0e4d5df 100644
--- a/shapely/tests/test_ragged_array.py
+++ b/shapely/tests/test_ragged_array.py
@@ -7,21 +7,38 @@ from shapely import MultiLineString, MultiPoint, MultiPolygon
 from shapely.testing import assert_geometries_equal
 from shapely.tests.common import (
     empty_line_string,
+    empty_line_string_m,
     empty_line_string_z,
+    empty_line_string_zm,
+    empty_multi_polygon_m,
+    empty_multi_polygon_z,
+    empty_multi_polygon_zm,
     geometry_collection,
     line_string,
+    line_string_m,
     line_string_z,
+    line_string_zm,
     linear_ring,
     multi_line_string,
+    multi_line_string_m,
     multi_line_string_z,
+    multi_line_string_zm,
     multi_point,
+    multi_point_m,
     multi_point_z,
+    multi_point_zm,
     multi_polygon,
+    multi_polygon_m,
     multi_polygon_z,
+    multi_polygon_zm,
     point,
+    point_m,
     point_z,
+    point_zm,
     polygon,
+    polygon_m,
     polygon_z,
+    polygon_zm,
 )
 
 all_types = (
@@ -42,6 +59,28 @@ all_types_z = (
     multi_polygon_z,
 )
 
+all_types_m = (
+    point_m,
+    line_string_m,
+    polygon_m,
+    multi_point_m,
+    multi_line_string_m,
+    multi_polygon_m,
+)
+
+all_types_zm = (
+    point_zm,
+    line_string_zm,
+    polygon_zm,
+    multi_point_zm,
+    multi_line_string_zm,
+    multi_polygon_zm,
+)
+
+all_types_dims_combos = all_types + all_types_z
+if shapely.geos_version >= (3, 12, 0):
+    all_types_dims_combos = all_types_dims_combos + all_types_m + all_types_zm
+
 all_types_not_supported = (
     linear_ring,
     geometry_collection,
@@ -54,24 +93,31 @@ def test_roundtrip(geom):
     assert_geometries_equal(actual, [geom, geom])
 
 
-@pytest.mark.parametrize("geom", all_types)
-def test_include_z(geom):
-    _, coords, _ = shapely.to_ragged_array([geom, geom], include_z=True)
-    # For 2D geoms, z coords are filled in with NaN
-    assert np.isnan(coords[:, 2]).all()
-
-
-@pytest.mark.parametrize("geom", all_types_z)
-def test_include_z_false(geom):
-    _, coords, _ = shapely.to_ragged_array([geom, geom], include_z=False)
-    # For 3D geoms, z coords are dropped
-    assert coords.shape[1] == 2
+@pytest.mark.parametrize("include_m", [None, True, False])
+@pytest.mark.parametrize("include_z", [None, True, False])
+@pytest.mark.parametrize("geom", all_types_dims_combos)
+def test_to_ragged_array(geom, include_z, include_m):
+    _, coords, _ = shapely.to_ragged_array(
+        [geom, geom], include_z=include_z, include_m=include_m
+    )
+    nan_dims = np.all(np.isnan(coords), axis=0).tolist()
+    expected = [False, False]  # XY
+    has_z = geom.has_z
+    if include_z or (include_z is None and has_z):
+        expected.append(not has_z)  # XYZ or XYZM
+    if shapely.geos_version >= (3, 12, 0):
+        has_m = geom.has_m
+    else:
+        has_m = False
+    if include_m or (include_m is None and has_m):
+        expected.append(not has_m)  # XYM or XYZM
+    assert nan_dims == expected
 
 
 def test_include_z_default():
     # corner cases for inferring dimensionality
 
-    # mixed 2D and 3D -> 3D
+    # mixed XY and XYZ -> XYZ
     _, coords, _ = shapely.to_ragged_array([line_string, line_string_z])
     assert coords.shape[1] == 3
 
@@ -81,7 +127,31 @@ def test_include_z_default():
     _, coords, _ = shapely.to_ragged_array([empty_line_string_z])
     assert coords.shape[1] == 2
     # empty collection -> GEOS indicates 2D
-    _, coords, _ = shapely.to_ragged_array(shapely.from_wkt(["MULTIPOLYGON Z EMPTY"]))
+    _, coords, _ = shapely.to_ragged_array([empty_multi_polygon_z])
+    assert coords.shape[1] == 2
+
+
+@pytest.mark.skipif(shapely.geos_version < (3, 12, 0), reason="GEOS < 3.12")
+def test_include_m_default():
+    # a few other corner cases for inferring dimensionality
+
+    # mixed XY and XYM -> XYM
+    _, coords, _ = shapely.to_ragged_array([line_string, line_string_m])
+    assert coords.shape[1] == 3
+
+    # mixed XY, XYM, and XYZM -> XYZM
+    _, coords, _ = shapely.to_ragged_array([line_string, line_string_m, line_string_zm])
+    assert coords.shape[1] == 4
+
+    # only empties -> always 2D
+    _, coords, _ = shapely.to_ragged_array([empty_line_string_m])
+    assert coords.shape[1] == 2
+    _, coords, _ = shapely.to_ragged_array([empty_line_string_zm])
+    assert coords.shape[1] == 2
+    # empty collection -> GEOS indicates 2D
+    _, coords, _ = shapely.to_ragged_array([empty_multi_polygon_m])
+    assert coords.shape[1] == 2
+    _, coords, _ = shapely.to_ragged_array([empty_multi_polygon_zm])
     assert coords.shape[1] == 2