GeoPandas/shapely/fiona ignore the M-values stored in GIS layers. How can we get around that?

Core issue

There are some files that have M-value data embedded inside them. However, GeoPandas/shapely/fiona ignore that information when reading it into memory (see here for a tiny compilation of info on the issue).

Real world example

A good example of this is the Texas Roadway Inventory. If we look at the 2022 vintage, specifically inside the TxDOT_Roadway_Linework layer stored in the 2022_Roadway_Inventory.gdb database, we'll see super long geometries, some of which span almost the entire state. The geometry data inside that layer have M-values stored for each point of every geometry, but given GeoPandas' upstream limitations, those M-values get ignored. This case in particular also doesn't have those beginning or ending mile marker columns inside them that are needed for linref to work. This means that even though the layer does indeed have all the ingredients needed, we just wouldn't be able to get that layer to work with linref.

It's also worth mentioning that besides cases like these where the beginning and ending columns are missing, there are other cases where linref's general assumption that the M-values scale linearly according to the X & Y coordinates doesn't hold true, so even when the beginning and ending columns do exist, they might lead us to missing a lot of granularity of how the M-values vary along each point inside the geometry.

Reproducible code

Here's a small example:

import geopandas as gpd
import pandas as pd
from osgeo import ogr, osr

# Setting up synthetic data
df = pd.DataFrame({'row_id':[1,2,3],
                   'route_id':['A','B','C'],
                   'wkt':["MULTILINESTRING M((0 0 10, 1 0 200, 2 0 300),(3 0 400, 4 0 500))",
                          "MULTILINESTRING M((10 10 100, 11 10 200, 12 10 300),(13 10 400, 14 10 500))",
                          "MULTILINESTRING M((20 20 100, 21 20 200, 22 20 300),(23 20 400, 24 20 5000))"
                       ],
                   'beg':[10,100,100],
                   'end':[500,500,5000]})
# Note that the M-values are super non-linear in a couple of locations. For example, 
# in the first geometry, the first point is at (0,0) with M-value 10, the second point is at (1,0) with 
# M-value 200 and the third point is at (3,0) with M-value 300. That means that the geometric 
# distance between the first two points is 1 and the M-value distance is 190, while the geometric
# distance between the second and third points is 1 and the M-value distance is 100. 

# Location where geodata will be written
out_gpkg = "c:/temp/temp_data.gpkg"
out_layer_name = "temp_data"

# Establishing the parameters needed to write to a new file using OGR
driver = ogr.GetDriverByName("GPKG")
data_source = driver.CreateDataSource(out_gpkg)
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)

layer = data_source.CreateLayer(out_layer_name, srs, ogr.wkbMultiLineStringM)

layer.CreateField(ogr.FieldDefn("row_id", ogr.OFTInteger))

field_data = ogr.FieldDefn("route_id", ogr.OFTString)
field_data.SetWidth(24)
layer.CreateField(field_data)

field_wkt = ogr.FieldDefn("wkt", ogr.OFTString)
field_wkt.SetWidth(1024)
layer.CreateField(field_wkt)

layer.CreateField(ogr.FieldDefn("beg", ogr.OFTReal))
layer.CreateField(ogr.FieldDefn("end", ogr.OFTReal))

for i,row in df.iterrows():
    # Create the feature
    feature = ogr.Feature(layer.GetLayerDefn())
    
    # Set the attributes using the values from the `df` object
    for col in df.columns: 
        feature.SetField(col, row[col])
    
    # Create the geometry from the wkt text
    ogr_geom = ogr.CreateGeometryFromWkt(row['wkt'])
    
    # Set the feature geometry using the newly-created `ogr_geom` object
    feature.SetGeometry(ogr_geom)
    
    # Add the feature to the layer
    layer.CreateFeature(feature)
    
    feature = None

data_source = None

# Reading in the newly-created GIS layer into a GeoDataFrame
gdf = gpd.read_file(out_gpkg, layer=out_layer_name)

# Printing WKTs
[print(geom.wkt) for geom in gdf.geometry]
# MULTILINESTRING ((0 0, 1 0, 2 0), (3 0, 4 0))
# MULTILINESTRING ((10 10, 11 10, 12 10), (13 10, 14 10))
# MULTILINESTRING ((20 20, 21 20, 22 20), (23 20, 24 20))

The process above creates a GDB file with a layer inside it that contains M-values that don't behave linearly.

For example. in the first geometry, the first point is at (0,0) with M-value 10, the second point is at (1,0) with M-value 200 and the third point is at (3,0) with M-value 300. That means that the geometric distance between the first two points is 1 and the M-value distance is 190, while the geometric distance between the second and third points is 1 and the M-value distance is 100. This nuance will be completely lost if we only rely on the beg and end columns in the data.

Also notice how, at the end, the M-values were simply dropped when they were read into the GeoDataFrame.

Addressing the issue (and my hesitance in doing so)

To address this, we might want to come up with an experimental method of reading in GIS layers using OGR. I've been toying around with this for a bit and have something we might be able to use, but it's super slow and probably prone to breaking.

Having said all of the above, I'd also like to mention that I don't think this is a feature that linref super needs. Frankly, I'd say that messing around with reading GIS layers feels WAY outside of linref's scope. So I'm fine if this stays eternally under an experimental branch. But I'm definitely going to take a crack at it, at the very least just for fun. It just bugs me so much that the data is right there, but GeoPandas/shapely/fiona simply ignores it.