2

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?

asked Apr 25, 2023 at 0:02

3 Answers 3

1

Why not make it simpler (with a list of wkt strings)

my_gdf = gpd.read_file(filename=inGDBPath,layer=inLayerName)
wkt =[]
for this_inFeature in inLayer:
 this_inGeom = this_inFeature.GetGeometryRef()
 this_inWkt = this_inGeom.ExportToIsoWkt()
 wkt.append(this_inWkt)
my_gdf['WKT'] = wkt
my_gdf.WKT.head()
0 MULTILINESTRING M ((-95.999231733 32.62464663 ...
1 MULTILINESTRING M ((-97.934444075 32.710178812...
2 MULTILINESTRING M ((-100.39850042 32.450091716...
3 MULTILINESTRING M ((-96.840064431 32.642864073...
4 MULTILINESTRING M ((-99.677214454 32.460975142...
my_gdf.geometry.head()
0 MULTILINESTRING ((-95.99923 32.62465, -95.9974...
1 MULTILINESTRING ((-97.93444 32.71018, -97.9264...
2 MULTILINESTRING ((-100.39850 32.45009, -100.39...
3 MULTILINESTRING ((-96.84006 32.64286, -96.8400...
4 MULTILINESTRING ((-99.67721 32.46098, -99.6751...

The M values are defined for each point in the geometry (multiple values per line).

 inLayer.ResetReading()
 # first feature
 first = inLayer.GetNextFeature()
 mline_geom = first.GetGeometryRef()
 print(mline_geom.ExportToIsoWkt())
 MULTILINESTRING M ((-95.999231733 32.62464663 516.996,-95.99740274 32.6240395400001 517.111,-95.9889207699999 32.6211937 517.642,-95.988217821 32.620957828 517.686))
 for j in range(mline_geom.GetGeometryCount()):
 geom = mline_geom.GetGeometryRef(j)
 for i, point in enumerate(geom.GetPoints()):
 print(geom.GetM(i), point)
 516.9959999999992 (-95.99923173299999, 32.62464663000003)
 517.1110000000044 (-95.99740273999998, 32.62403954000007)
 517.6420000000071 (-95.98892076999994, 32.62119370000005)
 517.6860000000015 (-95.98821782099998, 32.62095782800003)

LineString.GetM() will return the value for the first point( point=0), by default.

geomi = mline_geom.GetGeometryRef(0)
geomi.GetM()
516.9959999999992
answered Apr 25, 2023 at 16:37
0

Building on gene's answer, I've gone ahead and made a couple of more complete functions that do everything that's needed by just using the OGR geometries.

First, I created a function called extract_geometry_info that takes in one OGR geometry and extracts a bunch of important/valuable information from it. That includes a set of tuples that contain the (X,Y) coordinates and the M-values.

Then, I modified the get_geoms function so that it applies the get_geoms() function to every OGR geometry in a given layer.

def extract_geometry_info(ogr_geom):
 '''
 Extracts all of the coordinates and M-values of a given OGR geometry
 Parameters
 ----------
 ogr_geom : osgeo.ogr.geometry
 OGR Geometry to be analyzed.
 Returns
 -------
 out_dict : dict 
 Dictionary containing multiple elements:
 geom_type : STR that describes the type of geometry
 num_sub_geoms : INT that describes how many sub-geometries are in 
 the full geometry (i.e., how many parts are in the 
 multi-part geometry). For LineStrings, this will 
 default to zero. 
 num_points : INT or TUPLE that describes the number of points in 
 the full geometry. For LineStrings, it is just an INT.
 For MultiLineStrings, it is a TUPLE of INTs.
 geom_coords : TUPLE with the geometry's (x,y) coordinates
 geom_mvals : TUPLE with the geomety's M-values
 geom_coords_and_mvals : TUPLE with the geometry's (x, y, M-value) 
 data
 upcast_geom_data : TUPLE similar to the one in the 
 geom_coords_and_mvals variable, but when 
 LineStrings are found, this geometry is "upcast"
 to a MultiLineString containing just one 
 sub-geometry
 '''
 geom_type = ogr_geom.GetGeometryName()
 geom_wkt = ogr_geom.ExportToIsoWkt()
 geom_wkb = ogr_geom.ExportToIsoWkb()
 
 if geom_type.lower()[:len('multi')] == 'multi':
 geom_coords = []
 geom_mvals = []
 geom_coords_and_mvals = []
 num_sub_geoms = ogr_geom.GetGeometryCount()
 num_points = []
 for j in range(ogr_geom.GetGeometryCount()):
 sub_geom = ogr_geom.GetGeometryRef(j)
 sub_geom_coords = tuple([pts for pts in sub_geom.GetPoints()])
 sub_geom_mvals = tuple([sub_geom.GetM(i) for i in 
 range(len(sub_geom_coords))])
 sub_geom_coords_and_mvals = tuple([(*sub_geom_coords[i],
 sub_geom_mvals[i]) for i in 
 range(len(sub_geom_coords))])
 geom_coords.append(sub_geom_coords)
 geom_mvals.append(sub_geom_mvals)
 geom_coords_and_mvals.append(sub_geom_coords_and_mvals)
 num_points.append(len(sub_geom_mvals))
 geom_coords = tuple(geom_coords)
 geom_mvals = tuple(geom_mvals)
 geom_coords_and_mvals = tuple(geom_coords_and_mvals)
 num_points = tuple(num_points)
 upcast_geom_data = geom_coords_and_mvals
 else:
 geom_coords = tuple([pts for pts in ogr_geom.GetPoints()])
 geom_mvals = tuple([ogr_geom.GetM(i) for i in range(len(geom_coords))])
 geom_coords_and_mvals = tuple([(*geom_coords[i],
 geom_mvals[i]) for i in 
 range(len(geom_coords))])
 num_sub_geoms = 0
 num_points = len(geom_mvals)
 upcast_geom_data = tuple([geom_coords_and_mvals])
 
 out_dict = {'geom_type':geom_type,
 'WKB':geom_wkb,
 'WKT':geom_wkt,
 'num_sub_geoms':num_sub_geoms,
 'num_points':num_points,
 'geom_coords':geom_coords,
 'geom_mvals':geom_mvals,
 'geom_coords_and_mvals':geom_coords_and_mvals,
 'upcast_geom_data':upcast_geom_data}
 return(out_dict)
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.
 geom_type : a column of STR values that describes the type of 
 geometry
 num_sub_geoms : a column of INT values that describes how many 
 sub-geometries are in the full geometry
 num_points : a column of INT or TUPLE values that describes the 
 number of points in the full geometries. For 
 LineStrings, each row contains INTs. For 
 MultiLineStrings, each row contains a TUPLE of INTs.
 geom_coords : a column of tuples with the geometries' (x,y) 
 coordinates
 geom_mvals : a column of tuples with the geometries' M-values
 geom_coords_and_mvals : a column of tuples with the geometries' 
 (x, y, M-value) data
 upcast_geom_data : column of tuples similar to the ones in the 
 geom_coords_and_mvals variable, but when 
 LineStrings are found, this geometry is "upcast"
 to a MultiLineString containing just one 
 sub-geometry.
 '''
 # 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 and geometry
 this_FID = this_inFeature.GetFID()
 this_geom = this_inFeature.GetGeometryRef()
 
 # Extracting the rest of the geometry's data
 this_geom_info = extract_geometry_info(this_geom)
 
 # Appending this feature to the list 
 features_list.append({
 inLayerIDColname:this_FID,
 'OGR_Geom':ogr.CreateGeometryFromWkb(this_geom_info['WKB']),
 'WKB':this_geom_info['WKB'],
 'WKT':this_geom_info['WKT'],
 'geom_type':this_geom_info['geom_type'],
 'num_sub_geoms':this_geom_info['num_sub_geoms'],
 'num_points':this_geom_info['num_points'],
 'geom_coords':this_geom_info['geom_coords'],
 'geom_mvals':this_geom_info['geom_mvals'],
 'geom_coords_and_mvals':this_geom_info['geom_coords_and_mvals'],
 'upcast_geom_data':this_geom_info['upcast_geom_data']})
 
 # Releasing the input Feature
 this_inFeature = None
 
 # Releasing the input file
 inDataSource = None
 inLayer = None
 
 # Generating a Pandas DataFrame and setting the Feature ID as the index
 features_df = pd.DataFrame(features_list).set_index(inLayerIDColname)
 
 return(features_df)

Example of the extract_geometry_info function

Here's an example where I apply the extract_geometry_info described above to OGR geometries containing a LineString and a MultiLineString:

# Creating the OGR Geometry
ogr_geom = ogr.CreateGeometryFromWkt('LINESTRING M (100 100 25, 120 105 50, 140 98 55, 160 103 100)')
# Extracting all of the info from the OGR geometry 
extract_geometry_info(ogr_geom)
print(extract_geometry_info)
# {'geom_type': 'LINESTRING',
# 'WKB': bytearray(b'\x01\xd2\x07\x00\x00\x04\x00\x00\x00\x00\x00\x00\x00\x00\x00Y@\x00\x00\x00\x00\x00\x00Y@\x00\x00\x00\x00\x00\x009@\x00\x00\x00\x00\x00\x00^@\x00\x00\x00\x00\x00@Z@\x00\x00\x00\x00\x00\x00I@\x00\x00\x00\x00\x00\x80a@\x00\x00\x00\x00\x00\x80X@\x00\x00\x00\x00\x00\x80K@\x00\x00\x00\x00\x00\x00d@\x00\x00\x00\x00\x00\xc0Y@\x00\x00\x00\x00\x00\x00Y@'),
# 'WKT': 'LINESTRING M (100 100 25,120 105 50,140 98 55,160 103 100)',
# 'num_sub_geoms': 0,
# 'num_points': 4,
# 'geom_coords': ((100.0, 100.0),
# (120.0, 105.0),
# (140.0, 98.0),
# (160.0, 103.0)),
# 'geom_mvals': (25.0, 50.0, 55.0, 100.0),
# 'geom_coords_and_mvals': ((100.0, 100.0, 25.0),
# (120.0, 105.0, 50.0),
# (140.0, 98.0, 55.0),
# (160.0, 103.0, 100.0)),
# 'upcast_geom_data': (((100.0, 100.0, 25.0),
# (120.0, 105.0, 50.0),
# (140.0, 98.0, 55.0),
# (160.0, 103.0, 100.0)),)}

# Creating the OGR Geometry
ogr_geom = ogr.CreateGeometryFromWkt('MULTILINESTRING M ((100 100 25, 120 105 50), (140 98 55, 160 103 100))')
# Extracting all of the info from the OGR geometry 
extract_geometry_info(ogr_geom)
# {'geom_type': 'MULTILINESTRING',
# 'WKB': bytearray(b'\x01\xd5\x07\x00\x00\x02\x00\x00\x00\x01\xd2\x07\x00\x00\x02\x00\x00\x00\x00\x00\x00\x00\x00\x00Y@\x00\x00\x00\x00\x00\x00Y@\x00\x00\x00\x00\x00\x009@\x00\x00\x00\x00\x00\x00^@\x00\x00\x00\x00\x00@Z@\x00\x00\x00\x00\x00\x00I@\x01\xd2\x07\x00\x00\x02\x00\x00\x00\x00\x00\x00\x00\x00\x80a@\x00\x00\x00\x00\x00\x80X@\x00\x00\x00\x00\x00\x80K@\x00\x00\x00\x00\x00\x00d@\x00\x00\x00\x00\x00\xc0Y@\x00\x00\x00\x00\x00\x00Y@'),
# 'WKT': 'MULTILINESTRING M ((100 100 25,120 105 50),(140 98 55,160 103 100))',
# 'num_sub_geoms': 2,
# 'num_points': (2, 2),
# 'geom_coords': (((100.0, 100.0), 
# (120.0, 105.0)),
# ((140.0, 98.0), 
# (160.0, 103.0))),
# 'geom_mvals': ((25.0, 50.0), 
# (55.0, 100.0)),
# 'geom_coords_and_mvals': (((100.0, 100.0, 25.0), 
# (120.0, 105.0, 50.0)),
# ((140.0, 98.0, 55.0), 
# (160.0, 103.0, 100.0))),
# 'upcast_geom_data': (((100.0, 100.0, 25.0), 
# (120.0, 105.0, 50.0)),
# ((140.0, 98.0, 55.0), 
# (160.0, 103.0, 100.0)))}

This is still not super fast, but it gets the job done.

Warning: no checks are being carried out!

It's worth pointing out that the functions I provide here are not running ANY kind of checks. They simply ASSUME that the input data will be exactly the way it should. When using this in real world cases, we might need to consider expanding these a bit to check all the inputs. For example, we need to ensure that the geometries are all LineString or MultiLineString.

answered May 1, 2023 at 21:58
0

I was following this thread for a while and modified Felipe's code to fit my use case with GeoPandas and Shapely.

WARNING

this code works with my work at the moment, but I have not conducted any significant testing, might be buggy. I did my best to validate the input geometries are only LINESTRING or MULTILINESTRING

import geopandas as gpd
import numpy as np
import shapely
from osgeo import ogr
def get_linestrings_m(gdb_path: str,
 layername: str,
 file_driver: str = 'OpenFileGdb' # using the open-source filedriver, can use alternatives
 )->gpd.GeoDataFrame:
 
 # error checking
 driver = ogr.GetDriverByName(file_driver)
 if driver is None: 
 raise ValueError("Driver not found!")
 
 data_source = driver.Open(gdb_path, 0) # 0 is read-mode, 1 write
 feature_layer = data_source.GetLayerByName(layername)
 if feature_layer is None:
 raise ValueError(f"Feature layer named {layername} not found!")
 
 layer_ID_columnname = feature_layer.GetFIDColumn()
 layer_spatial_ref = feature_layer.GetSpatialRef() # fetch the spatial reference object for the file
 crs = f"EPSG:{layer_spatial_ref.GetAuthorityCode(None)}" # create the EPSG code
 feature_layer.ResetReading()
 features_list = []
 for this_feature in feature_layer:
 this_feature_id = this_feature.GetFID()
 this_geometry = this_feature.GetGeometryRef() # should be either LINESTRING or MULTILINESTRING
 this_line_id = this_feature.GetField('column_name') # update if this name has changed
 this_geom_info = get_geometry_info(ogr_geom=this_geometry)
 features_list.append({layer_ID_columnname:this_feature_id,
 'LINE_ID': this_line_id,
 'geometry': this_geom_info['geom_shapely_repr'],
 'geometry M-values': this_geom_info['geom_m_values'],
 'num_points': this_geom_info['num_points']
 })
 this_feature = None # reset variable for next iteration
 data_source = None # releasing resources
 feature_layer = None # releasing resources
 features_df = pd.DataFrame(features_list).explode(['geometry', 'geometry M-values']).reset_index(drop=True) # explode the geometry and geometry M-values, these should be the same shape
 features_gdf = gpd.GeoDataFrame(data=features_df, geometry=features_df['geometry'], crs=crs)
 return features_gdf
# called for each feature in a layer
def get_geometry_info(ogr_geom: ogr.Geometry)->dict:
 geom_type = ogr_geom.GetGeometryName() # Geometry Attributes
 
 # error checking
 if 'linestring' not in geom_type.lower():
 raise ValueError(f"Error! Found {geom_type}\n Only linestring/ multi-linestring supported at this time!")
 if 'multi' in geom_type.lower(): # MULTILINESTRING
 geom_m_values = []
 num_points = []
 geom_shapely_repr = []
 num_sub_geometries = ogr_geom.GetGeometryCount()
 # extend() rather than append() so we only explode once
 for i in range(num_sub_geometries):
 linestring_info_i = process_linestring_info(ogr_geom.GetGeometryRef(i))
 geom_m_values.extend(linestring_info_i['geom_m_values'])
 num_points.extend(linestring_info_i['num_points'])
 geom_shapely_repr.extend(linestring_info_i['geom_shapely_repr'])
 out_dict = {'num_points':num_points, 
 'geom_shapely_repr':geom_shapely_repr, 
 'geom_m_values': geom_m_values 
 } 
 else: # LINESTRING
 out_dict = process_linestring_info(ogr_geom)
 
 return out_dict
# called for each linestring in a feature (linestring or multilinestring)
def process_linestring_info(linestring_geom: ogr.Geometry)->dict:
 geom_n_points = linestring_geom.GetPointCount()
 geom_coords = np.array(linestring_geom.GetPoints(nCoordDimension=2)) # selecting only x,y dimensions and casting to numpy array (https://gdal.org/java/org/gdal/ogr/Geometry.html#GetPoints(int))
 geom_shapely_points = shapely.points(geom_coords) # list of shapely points (https://shapely.readthedocs.io/en/stable/reference/shapely.points.html#shapely.points)
 geom_m_values = [linestring_geom.GetM(idx) for idx in range(geom_n_points)] # list of M-valuess
 out_dict = {'num_points':[geom_n_points], # casting to a list for compatability with extend
 'geom_shapely_repr':geom_shapely_points,
 'geom_m_values': geom_m_values
 }
 return out_dict
answered May 21, 2024 at 20:25

Your Answer

Draft saved
Draft discarded

Sign up or log in

Sign up using Google
Sign up using Email and Password

Post as a guest

Required, but never shown

Post as a guest

Required, but never shown

By clicking "Post Your Answer", you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.