Python OGR to Shapely – How to Convert Geometries with M-Values to Those Without

m-valuesogrosgeopythonshapely

Main question

Suppose I have a few OGR geometries, all of them containing M-values.

Is there a generic solution that allows me to translate them to Shapely geometries without M-values?

Shapely and M-values

Sadly, shapely and M-values don't mix (see this GitHub issue). So in order to generate the Shapely geometries, I'm willing to drop the M-values. However, I'm not sure what the best approach is.

Small example

from osgeo import ogr

pt_2D = ogr.Geometry(ogr.wkbPointM)
pt_2D.AddPointM(0,0,100)

line1_2D = ogr.Geometry(ogr.wkbLineStringM)
line1_2D.AddPointM(0, 0, 50)
line1_2D.AddPointM(1, 1, 250)

line2_2D = ogr.Geometry(ogr.wkbLineStringM)
line2_2D.AddPointM(3, 1, 150)
line2_2D.AddPointM(4, 1, 25)

mline_2D = ogr.Geometry(ogr.wkbMultiLineStringM)
mline_2D.AddGeometry(line1_2D)
mline_2D.AddGeometry(line2_2D)

pt_3D = ogr.Geometry(ogr.wkbPointM)
pt_3D.AddPointZM(0,0,0,100)

line1_3D = ogr.Geometry(ogr.wkbLineStringM)
line1_3D.AddPointZM(0, 0, 0, 50)
line1_3D.AddPointZM(1, 1, 1, 250)

line2_3D = ogr.Geometry(ogr.wkbLineStringM)
line2_3D.AddPointZM(3, 1, 3, 150)
line2_3D.AddPointZM(4, 1, 6, 25)

mline_3D = ogr.Geometry(ogr.wkbMultiLineStringM)
mline_3D.AddGeometry(line1_3D)
mline_3D.AddGeometry(line2_3D)

In the code chunk above, I've created a few 2D and 3D M-value-enabled OGR geometries. How can I translate them into Shapely geometries without the M-values?

First attempt: Exporting to WKT

My first attempt was to convert the geometries using the WKT strings and to use the shapely.wkt.loads function. However, when I use this approach on 2D geometries, the M-value gets incorrectly interpreted as a Z-value.
Interestingly, though, this approach does work for 3D geometries!

import shapely.wkt

def ogr_to_shapely(input_ogr_geom):
    return shapely.wkt.loads(input_ogr_geom.ExportToIsoWkt())

print(ogr_to_shapely(pt_2D).wkt)
# POINT Z (0 0 100) # wrong

print(ogr_to_shapely(line1_2D).wkt)
# LINESTRING Z (0 0 50, 1 1 250) # wrong

print(ogr_to_shapely(mline_2D).wkt)
# MULTILINESTRING Z ((0 0 50, 1 1 250), (3 1 150, 4 1 25)) # wrong

print(ogr_to_shapely(pt_3D).wkt)
# POINT Z (0 0 0) # right!

print(ogr_to_shapely(line1_3D).wkt)
# LINESTRING Z (0 0 0, 1 1 1) # right!

print(ogr_to_shapely(mline_3D).wkt)
# MULTILINESTRING Z ((0 0 0, 1 1 1), (3 1 3, 4 1 6)) # right!

Second attempt: Exporting coordinates

My second attempt was to just get the list of raw coordinates (without the M-values!) and pass them to Shapely. The problem, though, is finding which type of geometry to create after I've extracted the coordinate values.

def get_ogr_coords(input_ogr_geom):
    geom_type = input_ogr_geom.GetGeometryName()
    
    if geom_type[:len('multi')].lower() == 'multi':
        nested  = True
        coords = []
        geoms = [input_ogr_geom.GetGeometryRef(j) for 
                 j in range(input_ogr_geom.GetGeometryCount())]
        for geom in geoms:
            coords.append(geom.GetPoints())
    else:
        nested = False
        coords = list(input_ogr_geom.GetPoints())

    # Got stuck here. I can't find a generic solution for all geometry types.
    # I feel like I have to write one solution for every geometry type.

Third attempt: Exporting to WKT and using shapely.ops.transform

Using this post from miekwatt as inspiration, I tried to use the shapely.ops.transform function in tandem with the WKT-based transformation:

# Inspiration taken from https://gis.stackexchange.com/a/415879/123828
def drop_m(input_ogr_geom):
    shapely_geom = shapely.wkt.loads(input_ogr_geom.ExportToIsoWkt())
    def _drop_m(x, y, z):
        return x, y
    if input_ogr_geom.Is3D():
        return shapely_geom
    else:
        return shapely.ops.transform(_drop_m, shapely_geom)

print(drop_m(pt_2D).wkt)
# POINT (0 0) # right!

print(drop_m(line1_2D).wkt)
# LINESTRING (0 0, 1 1) # right!

print(drop_m(mline_2D).wkt)
# MULTILINESTRING ((0 0, 1 1), (3 1, 4 1)) # right!

print(drop_m(pt_3D).wkt)
# POINT Z (0 0 0) # right!

print(drop_m(line1_3D).wkt)
# LINESTRING Z (0 0 0, 1 1 1) # right!

print(drop_m(mline_3D).wkt)
# MULTILINESTRING Z ((0 0 0, 1 1 1), (3 1 3, 4 1 6)) # right!

This is the best approach so far! It works for all types of geometries, it's flexible enough to handle both 2D and 3D geometries, and it's pretty concise. However, it still relies on a WKT-based approach, which can be quite slow.

Back to the main question

Given all of the above, I'm brought back to the main question: is there a simple, quick and generic solution to transform an OGR geometry with M-values into a Shapely geometry without M-values other than the ones I tried in the code chunks above?

Best Answer

Here's a potential solution following suggestions from jbalk's and user30184's comments. It uses OGR's SetMeasured function to drop the M-values and relies on WKB (instead of WKT) to do the actual conversion between OGR and Shapely.

def from_ogr_to_shapely(ogr_geom):
    # Creating a copy of the input OGR geometry. This is done in order to 
    # ensure that when we drop the M-values, we are only doing so in a 
    # local copy of the geometry, not in the actual original geometry. 
    ogr_geom_copy = ogr.CreateGeometryFromWkb(ogr_geom.ExportToIsoWkb())
    
    # Dropping the M-values
    ogr_geom_copy.SetMeasured(False)

    # Generating a new shapely geometry
    # Note: I needed to typecast the results from
    # `ogr_geom_copy.ExportToIsoWkb()` to bytes since, in newer versions 
    # of GDAL/OGR and shapely, the outputs of that method were incompatible 
    # with the expected inputs for the `shapely.wkb.loads` method, 
    # leading to an error.
    shapely_geom = shapely.wkb.loads(bytes(ogr_geom_copy.ExportToIsoWkb()))

    return shapely_geom

As mentioned in the commented section in the code above, we need to create a local copy of the input OGR geometry. This is because the function used to drop the M-values (the SetMeasured function) actually modifies the original input geometry. So in order to ensure that our original geometry remains intact, we need to create a local copy. Then we can modify the local copy to our heart's content without fear of changes creeping upstream to the original OGR geometry.

Here's the function being applied to the inputs in the question:

print(from_ogr_to_shapely(pt_2D).wkt)
# POINT (0 0) # right!

print(from_ogr_to_shapely(line1_2D).wkt)
# LINESTRING (0 0, 1 1) # right!

print(from_ogr_to_shapely(mline_2D).wkt)
# MULTILINESTRING ((0 0, 1 1), (3 1, 4 1)) # right!

print(from_ogr_to_shapely(pt_3D).wkt)
# POINT Z (0 0 0) # right!

print(from_ogr_to_shapely(line1_3D).wkt)
# LINESTRING Z (0 0 0, 1 1 1) # right!

print(from_ogr_to_shapely(mline_3D).wkt)
# MULTILINESTRING Z ((0 0 0, 1 1 1), (3 1 3, 4 1 6)) # right!

Edit from the future

As pointed out by @ermaure in their answer, it seems like more recent versions of GDAL/OGR and shapely have changed things a bit and now the outputs of the ogr_geom.ExportToIsoWkb method are incompatible with the expected inputs for the shapely.wkb.loads method, causing the shapely.wkb.loads(ogr_geom_copy.ExportToIsoWkb()) statement to yield an error. This can be fixed by forcing the result of the ExportToIsoWkb method to be typecast to bytes. Explicitly, the command needs to change to shapely.wkb.loads(bytes(ogr_geom_copy.ExportToIsoWkb())). I have edited the code above to reflect that change.

Special thanks to @ermaure. Please be sure to upvote their answer as well!

Related Question