Spatial Projection

This notebook covers projecting point and linear geometries onto a linearly referenced network.

import linref as lr

# Load roadways dataset with LRS pre-configured
roadways = lr.datasets.load('roadways', set_lrs=True)

Projecting Point Geometries

When you have point data with only spatial coordinates (no route or milepost information), you can project those points onto your linearly referenced network. This is useful for integrating external data sources like GPS points, sensor locations, or crash data that hasn’t been linearly referenced yet.

First, we’ll create some sample point geometries to simulate data without LRS information:

import geopandas as gpd
from shapely.geometry import Point

# First, ensure roadways have M-enabled geometries for projection
roadways_spatial = roadways.lr.add_geom_m()

# Create sample point features (e.g., sign locations) with only X/Y coordinates
# These points are near the roadway network but don't have route or milepost info
points = gpd.GeoDataFrame({
    'sign_id': [1, 2, 3],
    'sign_type': ['Speed Limit', 'Stop', 'Yield'],
    'geometry': [
        Point(1.5, 0.4),   # Near US-101
        Point(4.0, 6.1),   # Near SR-1
        Point(6.0, 11.8)   # Near I-5
    ]
}, crs=roadways.crs)

# Project points onto the roadway network
# buffer=2 searches for roads within 2 coordinate units of each point
projected = roadways_spatial.lr.project(points, buffer=2, nearest=True, dropna=True)

print("Points after projection (now have route and loc columns):")
print(projected[['sign_id', 'sign_type', 'route', 'loc', 'project_distance']])

Points after projection (now have route and loc columns):
   sign_id    sign_type   route       loc  project_distance
0        1  Speed Limit  US-101  1.519231          0.098058
1        2         Stop    SR-1  4.026450          0.127467
2        3        Yield     I-5  6.038462          0.114354

Projecting Linear Geometries

When you need to match linear geometries from one dataset to another (e.g., conflating road data from different sources), you can use parallel_project_hausdorff(). This function uses Hausdorff distance to find the best geometric match between linear features.

This is useful for:

• Conflating road inventories from different agencies

• Matching GPS traces to a road network

• Integrating pavement surveys onto a master road system

Here’s an example of creating alternative road segments and projecting them onto our master roadway network:

from linref.ext.spatial import parallel_project_hausdorff
from shapely.geometry import LineString

# Create sample linear segments from another source (e.g., county road inventory)
# These have geometries but no route IDs or mileposts from our system
other_segments = gpd.GeoDataFrame({
    'segment_id': ['A', 'B'],
    'data_source': ['County', 'County'],
    'geometry': [
        LineString([(1.0, 0.2), (3.0, 0.6)]),   # Should match US-101
        LineString([(4.0, 5.9), (6.0, 6.4)])     # Should match SR-1
    ]
}, crs=roadways.crs)

# Project these segments onto our master roadway network
# The function finds the best geometric match using Hausdorff distance
projected_segments = parallel_project_hausdorff(
    target=roadways_spatial,       # Our master roadway network with M-geometries
    projected=other_segments,      # Segments to be matched
    buffer=2,                      # Search within 2 coordinate units of segment endpoints
    max_distance=1,                # Only accept matches with similarity < 1 unit
    match=1,                       # Return best match only
    densify=0.1                    # Add vertices for accurate matching
)

print("Segments after projection (now have route, beg, and end from our LRS):")
print(projected_segments[['segment_id', 'data_source', 'route', 'beg', 'end']])

Segments after projection (now have route, beg, and end from our LRS):
  segment_id data_source   route       beg      end
0          A      County  US-101  1.000000  2.50000
1          B      County    SR-1  3.985852  6.00123