#!/usr/bin/env python3 """ GeoJSON Road Comparison Script - Optimized Version Compares two GeoJSON files containing road data and identifies: 1. Roads in file1 that don't have corresponding coverage in file2 (removed roads) 2. Roads in file2 that don't have corresponding coverage in file1 (added roads) Only reports differences that are significant (above minimum length threshold). Optimized for performance with parallel processing and spatial indexing. TODO: - ignore points outside of lines - put properties properly on removed roads, so they're visible in JOSM - handle polygons properly (on previous geojson step?) for circular roads - handle C 44a -> County Road 44A - handle Trce -> Trace/Terrace? """ import json import argparse from pathlib import Path from typing import List, Dict, Any, Tuple import geopandas as gpd from shapely.geometry import LineString, MultiLineString, Point, Polygon from shapely.ops import unary_union from shapely.strtree import STRtree import pandas as pd import warnings import multiprocessing as mp from functools import partial import numpy as np import gc import importlib qgisfunctions = importlib.import_module("qgis-functions") # Suppress warnings for cleaner output warnings.filterwarnings('ignore') import re def titlecase(s): return re.sub( r"[A-Za-z0-9]+('[A-Za-z0-9]+)?", lambda word: word.group(0).capitalize(), s) class RoadComparator: def __init__(self, tolerance_feet: float = 50.0, min_gap_length_feet: float = 100.0, n_jobs: int = None, chunk_size: int = 1000, exclude_unnamed: bool = False): """ Initialize the road comparator. Args: tolerance_feet: Distance tolerance for considering roads as overlapping (default: 50 feet) min_gap_length_feet: Minimum length of gap/extra to be considered significant (default: 100 feet) n_jobs: Number of parallel processes to use (default: CPU count - 1) chunk_size: Number of geometries to process per chunk (default: 1000) exclude_unnamed: Exclude features without name/highway tags from coverage (default: False) """ self.tolerance_feet = tolerance_feet self.min_gap_length_feet = min_gap_length_feet # Reduce worker count for Windows to prevent memory issues self.n_jobs = n_jobs or min(2, max(1, mp.cpu_count() // 2)) self.chunk_size = chunk_size self.exclude_unnamed = exclude_unnamed # Convert feet to degrees (approximate conversion for continental US) # 1 degree latitude ≈ 364,000 feet # 1 degree longitude ≈ 288,000 feet (at 40° latitude) self.tolerance_deg = tolerance_feet / 364000.0 self.min_gap_length_deg = min_gap_length_feet / 364000.0 if self.exclude_unnamed: print("Excluding unnamed features from coverage calculation") def _has_name(self, row) -> bool: """Check if a feature has a name tag (for OSM data filtering).""" # Check for OSM-style tags (stored as JSON string) if 'tags' in row.index: tags = row.get('tags') if isinstance(tags, dict): return bool(tags.get('name')) elif isinstance(tags, str): # Tags stored as JSON string try: tags_dict = json.loads(tags) return bool(tags_dict.get('name')) except (json.JSONDecodeError, TypeError): return False return False # Check for direct name properties name = row.get('name') or row.get('NAME') or row.get('FULLNAME') return bool(name) def load_geojson(self, filepath: str, filter_unnamed: bool = False) -> gpd.GeoDataFrame: """Load and validate GeoJSON file with optimizations.""" try: # Use pyogr engine for faster loading of large files gdf = gpd.read_file(filepath, engine='pyogrio') # Filter only LineString, MultiLineString, and Polygon geometries line_types = ['LineString', 'MultiLineString', 'Polygon'] gdf = gdf[gdf.geometry.type.isin(line_types)].copy() if len(gdf) == 0: raise ValueError(f"No line geometries found in {filepath}") # Reset index for efficient processing gdf = gdf.reset_index(drop=True) # Ensure geometry is valid and fix simple issues invalid_mask = ~gdf.geometry.is_valid if invalid_mask.any(): print(f"Fixing {invalid_mask.sum()} invalid geometries...") gdf.loc[invalid_mask, 'geometry'] = gdf.loc[invalid_mask, 'geometry'].buffer(0) # Filter unnamed features if requested if filter_unnamed: original_count = len(gdf) named_mask = gdf.apply(self._has_name, axis=1) gdf = gdf[named_mask].copy() gdf = gdf.reset_index(drop=True) filtered_count = original_count - len(gdf) print(f"Filtered out {filtered_count} unnamed features") print(f"Loaded {len(gdf)} road features from {filepath}") return gdf except Exception as e: raise Exception(f"Error loading {filepath}: {str(e)}") def create_buffer_index(self, gdf: gpd.GeoDataFrame) -> Tuple[np.ndarray, STRtree]: """Buffer every geometry and index the buffers. Each road is later compared only against the union of the buffers that touch it (see _local_cover), which gives the same difference as using one county-wide union - buffers that don't touch a line can't change its difference - without every comparison paying for the size of the whole county.""" print("Buffering and indexing geometries...") buffers = np.asarray(gdf.geometry.buffer(self.tolerance_deg).values, dtype=object) return buffers, STRtree(buffers) @staticmethod def _local_cover(geom, buffers: np.ndarray, tree: STRtree): """Union of the target buffers that touch geom, or None if none do.""" idx = tree.query(geom, predicate='intersects') if len(idx) == 0: return None if len(idx) == 1: return buffers[idx[0]] return unary_union(list(buffers[idx])) def create_spatial_index(self, gdf: gpd.GeoDataFrame) -> STRtree: """Create spatial index for fast intersection queries.""" print("Creating spatial index...") # Create STRtree for fast spatial queries geometries = gdf.geometry.tolist() return STRtree(geometries) def find_removed_segments_optimized(self, source_gdf: gpd.GeoDataFrame, target: Tuple[np.ndarray, STRtree]) -> List[Dict[str, Any]]: """ Find segments in source_gdf that are not covered by the target buffers (removed roads). """ print("Finding removed segments...") rows = self._rows_with_cover(source_gdf, target, "removed") return self._process_removed_chunk(rows, self.min_gap_length_deg) def _rows_with_cover(self, gdf: gpd.GeoDataFrame, target: Tuple[np.ndarray, STRtree], label: str) -> List[Dict[str, Any]]: buffers, tree = target geoms = gdf.geometry.tolist() props = gdf.drop(columns='geometry').to_dict('records') rows = [] for i, (geom, p) in enumerate(zip(geoms, props)): rows.append({'geometry': geom, 'properties': p, 'cover': self._local_cover(geom, buffers, tree)}) if (i + 1) % 10000 == 0: print(f" {label}: compared {i + 1}/{len(geoms)}") return rows @staticmethod def _process_removed_chunk(chunk_data: List[Dict], min_length_deg: float) -> List[Dict[str, Any]]: """Process a chunk of geometries to find removed segments.""" removed_segments = [] for row_data in chunk_data: geom = row_data['geometry'] properties = row_data['properties'] # Handle MultiLineString by processing each component if isinstance(geom, MultiLineString): lines = list(geom.geoms) else: lines = [geom] # Polygon and Line can be accessed directly for line in lines: try: # Find parts of the line that don't intersect with target_union cover = row_data['cover'] uncovered = line.difference(cover) if cover is not None else line if uncovered.is_empty: continue # Handle different geometry types returned by difference uncovered_lines = [] if hasattr(uncovered, 'geoms'): for geom_part in uncovered.geoms: if isinstance(geom_part, LineString): uncovered_lines.append(geom_part) elif isinstance(uncovered, LineString): uncovered_lines.append(uncovered) # Check each uncovered line segment for uncovered_line in uncovered_lines: if uncovered_line.length >= min_length_deg: # Create properties dict with original metadata plus 'removed: true' result_properties = properties.copy() result_properties['removed'] = True removed_segments.append({ 'geometry': uncovered_line, **result_properties }) except Exception as e: continue # Skip problematic geometries return removed_segments def find_added_roads_optimized(self, source_gdf: gpd.GeoDataFrame, target: Tuple[np.ndarray, STRtree]) -> List[Dict[str, Any]]: """ Find entire roads in source_gdf that don't significantly overlap with the target buffers. """ print("Finding added roads...") rows = self._rows_with_cover(source_gdf, target, "added") return self._process_added_chunk(rows, self.min_gap_length_deg) @staticmethod def _process_added_chunk(chunk_data: List[Dict], min_length_deg: float) -> List[Dict[str, Any]]: """Process a chunk of geometries to find added roads.""" added_roads = [] for row_data in chunk_data: geom = row_data['geometry'] original_properties = row_data['properties'] try: # Check what portion of the road is not covered cover = row_data['cover'] uncovered = geom.difference(cover) if cover is not None else geom if not uncovered.is_empty: # Calculate what percentage of the original road is uncovered uncovered_length = 0 if hasattr(uncovered, 'geoms'): for geom_part in uncovered.geoms: if isinstance(geom_part, LineString): uncovered_length += geom_part.length elif isinstance(uncovered, LineString): uncovered_length = uncovered.length original_length = geom.length uncovered_ratio = uncovered_length / original_length if original_length > 0 else 0 # Include the entire road if: # 1. The uncovered portion is above minimum threshold, AND # 2. More than 10% of the road is uncovered if uncovered_ratio > 0.1: #uncovered_length >= min_length_deg and # Include entire original road with all original metadata properties = { 'surface': 'asphalt' } # Detect county format based on available fields is_lake_county = 'FullStreet' in original_properties is_sumter_county = 'NAME' in original_properties and 'RoadClass' in original_properties if is_lake_county: # Lake County field mappings for key, value in original_properties.items(): if key == 'FullStreet': properties['name'] = titlecase(qgisfunctions.formatstreet(value,None,None)) if value is not None else None elif key == 'SpeedLimit': properties['maxspeed'] = f"{value} mph" if value is not None else None elif key == 'NumberOfLa': try: num_value = int(float(value)) if value is not None else 0 if num_value > 0: properties['lanes'] = str(num_value) except (ValueError, TypeError): pass elif key == 'StreetClas': highway_type = qgisfunctions.gethighwaytype(value, None, None) properties['highway'] = highway_type if highway_type else 'residential' elif key == 'Oneway': if value == 'FT': properties['oneway'] = 'yes' elif value == 'TF': properties['oneway'] = '-1' elif is_sumter_county: # Sumter County field mappings for key, value in original_properties.items(): if key == 'NAME': properties['name'] = titlecase(qgisfunctions.formatstreet(value,None,None)) if value is not None else None elif key == 'SpeedLimit': properties['maxspeed'] = f"{value} mph" if value is not None else None elif key == 'RoadClass': value_str = str(value) if value is not None else None if value_str is None: properties['highway'] = 'residential' elif value_str.startswith('PRIMARY'): properties['highway'] = 'trunk' elif value_str.startswith('MAJOR'): properties['highway'] = 'primary' else: properties['highway'] = 'residential' elif key == 'OneWayCode': if value == 'F': properties['oneway'] = 'yes' elif value == 'T': properties['oneway'] = '-1' elif key == 'LANES': try: num_value = int(float(value)) if value is not None else 0 if num_value > 0: properties['lanes'] = str(num_value) except (ValueError, TypeError): pass else: # Unknown format - try common field names name = original_properties.get('NAME') or original_properties.get('FullStreet') or original_properties.get('name') if name: properties['name'] = titlecase(qgisfunctions.formatstreet(name,None,None)) speed = original_properties.get('SpeedLimit') if speed: properties['maxspeed'] = f"{speed} mph" properties['highway'] = 'residential' added_roads.append({ 'geometry': geom, **properties }) except Exception as e: print(e) continue # Skip problematic geometries return added_roads def compare_roads(self, file1_path: str, file2_path: str) -> Tuple[List[Dict], List[Dict]]: """ Compare two GeoJSON files and find significant differences. Optimized version with parallel processing. Returns: Tuple of (removed_roads, added_roads) """ print(f"Comparing {file1_path} and {file2_path}") print(f"Tolerance: {self.tolerance_feet} feet") print(f"Minimum significant length: {self.min_gap_length_feet} feet") print("-" * 50) # Load both files # Filter unnamed features from file1 (OSM data) if exclude_unnamed is set gdf1 = self.load_geojson(file1_path, filter_unnamed=self.exclude_unnamed) gdf2 = self.load_geojson(file2_path) # Filter county data (gdf2) to active/current roads only. # Sumter uses LIFECYCLE='ACTV'. Lake has no lifecycle field — skip silently. LIFECYCLE_FIELD = 'LIFECYCLE' ACTIVE_VALUES = {'ACTV', 'Active', 'ACTIVE', 'Current', 'CURRENT'} if LIFECYCLE_FIELD in gdf2.columns: before = len(gdf2) gdf2 = gdf2[gdf2[LIFECYCLE_FIELD].isin(ACTIVE_VALUES)].copy().reset_index(drop=True) filtered = before - len(gdf2) if filtered: print(f"Filtered out {filtered} non-active county roads (LIFECYCLE not in {ACTIVE_VALUES})") else: print(f"All {len(gdf2)} county roads are active (LIFECYCLE={gdf2[LIFECYCLE_FIELD].iloc[0] if len(gdf2) else 'n/a'})") else: print(f"Note: county data has no {LIFECYCLE_FIELD!r} field — no lifecycle filtering applied") # Ensure both are in the same CRS if gdf1.crs != gdf2.crs: print(f"Warning: CRS mismatch. Converting {file2_path} to match {file1_path}") gdf2 = gdf2.to_crs(gdf1.crs) index1 = self.create_buffer_index(gdf1) index2 = self.create_buffer_index(gdf2) removed_roads = self.find_removed_segments_optimized(gdf1, index2) added_roads = self.find_added_roads_optimized(gdf2, index1) # Clean up memory del gdf1, gdf2, index1, index2 gc.collect() return removed_roads, added_roads def save_results(self, removed: List[Dict], added: List[Dict], output_path: str): """Save results to GeoJSON file.""" all_results = removed + added if not all_results: print("No significant differences found!") return # Create GeoDataFrame efficiently print("Saving results...") results_gdf = gpd.GeoDataFrame(all_results) # Save to file with optimization, with fallback for locked files try: results_gdf.to_file(output_path, driver='GeoJSON', engine='pyogrio') print(f"Results saved to: {output_path}") except (PermissionError, OSError) as e: # File is locked, try with a timestamp suffix from datetime import datetime timestamp = datetime.now().strftime("%H%M%S") base = Path(output_path) fallback_path = str(base.parent / f"{base.stem}_{timestamp}{base.suffix}") print(f"Warning: Could not save to {output_path}: {e}") print(f"Saving to fallback: {fallback_path}") results_gdf.to_file(fallback_path, driver='GeoJSON', engine='pyogrio') print(f"Results saved to: {fallback_path}") def print_summary(self, removed: List[Dict], added: List[Dict], file1_name: str, file2_name: str): """Print a summary of the comparison results.""" print("\n" + "="*60) print("COMPARISON SUMMARY") print("="*60) print(f"\nFile 1: {file1_name}") print(f"File 2: {file2_name}") print(f"Tolerance: {self.tolerance_feet} feet") print(f"Minimum significant length: {self.min_gap_length_feet} feet") if removed: print(f"\nREMOVED ROADS ({len(removed)} segments):") print("These road segments exist in File 1 but are missing or incomplete in File 2:") # Calculate total length of removed segments total_removed_length = 0 removed_by_road = {} for segment in removed: geom = segment['geometry'] length_feet = geom.length * 364000.0 # Convert to feet total_removed_length += length_feet # Get road name road_name = "Unknown" name_fields = ['name', 'NAME', 'road_name', 'street_name', 'FULLNAME'] for field in name_fields: if field in segment and pd.notna(segment[field]): road_name = str(segment[field]) break if road_name not in removed_by_road: removed_by_road[road_name] = [] removed_by_road[road_name].append(length_feet) print(f"Total removed length: {total_removed_length:,.1f} feet ({total_removed_length/5280:.2f} miles)") for road, lengths in sorted(removed_by_road.items()): road_total = sum(lengths) print(f" • {road}: {len(lengths)} segment(s), {road_total:,.1f} feet") if added: print(f"\nADDED ROADS ({len(added)} roads):") print("These roads exist in File 2 but are missing or incomplete in File 1:") # Calculate total length of added roads total_added_length = 0 added_by_road = {} for road in added: geom = road['geometry'] length_feet = geom.length * 364000.0 # Convert to feet total_added_length += length_feet # Get road name road_name = "Unknown" name_fields = ['name', 'NAME', 'road_name', 'street_name', 'FULLNAME'] for field in name_fields: if field in road and pd.notna(road[field]): road_name = str(road[field]) break if road_name not in added_by_road: added_by_road[road_name] = 0 added_by_road[road_name] += length_feet print(f"Total added length: {total_added_length:,.1f} feet ({total_added_length/5280:.2f} miles)") for road, length in sorted(added_by_road.items()): print(f" • {road}: {length:,.1f} feet") if not removed and not added: print("\nNo significant differences found!") print("The road networks have good coverage overlap within the specified tolerance.") def main(): parser = argparse.ArgumentParser( description="Compare two GeoJSON files containing roads and find significant gaps or extras (Optimized)", formatter_class=argparse.RawDescriptionHelpFormatter, epilog=""" Examples: python diff_highways.py roads1.geojson roads2.geojson python diff_highways.py roads1.geojson roads2.geojson --tolerance 100 --min-length 200 python diff_highways.py roads1.geojson roads2.geojson --output differences.geojson python diff_highways.py roads1.geojson roads2.geojson --jobs 8 --chunk-size 2000 """ ) parser.add_argument('file1', help='First GeoJSON file') parser.add_argument('file2', help='Second GeoJSON file') parser.add_argument('--tolerance', '-t', type=float, default=50.0, help='Distance tolerance in feet for considering roads as overlapping (default: 50)') parser.add_argument('--min-length', '-m', type=float, default=100.0, help='Minimum length in feet for gaps/extras to be considered significant (default: 100)') parser.add_argument('--output', '-o', help='Output GeoJSON file for results (optional)') parser.add_argument('--jobs', '-j', type=int, default=None, help='Ignored (kept for compatibility)') parser.add_argument('--chunk-size', '-c', type=int, default=1000, help='Ignored (kept for compatibility)') parser.add_argument('--exclude-unnamed', '-e', action='store_true', help='Exclude features without name tags from coverage calculation (helps detect roads covered by unnamed geometry)') args = parser.parse_args() # Validate input files if not Path(args.file1).exists(): print(f"Error: File {args.file1} does not exist") return 1 if not Path(args.file2).exists(): print(f"Error: File {args.file2} does not exist") return 1 try: # Create comparator and run comparison comparator = RoadComparator( tolerance_feet=args.tolerance, min_gap_length_feet=args.min_length, n_jobs=args.jobs, chunk_size=args.chunk_size, exclude_unnamed=args.exclude_unnamed ) removed, added = comparator.compare_roads(args.file1, args.file2) # Print summary comparator.print_summary(removed, added, args.file1, args.file2) # Save results if output file specified if args.output: comparator.save_results(removed, added, args.output) elif removed or added: # Auto-generate output filename if differences found output_file = f"road_differences_{Path(args.file1).stem}_vs_{Path(args.file2).stem}.geojson" comparator.save_results(removed, added, output_file) return 0 except Exception as e: print(f"Error: {str(e)}") return 1 if __name__ == "__main__": exit(main())