WIP: ENH: add support for Measured geometry types

pyogrio/_compat.py

@@ -49,5 +49,7 @@
 
 HAS_GDAL_GEOS = __gdal_geos_version__ is not None
 
+GEOS_GE_312 = shapely is not None and shapely.geos_version >= (3, 12, 0)
+
 HAS_SHAPELY = shapely is not None and Version(shapely.__version__) >= Version("2.0.0")
 SHAPELY_GE_21 = shapely is not None and Version(shapely.__version__) >= Version("2.1.0")

pyogrio/_geometry.pyx

@@ -96,6 +96,7 @@ cdef str get_geometry_type(void *ogr_layer):
     if ogr_type not in GEOMETRY_TYPES:
         raise GeometryError(f"Geometry type is not supported: {ogr_type}")
 
+    """
     if OGR_GT_HasM(ogr_type):
         original_type = GEOMETRY_TYPES[ogr_type]
 
@@ -108,6 +109,8 @@ cdef str get_geometry_type(void *ogr_layer):
             f"Original type '{original_type}' "
             f"is converted to '{GEOMETRY_TYPES[ogr_type]}'")
 
+    """
+
     return GEOMETRY_TYPES[ogr_type]
 
 

pyogrio/_io.pyx

@@ -866,7 +866,9 @@ cdef validate_feature_range(
 
 @cython.boundscheck(False)  # Deactivate bounds checking
 @cython.wraparound(False)   # Deactivate negative indexing.
-cdef process_geometry(OGRFeatureH ogr_feature, int i, geom_view, uint8_t force_2d):
+cdef process_geometry(
+    OGRFeatureH ogr_feature, int i, geom_view, uint8_t force_2d, uint8_t keep_m
+):
 
     cdef OGRGeometryH ogr_geometry = NULL
     cdef OGRwkbGeometryType ogr_geometry_type
@@ -883,8 +885,9 @@ cdef process_geometry(OGRFeatureH ogr_feature, int i, geom_view, uint8_t force_2
             ogr_geometry_type = OGR_G_GetGeometryType(ogr_geometry)
 
             # if geometry has M values, these need to be removed first
-            if (OGR_G_IsMeasured(ogr_geometry)):
-                OGR_G_SetMeasured(ogr_geometry, 0)
+            if not keep_m:
+                if (OGR_G_IsMeasured(ogr_geometry)):
+                    OGR_G_SetMeasured(ogr_geometry, 0)
 
             if force_2d and OGR_G_Is3D(ogr_geometry):
                 OGR_G_Set3D(ogr_geometry, 0)
@@ -1074,7 +1077,8 @@ cdef get_features(
     int skip_features,
     int num_features,
     uint8_t return_fids,
-    bint datetime_as_string
+    bint datetime_as_string,
+    bint keep_m,
 ):
 
     cdef OGRFeatureH ogr_feature = NULL
@@ -1150,7 +1154,7 @@ cdef get_features(
                 fid_view[i] = OGR_F_GetFID(ogr_feature)
 
             if read_geometry:
-                process_geometry(ogr_feature, i, geom_view, force_2d)
+                process_geometry(ogr_feature, i, geom_view, force_2d, keep_m)
 
             process_fields(
                 ogr_feature, i, n_fields, field_data, field_data_view,
@@ -1185,7 +1189,8 @@ cdef get_features_by_fid(
     encoding,
     uint8_t read_geometry,
     uint8_t force_2d,
-    bint datetime_as_string
+    bint datetime_as_string,
+    bint keep_m,
 ):
 
     cdef OGRFeatureH ogr_feature = NULL
@@ -1235,7 +1240,7 @@ cdef get_features_by_fid(
                 raise FeatureError(str(exc))
 
             if read_geometry:
-                process_geometry(ogr_feature, i, geom_view, force_2d)
+                process_geometry(ogr_feature, i, geom_view, force_2d, keep_m)
 
             process_fields(
                 ogr_feature, i, n_fields, field_data, field_data_view,
@@ -1337,6 +1342,7 @@ def ogr_read(
     str sql_dialect=None,
     int return_fids=False,
     bint datetime_as_string=False,
+    bint keep_m=False,
 ):
 
     cdef int err = 0
@@ -1470,6 +1476,7 @@ def ogr_read(
                 read_geometry=read_geometry and geometry_type is not None,
                 force_2d=force_2d,
                 datetime_as_string=datetime_as_string,
+                keep_m=keep_m,
             )
 
             # bypass reading fids since these should match fids used for read
@@ -1503,7 +1510,8 @@ def ogr_read(
                 skip_features=skip_features,
                 num_features=num_features,
                 return_fids=return_fids,
-                datetime_as_string=datetime_as_string
+                datetime_as_string=datetime_as_string,
+                keep_m=keep_m,
             )
 
         ogr_types = [FIELD_TYPE_NAMES.get(field[1], "Unknown") for field in fields]
@@ -1602,6 +1610,7 @@ def ogr_open_arrow(
     int return_fids=False,
     int batch_size=0,
     use_pyarrow=False,
+    bint keep_m=False,
 ):
 
     cdef int err = 0

pyogrio/geopandas.py

@@ -7,12 +7,14 @@
 import numpy as np
 
 from pyogrio._compat import (
+    GEOS_GE_312,
     HAS_GEOPANDAS,
     PANDAS_GE_15,
     PANDAS_GE_20,
     PANDAS_GE_22,
     PANDAS_GE_30,
     PYARROW_GE_19,
+    SHAPELY_GE_21,
 )
 from pyogrio.errors import DataSourceError
 from pyogrio.raw import (
@@ -268,6 +270,8 @@ def read_dataframe(
 
     read_func = read_arrow if use_arrow else read
     gdal_force_2d = False if use_arrow else force_2d
+    keep_m = True if SHAPELY_GE_21 and GEOS_GE_312 else False
+    keep_m = True
     if not use_arrow:
         # For arrow, datetimes are read as is.
         # For numpy IO, datetimes are read as string values to preserve timezone info
@@ -289,6 +293,7 @@ def read_dataframe(
         sql=sql,
         sql_dialect=sql_dialect,
         return_fids=fid_as_index,
+        keep_m=keep_m,
         **kwargs,
     )
 

pyogrio/raw.py

@@ -54,6 +54,7 @@ def read(
     sql_dialect=None,
     return_fids=False,
     datetime_as_string=False,
+    keep_m=False,
     **kwargs,
 ):
     """Read OGR data source into numpy arrays.
@@ -153,6 +154,8 @@ def read(
         If True, will return datetime dtypes as detected by GDAL as a string
         array (which can be used to extract timezone info), instead of
         a datetime64 array.
+    keep_m : bool, optional (default: False)
+        If True, will keep the M dimension in geometries if present.
 
     **kwargs
         Additional driver-specific dataset open options passed to OGR.  Invalid
@@ -215,6 +218,7 @@ def read(
         return_fids=return_fids,
         dataset_kwargs=dataset_kwargs,
         datetime_as_string=datetime_as_string,
+        keep_m=keep_m,
     )
 
 
@@ -235,6 +239,7 @@ def read_arrow(
     sql=None,
     sql_dialect=None,
     return_fids=False,
+    keep_m=False,
     **kwargs,
 ):
     """Read OGR data source into a pyarrow Table.
@@ -309,6 +314,7 @@ def read_arrow(
         skip_features=gdal_skip_features,
         batch_size=batch_size,
         use_pyarrow=True,
+        keep_m=keep_m,
         **kwargs,
     ) as source:
         meta, reader = source
@@ -364,6 +370,7 @@ def open_arrow(
     return_fids=False,
     batch_size=65_536,
     use_pyarrow=False,
+    keep_m=False,
     **kwargs,
 ):
     """Open OGR data source as a stream of Arrow record batches.
@@ -456,6 +463,7 @@ def open_arrow(
         dataset_kwargs=dataset_kwargs,
         batch_size=batch_size,
         use_pyarrow=use_pyarrow,
+        keep_m=keep_m,
     )
 
 

pyogrio/tests/test_core.py

@@ -183,10 +183,9 @@ def test_list_layers(
 
     # Measured 3D is downgraded to plain 3D during read
     # Make sure this warning is raised
-    with pytest.warns(
-        UserWarning, match=r"Measured \(M\) geometry types are not supported"
-    ):
-        assert array_equal(list_layers(line_zm_file), [["line_zm", "LineString Z"]])
+    assert array_equal(
+        list_layers(line_zm_file), [["line_zm", "Measured 3D LineString"]]
+    )
 
     # Curve / surface types are downgraded to plain types
     assert array_equal(list_layers(curve_file), [["curve", "LineString"]])
@@ -597,6 +596,10 @@ def test_read_info_without_geometry(no_geometry_file):
     assert read_info(no_geometry_file)["total_bounds"] is None
 
 
+def test_read_info_zm(line_zm_file):
+    assert read_info(line_zm_file)["geometry_type"] == "Measured 3D LineString"
+
+
 @pytest.mark.parametrize(
     "name,value,expected",
     [

pyogrio/tests/test_geopandas_io.py

@@ -20,6 +20,7 @@
 from pyogrio._compat import (
     GDAL_GE_37,
     GDAL_GE_311,
+    GEOS_GE_312,
     HAS_ARROW_WRITE_API,
     HAS_PYPROJ,
     PANDAS_GE_15,
@@ -1106,6 +1107,18 @@ def test_read_sql_dialect_sqlite_gpkg(naturalearth_lowres, use_arrow):
     assert df.iloc[0].geometry.area > area_canada
 
 
+def test_read_zm(line_zm_file, use_arrow):
+    df = read_dataframe(line_zm_file, use_arrow=use_arrow)
+
+    if SHAPELY_GE_21 and GEOS_GE_312:
+        assert df.geometry.iloc[0].has_z
+        assert df.geometry.iloc[0].has_m
+        assert df.geometry.iloc[0].wkt.startswith("LINESTRING ZM (")
+    else:
+        assert df.geometry.iloc[0].has_z
+        assert df.geometry.iloc[0].wkt.startswith("LINESTRING Z (")
+
+
 @pytest.mark.parametrize(
     "encoding, arrow",
     [