Problem statement
I am aware that geopandas (due to the limitations of shapely and pygeos) does not read geometries with M-values - when reading a layer that has them, geopandas now just drops those values.
(For historical reference: Until recently, when reading 2D features with M-values, geopandas would mistake the M-value for Z coordinate values - it was a known shortcoming related to a bug way upstream in the dependency tree: GEOS->pygeos->shapely->geopandas. That bug has been recently fixed in GEOS and the update has already been pushed downstream to pygeos. The folks at shapely are still working on updates to make M-values work with shapely, as seen here and here).
That being the case, how can I quickly extract the geometries from a feature layer in a way that retains the M-values? My goal is to be able to easily access the geometry's coordinates and M-values. I'm working on linear referencing-related stuff, in case that matters.
To put some constraints on the problem: suppose that I'm working exclusively with Feature Layers inside a GeoDataBase (.gdb) and that the layers contain exclusively either LineString or MultiLineString geometries.
What I've done so far
I've developed a function/method called get_geoms (described below) that uses the osgeo/ogr library to read in a layer's rows and store their WKT, WKB and OGR geometries in a pandas DataFrame:
from osgeo import ogr
import pandas as pd
def get_geoms(inGDBPath,
inLayerName):
'''
Parameters
----------
inGDBPath : str or pathlib.Path
Path to the GDB file where the feature layer is stored.
inLayerName : str
name of the layer where the features are stored.
Returns
-------
features_df : pandas.DataFrame
DataFrame that contains four columns:
FID : a "Field ID" column. The name of this column actually depends
on the column name used in the original layer. The column
will typically contain integers from 1 to N.
WKT : a column of STR values containing the well-known-text (WKT)
version of the geometry.
WKB : a column of BYTEARRAY values containing the well-known-binary
(WKB) version of the geometry.
OGR_Geom : a column of osgeo.ogr.Geometry values containing the
OGR geometries of each feature in the layer. This also
includes the M-values in cases where they are present.
'''
# use OGR specific exceptions
ogr.UseExceptions()
# Definitions for input file name and layer name
inDriverName = "OpenFileGDB"
inDriver = ogr.GetDriverByName(inDriverName)
inDataSource = inDriver.Open(inGDBPath, 0)
inLayer = inDataSource.GetLayerByName(inLayerName)
inLayerIDColname = inLayer.GetFIDColumn()
# Making sure the feature reader is reset
inLayer.ResetReading()
# Generating an empty list to hold all the features
features_list = []
# Iterating over every feature in the input layer
for this_inFeature in inLayer:
# Extracting the input feature's ID, geometry, WKT and WKB
this_FID = this_inFeature.GetFID()
this_inGeom = this_inFeature.GetGeometryRef()
this_inWkt = this_inGeom.ExportToIsoWkt()
this_inWkb = this_inGeom.ExportToIsoWkb()
# Copying the geometry
this_inGeom_copy = ogr.CreateGeometryFromWkb(this_inWkb)
# Appending this feature to the list
features_list.append({inLayerIDColname:this_FID,
'WKT':this_inWkt,
'WKB':this_inWkb,
'OGR_Geom':this_inGeom_copy})
# Releasing the input Feature
this_inFeature = None
# Releasing the input file
inDataSource = None
inLayer = None
features_df = pd.DataFrame(features_list).set_index(inLayerIDColname)
return(features_df)
Applying my method to test data
Here is the get_geoms function being applied to this data (which is a subset of this significantly larger dataset):
import geopandas as gpd
inGDBPath = 'C:/Temp/mytemp.gdb'
inLayerName = 'I20_Temp'
# Loading the GIS layer into a GeoPandas GeoDataFrame
my_gdf = gpd.read_file(filename=inGDBPath,
layer=inLayerName)
# Loading the extra data (WKT, WKB and OGR geometries)
extra_data = get_geoms(inGDBPath=inGDBPath,
inLayerName=inLayerName)
# Merging them into one big GeoDataFrame
full_gdf = my_gdf.merge(extra_data,
how='left',
left_index=True,
right_index=True)
The problems with my approach
First of all, the approach I show above is slow. For this small layer which has 1,354 rows and 202 columns, it takes about 114 ms (timed with the %timeit magic command). But when applied to the full original dataset which contains 883,837 rows and 202 columns, it takes almost 100 seconds. This is because my implementation relies on slow Python-based loops to read in the data from the GIS layer.
The second problem is that there is no great way to manipulate the data that I've extracted. WKTs and WKBs are only really useful to convert geometries (for example, from OGR geometry to shapely geometry), but they aren't very useful on their own. And it's a pain to use the OGR geometries because they're structured in such an odd way (at least to me). For example, navigating the coordinates of multi-part features is annoying as heck. I even created a function that extracts the (X,Y) coordinates and M-values from the OGR geometries, but my implementation is super clunky.
Back to the main question
So, is there a well-established Python library I could use to quickly extract and easily manipulate geometries from a GIS layer in a way that preserves its M-values?
