zyphlar
305 lines
12 KiB
Python
305 lines
12 KiB
Python
#!/usr/bin/env python3
|
|
"""
|
|
Compare OSM-formatted county address GeoJSON against OpenStreetMap addresses.
|
|
|
|
Expects both inputs to already be in OSM field format (addr:housenumber, addr:street, etc.).
|
|
Use convert-addresses.py to prepare the county file and download-overpass.py for the OSM file.
|
|
|
|
Usage:
|
|
python compare-addresses.py \
|
|
--local-file /data/sumter/county-addresses.geojson \
|
|
--osm-file /data/sumter/osm-addresses.geojson \
|
|
--output-dir /data/sumter
|
|
"""
|
|
|
|
import argparse
|
|
import re
|
|
import subprocess
|
|
import sys
|
|
from pathlib import Path
|
|
from typing import Dict, List, Tuple
|
|
|
|
import geopandas as gpd
|
|
import pandas as pd
|
|
from shapely.strtree import STRtree
|
|
import warnings
|
|
|
|
warnings.filterwarnings('ignore')
|
|
|
|
|
|
CONFLICT_RADIUS_METERS = 5.0
|
|
|
|
|
|
class AddressComparator:
|
|
def __init__(self, tolerance_meters: float = 500.0):
|
|
self.tolerance_meters = tolerance_meters
|
|
# 1 degree latitude ≈ 111,000 metres
|
|
self.tolerance_deg = tolerance_meters / 111000.0
|
|
|
|
@staticmethod
|
|
def _normalize_multi(value) -> str:
|
|
"""Normalise a semicolon-delimited field by sorting its parts.
|
|
|
|
Treats '101;201' and '201;101' as equal so ordering differences
|
|
don't create false positives.
|
|
"""
|
|
if value is None or (isinstance(value, float) and pd.isna(value)):
|
|
return ''
|
|
s = str(value).strip()
|
|
if not s or s == 'nan':
|
|
return ''
|
|
parts = sorted(p.strip().lower() for p in s.split(';') if p.strip())
|
|
return ';'.join(parts)
|
|
|
|
def _normalize_street_name(self, street: str) -> str:
|
|
"""Normalise street names for fuzzy matching across county/OSM formatting differences."""
|
|
if not street or street == 'nan':
|
|
return ''
|
|
s = street.strip().lower()
|
|
s = re.sub(r'\bstate road\b', 'sr', s)
|
|
s = re.sub(r'\bstate route\b', 'sr', s)
|
|
s = re.sub(r'\bcounty road\b', 'cr', s)
|
|
s = re.sub(r'\bc\b', 'cr', s) # "C 44a" -> "cr 44a"
|
|
# normalise CR sub-segment separators: "cr 109d 1" -> "cr 109d-1"
|
|
s = re.sub(r'\b(cr\s+\d+[a-z])\s+(\d+)', r'\1-\2', s)
|
|
s = re.sub(r'\s+', ' ', s).strip()
|
|
return s
|
|
|
|
@staticmethod
|
|
def _flag_close_conflicts(addresses: List[Dict], radius_meters: float = CONFLICT_RADIUS_METERS) -> None:
|
|
"""Mark each address dict's 'conflict' key True when another point in the
|
|
same list falls within radius_meters. Doesn't dedupe or exclude anything -
|
|
just surfaces tight clusters (e.g. ambiguous source data reusing the same
|
|
unit label across distinct points) for manual review in the map UI."""
|
|
if not addresses:
|
|
return
|
|
geoms = [a['geometry'] for a in addresses]
|
|
tree = STRtree(geoms)
|
|
radius_deg = radius_meters / 111000.0
|
|
for i, addr in enumerate(addresses):
|
|
conflict = False
|
|
for j in tree.query(geoms[i].buffer(radius_deg)):
|
|
if j == i:
|
|
continue
|
|
if geoms[i].distance(geoms[j]) * 111000.0 <= radius_meters:
|
|
conflict = True
|
|
break
|
|
addr['conflict'] = conflict
|
|
|
|
def compare_addresses(self, local_file: str, osm_file: str) -> Tuple[List[Dict], List[Dict], List[Dict]]:
|
|
"""Compare local and OSM address data.
|
|
|
|
Returns:
|
|
(new_addresses, existing_addresses, removed_addresses)
|
|
"""
|
|
print(f"Comparing addresses: {local_file} vs {osm_file}")
|
|
|
|
local_gdf = gpd.read_file(local_file)
|
|
osm_gdf = gpd.read_file(osm_file)
|
|
|
|
print(f"Loaded local addresses: {len(local_gdf)}")
|
|
print(f"Loaded OSM addresses: {len(osm_gdf)}")
|
|
|
|
if hasattr(self, 'sample_size') and self.sample_size:
|
|
if len(local_gdf) > self.sample_size:
|
|
local_gdf = local_gdf.sample(n=self.sample_size, random_state=42).reset_index(drop=True)
|
|
print(f"Sampled local addresses to: {len(local_gdf)}")
|
|
|
|
if hasattr(self, 'max_osm') and self.max_osm:
|
|
if len(osm_gdf) > self.max_osm:
|
|
osm_gdf = osm_gdf.sample(n=self.max_osm, random_state=42).reset_index(drop=True)
|
|
print(f"Sampled OSM addresses to: {len(osm_gdf)}")
|
|
|
|
if local_gdf.crs != osm_gdf.crs:
|
|
osm_gdf = osm_gdf.to_crs(local_gdf.crs)
|
|
|
|
osm_index = STRtree(osm_gdf.geometry.tolist())
|
|
|
|
existing_addresses = []
|
|
new_addresses = []
|
|
|
|
for _, local_row in local_gdf.iterrows():
|
|
local_point = local_row.geometry
|
|
nearby = osm_index.query(local_point.buffer(self.tolerance_deg))
|
|
|
|
best_match = None
|
|
min_distance = float('inf')
|
|
|
|
for osm_idx in nearby:
|
|
osm_row = osm_gdf.iloc[osm_idx]
|
|
distance = local_point.distance(osm_row.geometry)
|
|
|
|
local_num = self._normalize_multi(local_row.get('addr:housenumber', ''))
|
|
osm_num = self._normalize_multi(osm_row.get('addr:housenumber', ''))
|
|
|
|
local_street = self._normalize_street_name(str(local_row.get('addr:street', '')))
|
|
osm_street = self._normalize_street_name(str(osm_row.get('addr:street', '')))
|
|
street_match = local_street == osm_street and local_street != ''
|
|
|
|
local_unit = self._normalize_multi(local_row.get('addr:unit'))
|
|
osm_unit = self._normalize_multi(osm_row.get('addr:unit'))
|
|
|
|
if local_unit and osm_unit:
|
|
unit_match = local_unit == osm_unit
|
|
elif local_unit or osm_unit:
|
|
unit_match = False
|
|
else:
|
|
unit_match = True
|
|
|
|
if local_num == osm_num and street_match and unit_match and distance < min_distance:
|
|
min_distance = distance
|
|
best_match = osm_idx
|
|
|
|
distance_m = min_distance * 111000.0
|
|
|
|
if best_match is not None and distance_m <= self.tolerance_meters:
|
|
existing_addresses.append({
|
|
'geometry': local_point,
|
|
'local_data': dict(local_row.drop('geometry')),
|
|
'osm_data': dict(osm_gdf.iloc[best_match].drop('geometry')),
|
|
'distance_meters': distance_m,
|
|
})
|
|
else:
|
|
props = dict(local_row.drop('geometry'))
|
|
props['status'] = 'new'
|
|
new_addresses.append({'geometry': local_point, **props})
|
|
|
|
# Second pass: find OSM addresses with no local match (potentially removed)
|
|
matched_osm = set()
|
|
for _, local_row in local_gdf.iterrows():
|
|
local_point = local_row.geometry
|
|
nearby = osm_index.query(local_point.buffer(self.tolerance_deg))
|
|
|
|
local_num = self._normalize_multi(local_row.get('addr:housenumber', ''))
|
|
local_street = self._normalize_street_name(str(local_row.get('addr:street', '')))
|
|
local_unit = self._normalize_multi(local_row.get('addr:unit'))
|
|
|
|
for osm_idx in nearby:
|
|
osm_row = osm_gdf.iloc[osm_idx]
|
|
dist_m = local_point.distance(osm_row.geometry) * 111000.0
|
|
osm_num = self._normalize_multi(osm_row.get('addr:housenumber', ''))
|
|
osm_st = self._normalize_street_name(str(osm_row.get('addr:street', '')))
|
|
osm_unit = self._normalize_multi(osm_row.get('addr:unit'))
|
|
|
|
unit_match = True
|
|
if local_unit and osm_unit:
|
|
unit_match = local_unit == osm_unit
|
|
|
|
if (dist_m <= self.tolerance_meters and
|
|
local_num == osm_num and
|
|
local_street == osm_st and local_street != '' and
|
|
unit_match):
|
|
matched_osm.add(osm_idx)
|
|
break
|
|
|
|
removed_addresses = []
|
|
for idx, osm_row in osm_gdf.iterrows():
|
|
if idx not in matched_osm:
|
|
props = dict(osm_row.drop('geometry'))
|
|
props['status'] = 'removed'
|
|
removed_addresses.append({'geometry': osm_row.geometry, **props})
|
|
|
|
self._flag_close_conflicts(new_addresses)
|
|
|
|
return new_addresses, existing_addresses, removed_addresses
|
|
|
|
def save_results(self, new_addresses, existing_addresses, removed_addresses, output_dir):
|
|
output_dir = Path(output_dir)
|
|
output_dir.mkdir(parents=True, exist_ok=True)
|
|
|
|
if new_addresses:
|
|
gpd.GeoDataFrame(new_addresses).to_file(
|
|
output_dir / "addresses-to-add.geojson", driver='GeoJSON')
|
|
print(f"Saved {len(new_addresses)} new addresses to {output_dir}/addresses-to-add.geojson")
|
|
|
|
if removed_addresses:
|
|
gpd.GeoDataFrame(removed_addresses).to_file(
|
|
output_dir / "addresses-potentially-removed.geojson", driver='GeoJSON')
|
|
print(f"Saved {len(removed_addresses)} potentially removed to {output_dir}/addresses-potentially-removed.geojson")
|
|
|
|
if existing_addresses:
|
|
existing_simple = [
|
|
{'geometry': a['geometry'], 'distance_meters': a['distance_meters'], 'status': 'existing'}
|
|
for a in existing_addresses
|
|
]
|
|
gpd.GeoDataFrame(existing_simple).to_file(
|
|
output_dir / "addresses-existing.geojson", driver='GeoJSON')
|
|
print(f"Saved {len(existing_addresses)} existing addresses to {output_dir}/addresses-existing.geojson")
|
|
|
|
def print_summary(self, new_addresses, existing_addresses, removed_addresses):
|
|
print("\n" + "=" * 60)
|
|
print("ADDRESS COMPARISON SUMMARY")
|
|
print("=" * 60)
|
|
print(f"\n New (to add to OSM): {len(new_addresses)}")
|
|
print(f" Existing (matched): {len(existing_addresses)}")
|
|
print(f" Potentially removed: {len(removed_addresses)}")
|
|
|
|
if existing_addresses:
|
|
distances = [a['distance_meters'] for a in existing_addresses]
|
|
print(f"\n Avg match distance: {sum(distances)/len(distances):.1f} m")
|
|
print(f" Max match distance: {max(distances):.1f} m")
|
|
|
|
if new_addresses:
|
|
print(f"\n Sample new addresses by street:")
|
|
streets: Dict[str, int] = {}
|
|
for addr in new_addresses[:10]:
|
|
s = addr.get('addr:street', 'Unknown')
|
|
streets[s] = streets.get(s, 0) + 1
|
|
for s, n in sorted(streets.items()):
|
|
print(f" {s}: {n}")
|
|
if len(new_addresses) > 10:
|
|
print(f" ... and {len(new_addresses) - 10} more")
|
|
|
|
|
|
def main():
|
|
parser = argparse.ArgumentParser(
|
|
description="Compare OSM-formatted county addresses against OpenStreetMap",
|
|
epilog="""
|
|
Run convert-addresses.py first to produce county-addresses.geojson,
|
|
and download-overpass.py to produce osm-addresses.geojson.
|
|
"""
|
|
)
|
|
parser.add_argument('--local-file', required=True,
|
|
help='Pre-converted county addresses GeoJSON (from convert-addresses.py)')
|
|
parser.add_argument('--osm-file', required=True,
|
|
help='Downloaded OSM addresses GeoJSON (from download-overpass.py)')
|
|
parser.add_argument('--output-dir', '-o', required=True,
|
|
help='Output directory for result GeoJSON files')
|
|
parser.add_argument('--tolerance', '-t', type=float, default=500.0,
|
|
help='Spatial tolerance in metres for matching (default: 500)')
|
|
parser.add_argument('--sample', '-s', type=int,
|
|
help='Process only N local addresses (for testing)')
|
|
parser.add_argument('--max-osm', type=int,
|
|
help='Cap OSM addresses loaded (for testing)')
|
|
|
|
args = parser.parse_args()
|
|
|
|
local_file = Path(args.local_file)
|
|
osm_file = Path(args.osm_file)
|
|
output_dir = Path(args.output_dir)
|
|
|
|
if not local_file.exists():
|
|
print(f"Error: local file not found: {local_file}", file=sys.stderr)
|
|
print("Run convert-addresses.py first.", file=sys.stderr)
|
|
return 1
|
|
|
|
if not osm_file.exists():
|
|
print(f"Error: OSM file not found: {osm_file}", file=sys.stderr)
|
|
print("Run download-overpass.py first.", file=sys.stderr)
|
|
return 1
|
|
|
|
comparator = AddressComparator(tolerance_meters=args.tolerance)
|
|
if args.sample:
|
|
comparator.sample_size = args.sample
|
|
if args.max_osm:
|
|
comparator.max_osm = args.max_osm
|
|
|
|
new, existing, removed = comparator.compare_addresses(str(local_file), str(osm_file))
|
|
comparator.save_results(new, existing, removed, output_dir)
|
|
comparator.print_summary(new, existing, removed)
|
|
return 0
|
|
|
|
|
|
if __name__ == '__main__':
|
|
sys.exit(main())
|