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.
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
# 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?