[Python] Route Geometry for Road Detection is wrongly detecting nearby Streets in GPS analyzer...

I have developed a GPS analyzer pipeline where a csv with Lat,Lon,Speed, and Time is fed in the GPS analyzer pipeline. The pipeline generates a road compliance report of driver. The original GPS coordinates have noise upt0 50 meters. So, the backbone of the model is detection of road called route_geometry in my system. The problem is route_geometry of overly magnetic and it stick to nearby roads at the roundabouts or turns. How can I make this smooth. Please see the attached image. Orange line and red dots are original gps coordinates from tracker, the green line is the route_geometry.





Following is the code. perform_map_matching() function creates route_geometry which is a LineString:

def _build_optimized_graph(self, corridor_polygon, area_hash):
print(f"DEBUG: Starting optimized graph building with area_hash: {area_hash}")

import networkx as nx

G = nx.MultiDiGraph()
node_coords = {} # Will map coord_key -> OSM node ID

print(f" Processing {len(self.map_manager.map_data['roads'])} roads...")

# Pre-filter roads by bounding box (same as before)
bbox = corridor_polygon.bounds
candidate_roads = []

for road in self.map_manager.map_data['roads']:
geometry = road['geometry']
if len(geometry) < 2:
continue

# Quick bbox check first
road_minx = min(p['lon'] for p in geometry)
road_maxx = max(p['lon'] for p in geometry)
road_miny = min(p['lat'] for p in geometry)
road_maxy = max(p['lat'] for p in geometry)

if (road_maxx >= bbox[0] and road_minx <= bbox[2] and
road_maxy >= bbox[1] and road_miny <= bbox[3]):
candidate_roads.append(road)

print(f" {len(candidate_roads)} roads in bounding box, testing intersection...")

# Now test actual intersection only on candidates
relevant_roads = []
for road in candidate_roads:
road_line = LineString([(p['lon'], p['lat']) for p in road['geometry']])
if corridor_polygon.intersects(road_line):
relevant_roads.append(road)

print(f" ✅ {len(relevant_roads)} roads intersect corridor")

# Build graph from relevant roads using ORIGINAL OSM NODE IDs
for road in relevant_roads:
geometry = road['geometry']
node_refs = road.get('node_refs', [])

# Skip if we don't have node references
if len(node_refs) != len(geometry):
print(f" ⚠️ Skipping road {road.get('id')}: geometry/node_refs mismatch")
continue

road_nodes = []

for i, point in enumerate(geometry):
coord_key = (point['lat'], point['lon'])
osm_node_id = node_refs # Use ORIGINAL OSM node ID

if coord_key not in node_coords:
node_coords[coord_key] = osm_node_id
G.add_node(osm_node_id, y=point['lat'], x=point['lon'])

road_nodes.append(node_coords[coord_key])

# Create edges
for i in range(len(road_nodes) - 1):
u, v = road_nodes, road_nodes[i + 1]
G.add_edge(u, v, 0, # key=0
osmid=road['id'],
highway=road['tags'].get('highway', 'unclassified'),
maxspeed=road.get('speed_limit', 50))

G.graph['crs'] = 'EPSG:4326'
self.road_graph = G

# Cache the graph
cache_file = os.path.join(self.cache_dir, f"road_graph_{area_hash}.pkl")
os.makedirs(self.cache_dir, exist_ok=True)
try:
with open(cache_file, 'wb') as f:
pickle.dump({'graph': G}, f)
print(f" Cached graph to {cache_file}")
except Exception:
pass

return G

def _perform_map_matching(self, df):

# CLEAR previous week's edge data to avoid KeyError with new graph
"""
Optimized map matching that combines PBF data with old version's efficiency
"""
try:
# 1) Build corridor (same as old version)
coords = list(zip(df['lon'].tolist(), df['lat'].tolist()))
raw_line = LineString(coords)

# Use simple buffer in lat/lon (faster than projection)
buffer_degrees = 0.0003 # ~33 meters at France latitude
corridor_ll = raw_line.buffer(buffer_degrees)

print(" Building optimized road network from PBF data...")

# 2) Get or build cached graph for this area
bbox = corridor_ll.bounds
area_hash = f"{bbox[0]:.4f}_{bbox[1]:.4f}_{bbox[2]:.4f}_{bbox[3]:.4f}"

if not self._load_cached_graph(area_hash):
self.road_graph = self._build_optimized_graph(corridor_ll, area_hash)

# 3) Use OSMnx's optimized nearest_edges (same as old version)
print(" Matching GPS points to roads...")
X = df['lon'].to_numpy()
Y = df['lat'].to_numpy()

# This is the same efficient call as old version
edge_keys, dists = ox.distance.nearest_edges(
self.road_graph,
X,
Y,
return_dist=True
)

# 4) Rest of processing (identical to old version)
matched_edges = []
route_points = []
nodes = self.road_graph.nodes
get_node_xy = lambda n: (nodes[n]['x'], nodes[n]['y'])

current_road_id = None
previous_edge = None

for i, ek in enumerate(edge_keys):
if dists is not None and dists <= 50:
# Use simplified edge selection (avoid complex helper methods)
selected_edge = self._select_best_edge(
ek, current_road_id, previous_edge
)

matched_edges.append(selected_edge)
if selected_edge:
self.driven_edges.add(selected_edge)
current_road_id = self._get_road_id(selected_edge)
previous_edge = selected_edge

# Geometry handling (same as old version)
geom = self._get_edge_geometry_fast(selected_edge)

if geom is not None:
pt = Point(X, Y)
try:
projected = geom.interpolate(geom.project(pt))
route_points.append((projected.y, projected.x))
except Exception:
route_points.append((Y, X))
else:
route_points.append((Y, X))
else:
matched_edges.append(None)
route_points.append((Y, X))
current_road_id = None
previous_edge = None

# 5) Write results (same as old version)
df['matched_edge'] = matched_edges
df['road_matched_lat'] = [p[0] for p in route_points]
df['road_matched_lon'] = [p[1] for p in route_points]

# 6) Route geometry (same as old version)
valid_points = []
for lat, lon in route_points:
if lat and lon and not (lat == 0 or lon == 0):
valid_points.append((lon, lat))
print(f" DEBUG: route_points sample: {route_points[:5]}")
print(f" DEBUG: valid_points count: {len(valid_points)}")
if len(valid_points) <= 1:
print(f" DEBUG: All route_points: {route_points}")
print(f" DEBUG: Y array sample: {Y[:5]}")
print(f" DEBUG: X array sample: {X[:5]}")

if len(valid_points) <= 1:
print(" ⚠️ No valid route points found, using original GPS coordinates")
fallback_points = []
for i in range(len(df)):
if pd.notna(df.iloc['lat']) and pd.notna(df.iloc['lon']):
fallback_points.append((df.iloc['lon'], df.iloc['lat']))

if len(fallback_points) > 1:
self.route_geometry = LineString(fallback_points)
print(f" ✅ Created fallback route geometry with {len(fallback_points)} GPS points")

if len(valid_points) > 1:
try:
self.route_geometry = LineString(valid_points)
print(f" ✅ Created route geometry with {len(valid_points)} points")
self.route_geometry = self.route_geometry.simplify(0.000001, preserve_topology=True)

# Update road context
if hasattr(self, 'map_manager') and self.map_manager:
driven_road_ids = {self.road_graph.edges[edge].get('osmid')
for edge in self.driven_edges if edge}
driven_road_ids.discard(None)
if driven_road_ids:
self.map_manager.set_road_context(driven_road_ids, self.route_geometry)

except Exception as e:
print(f" ⚠️ Could not create route geometry: {type(e).__name__}: {str(e)}")
import traceback
print(f" Full error: {traceback.format_exc()}")

self.route_geometry = None

except Exception as e:
print(f" ⚠️ Map matching failed: {e}. Using GPS coordinates.")
df['matched_edge'] = [None] * len(df)
df['road_matched_lat'] = df['lat']
df['road_matched_lon'] = df['lon']

def _load_cached_graph(self, area_hash):
"""Load cached road graph if available"""
cache_file = os.path.join(self.cache_dir, f"road_graph_{area_hash}.pkl")

if False:
print(" Loading cached road graph...")
try:
with open(cache_file, 'rb') as f:
graph_data = pickle.load(f)
self.road_graph = graph_data['graph']
# Restore any cached indices
if 'geometry_cache' in graph_data:
self.geometry_cache = graph_data['geometry_cache']
return True
except Exception as e:
print(f" ⚠️ Cache load failed: {e}")
return False
return False

def _select_best_edge(self, nearest_edge, current_road_id, previous_edge):
"""Simplified edge selection without expensive lookups"""
if not current_road_id:
return nearest_edge

# Check if nearest edge has same road ID
nearest_road_id = self._get_road_id(nearest_edge)
if nearest_road_id == current_road_id:
return nearest_edge

# For performance, just use nearest edge instead of expensive search
return nearest_edge

def _get_edge_geometry_fast(self, edge_key):
"""Fast geometry lookup with caching"""
if not edge_key:
return None

try:
edge_data = self.road_graph.edges[edge_key]
if 'geometry' in edge_data and edge_data['geometry'] is not None:
return edge_data['geometry']
else:
# Create simple line between nodes
u, v, k = edge_key
nodes = self.road_graph.nodes
ux, uy = nodes['x'], nodes['y']
vx, vy = nodes[v]['x'], nodes[v]['y']
return LineString([(ux, uy), (vx, vy)])
except Exception:
return None

def _is_likely_roundabout(self, edge_key, current_idx, X, Y, look_ahead=4):
"""
Detect if we're likely in a roundabout based on curvature and road type
"""
try:
# Check edge attributes for roundabout indicators
edge_data = self.road_graph.edges[edge_key]

# Check for roundabout tags
if 'junction' in edge_data and edge_data['junction'] == 'roundabout':
return True

# Check for circular geometry pattern
if current_idx >= 2 and current_idx < len(X) - look_ahead:
# Calculate curvature over several points
curvature = self._calculate_curvature(current_idx, X, Y, window=look_ahead)
# High curvature + short segments often indicate roundabouts
return curvature > 0.3 # Adjust threshold as needed

except Exception:
pass

return False

def _select_roundabout_edge(self, nearest_edge, current_idx, X, Y, previous_edge):
"""
Special logic for selecting edges in roundabouts
"""
if not previous_edge:
return nearest_edge

try:
# Get direction of travel
if current_idx >= 2:
travel_bearing = self._calculate_bearing(
X[current_idx-2], Y[current_idx-2],
X[current_idx], Y[current_idx]
)

# Find edges that align with travel direction
candidates = self._get_edges_near_point(X[current_idx], Y[current_idx], max_dist=30)

best_edge = nearest_edge
best_alignment = -1

for candidate_edge, dist in candidates[:5]: # Check top 5 closest
edge_bearing = self._get_edge_bearing(candidate_edge)
if edge_bearing is not None:
alignment = self._bearing_similarity(travel_bearing, edge_bearing)
if alignment > best_alignment:
best_alignment = alignment
best_edge = candidate_edge

return best_edge

except Exception:
pass

return nearest_edge

def _calculate_curvature(self, current_idx, X, Y, window=4):
"""
Calculate path curvature around current point
"""
try:
if current_idx < window or current_idx >= len(X) - window:
return 0

# Get points before and after
x1, y1 = X[current_idx - window], Y[current_idx - window]
x2, y2 = X[current_idx], Y[current_idx]
x3, y3 = X[current_idx + window], Y[current_idx + window]

# Calculate angle change
bearing1 = self._calculate_bearing(x1, y1, x2, y2)
bearing2 = self._calculate_bearing(x2, y2, x3, y3)

angle_change = abs(bearing2 - bearing1)
if angle_change > 180:
angle_change = 360 - angle_change

# Normalize by distance to get curvature
dist = ((x3-x1)**2 + (y3-y1)**2)**0.5
return angle_change / (dist * 111000) if dist > 0 else 0 # Convert to curvature per meter

except Exception:
return 0

def _calculate_bearing(self, x1, y1, x2, y2):
"""
Calculate bearing between two points in degrees
"""
dx = x2 - x1
dy = y2 - y1
bearing = math.degrees(math.atan2(dy, dx))
return (bearing + 360) % 360

def _get_edge_bearing(self, edge_key):
"""Fast geometry lookup with caching"""
if not edge_key:
return None

try:
edge_data = self.road_graph.edges[edge_key]
if 'geometry' in edge_data and edge_data['geometry'] is not None:
return edge_data['geometry']
else:
# Create simple line between nodes
u, v, k = edge_key
nodes = self.road_graph.nodes
ux, uy = nodes['x'], nodes['y']
vx, vy = nodes[v]['x'], nodes[v]['y']
return LineString([(ux, uy), (vx, vy)])
except Exception:
return None

def _bearing_similarity(self, bearing1, bearing2):
"""
Calculate similarity between two bearings (0-1, where 1 is identical)
"""
diff = abs(bearing1 - bearing2)
if diff > 180:
diff = 360 - diff
return 1 - (diff / 180)

def _get_distance_to_edge(self, lon, lat, edge_key):
"""
Calculate distance from point to edge in meters
"""
try:
edge_data = self.road_graph.edges[edge_key]
point = Point(lon, lat)

if 'geometry' in edge_data and edge_data['geometry']:
geom = edge_data['geometry']
else:
u, v, k = edge_key
u_data = self.road_graph.nodes[u]
v_data = self.road_graph.nodes[v]
geom = LineString([(u_data['x'], u_data['y']), (v_data['x'], v_data['y'])])

dist_deg = point.distance(geom)
return dist_deg * 111000 # Convert to meters

except Exception:
return float('inf')

# Keep your existing helper methods
def _get_road_id(self, edge_key):
"""Extract road identifier (name or osmid) from edge"""
if not edge_key:
return None
try:
edge_data = self.road_graph.edges[edge_key]
# Try road name first, then osmid
if 'name' in edge_data and edge_data['name']:
return edge_data['name']
elif 'osmid' in edge_data and edge_data['osmid']:
return edge_data['osmid']
else:
return f"edge_{edge_key[0]}_{edge_key[1]}" # fallback ID
except:
return None

def _find_same_road_nearby(self, lon, lat, target_road_id, max_dist=40):
"""Find edge with same road ID near the point"""
try:
candidates = self._get_edges_near_point(lon, lat, max_dist)

for edge_key, dist in candidates:
road_id = self._get_road_id(edge_key)
if road_id == target_road_id:
return edge_key
except Exception:
pass
return None

def _get_edges_near_point(self, lon, lat, max_dist=50):
"""Get all edges within distance of point"""
candidates = []
try:
point = Point(lon, lat)

# Get nearby nodes first
nearby_nodes = []
for node_id, node_data in self.road_graph.nodes(data=True):
node_point = Point(node_data['x'], node_data['y'])
dist_deg = point.distance(node_point)
dist_m = dist_deg * 111000 # rough conversion
if dist_m <= max_dist * 2: # wider search for nodes
nearby_nodes.append(node_id)

# Get edges from nearby nodes
for node in nearby_nodes:
for neighbor in self.road_graph.neighbors(node):
for key in self.road_graph[node][neighbor].keys():
edge_key = (node, neighbor, key)

# Calculate distance to edge
try:
edge_data = self.road_graph.edges[edge_key]
if 'geometry' in edge_data and edge_data['geometry']:
geom = edge_data['geometry']
else:
u_coord = (self.road_graph.nodes[node]['x'], self.road_graph.nodes[node]['y'])
v_coord = (self.road_graph.nodes[neighbor]['x'], self.road_graph.nodes[neighbor]['y'])
geom = LineString([u_coord, v_coord])

dist_deg = point.distance(geom)
dist_m = dist_deg * 111000

if dist_m <= max_dist:
candidates.append((edge_key, dist_m))
except Exception:
continue

except Exception:
pass

[url="https://stackoverflow.com/questions/79767140/route-geometry-for-road-detection-is-wrongly-detecting-nearby-streets-in-gps-ana"]Continue reading...[/url][/u]