[Request] Add support for "Measured 3D LineString"

[DeflateAwning]

When I import a dataset with a "Measured 3D LineString" geometry, I receive a UserWarning: Measured (M) geometry types are not supported. Original type 'Measured 3D LineString' is converted to 'LineString Z'.

I would like to request that support for Measured 3D LineStrings be added, as they are useful for a variety of applications.

The userwarning is generated in pyogrio/raw.py:194

An example dataset (albeit a rather large one: 70 MB and 112,995 rows) is here: https://static.aer.ca/prd/documents/sts/st37/Township_1-25.zip (specifically, the ST37_WG_GCS_NAD83_TWP_01_25.shp shapefile).

[brendan-ward]

As I recall, up until GEOS 3.12, GEOS did not support ZM coordinates, so we had no way to do anything with them and opted to force them to drop their M value when reading using GDAL. In order for us to use in a meaningful way, ZM support needs to be available within Shapely, and this functionality would be limited to GEOS >= 3.12. We could probably add a parameter to allow return of ZM WKB instead of forcing the M value to be dropped, but that doesn't seem as useful on its own.

Related: shapely/shapely#1648, shapely/shapely#1808

[DeflateAwning]

The use case for me is just doing basic pandas transformations using geopandas, which depends on this lib, and then exporting them as WKT/WKB afterwards.

[brendan-ward]

In order to use in GeoPandas, we'd need to be able to handle the geometries in Shapely first. However, if you are only doing operations in Pandas (nothing involving the geometry), then it would be possible (once we've enabled it at a lower level) to use the pyogrio.raw.read interface and use the results to construct a Pandas DataFrame that has a WKB column.

[theroggy]

@DeflateAwning

In the mean time basic support for Measured geometries has been added in shapely 2.1.

If you have shapely >= 2.1 installed, and you use pyogrio with use_arrow=True, you can read the Measured information already.

The warning is still (incorrectly) given, but the data is returned with Z and M values.

E.g.:

import geopandas as gpd

path = "/vsizip/c:/temp/polygonm/Township_1-25.zip/ST37_WG_TWP_001_025/ST37_WG_GCS_NAD83_TWP_01_25.shp"
gdf = gpd.read_file(path, use_arrow=True)

print(gdf.head())

Output:

  Well_LicNo            UWI  ...    KBE                                           geometry
0    0000025  0104151907000  ...  791.0  LINESTRING ZM (-112.02404 49.8359 791 0, -112....     
1    0000055  0174131603020  ...  750.4  LINESTRING ZM (-111.73793 50.42805 750.4 0, -1...     
2    0000059  0154091505000  ...  763.2  LINESTRING ZM (-111.17012 50.25728 763.2 0, -1...     
3    0000060  0154091610000  ...  759.9  LINESTRING ZM (-111.18168 50.2609 759.9 0, -11...     
4    0000061  0154091006000  ...  760.8  LINESTRING ZM (-111.16419 50.24265 760.8 0, -1... 

[DeflateAwning]

What? Do you think that's good or right?

I'm not convinced

[theroggy]

What? Do you think that's good or right?

I'm not convinced

Not sure what you mean. Because the m values are all 0?

I don't know the data, but this passes correctly: the M value = 4 is written correctly and read correctly again:

wkt = "POINT ZM (1 2 3 4)"
gdf = gp.GeoDataFrame(geometry=shapely.from_wkt([wkt]), crs="EPSG:4326")
filename = tmp_path / "test.gpkg"
write_dataframe(gdf, filename, use_arrow=True)
result_gdf = read_dataframe(filename, use_arrow=True)
assert_geodataframe_equal(gdf, result_gdf)