worlddiff.py

#!/usr/bin/env python3
"""Exact, dependency-free change detection for GeoJSON FeatureCollections."""

from __future__ import annotations

import argparse
import json
from pathlib import Path
from typing import Any


def load_features(path: str | Path, key: str = "id") -> dict[str, dict[str, Any]]:
    data = json.loads(Path(path).read_text(encoding="utf-8"))
    if data.get("type") != "FeatureCollection" or not isinstance(data.get("features"), list):
        raise ValueError(f"{path} is not a GeoJSON FeatureCollection")
    indexed: dict[str, dict[str, Any]] = {}
    for position, feature in enumerate(data["features"], 1):
        raw_id = feature.get("id") if key == "id" else feature.get("properties", {}).get(key)
        if raw_id is None:
            raise ValueError(f"feature {position} in {path} has no {key!r}")
        feature_id = str(raw_id)
        if feature_id in indexed:
            raise ValueError(f"duplicate {key}={feature_id!r} in {path}")
        indexed[feature_id] = feature
    return indexed


def compare(old: dict[str, dict[str, Any]], new: dict[str, dict[str, Any]]) -> dict[str, Any]:
    old_ids, new_ids = set(old), set(new)
    changed = []
    for feature_id in sorted(old_ids & new_ids):
        before, after = old[feature_id], new[feature_id]
        old_props = before.get("properties") or {}
        new_props = after.get("properties") or {}
        fields = sorted(
            name for name in set(old_props) | set(new_props) if old_props.get(name) != new_props.get(name)
        )
        geometry_changed = before.get("geometry") != after.get("geometry")
        if fields or geometry_changed:
            changed.append(
                {
                    "id": feature_id,
                    "geometry_changed": geometry_changed,
                    "properties_changed": fields,
                }
            )
    return {
        "added": sorted(new_ids - old_ids),
        "removed": sorted(old_ids - new_ids),
        "changed": changed,
        "unchanged_count": len(old_ids & new_ids) - len(changed),
    }


def change_layer(
    old: dict[str, dict[str, Any]], new: dict[str, dict[str, Any]], summary: dict[str, Any]
) -> dict[str, Any]:
    features = []
    selections = (
        ("added", summary["added"], new),
        ("removed", summary["removed"], old),
        ("changed", [item["id"] for item in summary["changed"]], new),
    )
    for label, ids, source in selections:
        for feature_id in ids:
            feature = json.loads(json.dumps(source[feature_id]))
            feature.setdefault("properties", {})["_worlddiff_change"] = label
            features.append(feature)
    return {"type": "FeatureCollection", "features": features}


def svg_from_layer(layer: dict[str, Any]) -> str:
    """Render a change_layer FeatureCollection as an SVG map."""
    features = layer.get("features", [])
    if not features:
        return '<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 800 600"><text x="400" y="300" text-anchor="middle">No features</text></svg>'

    # Collect all coordinates to compute bbox
    all_coords = []
    for feature in features:
        geom = feature.get("geometry", {})
        coords = geom.get("coordinates", [])
        if geom.get("type") == "Point":
            all_coords.append(coords)
        elif geom.get("type") == "LineString":
            all_coords.extend(coords)
        elif geom.get("type") == "Polygon":
            for ring in coords:
                all_coords.extend(ring)

    if not all_coords:
        return '<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 800 600"><text x="400" y="300" text-anchor="middle">No coordinates</text></svg>'

    # Compute bbox
    lons = [c[0] for c in all_coords]
    lats = [c[1] for c in all_coords]
    min_lon, max_lon = min(lons), max(lons)
    min_lat, max_lat = min(lats), max(lats)

    # Add 5% padding
    lon_range = max_lon - min_lon or 0.01
    lat_range = max_lat - min_lat or 0.01
    padding = 0.05
    min_lon -= lon_range * padding
    max_lon += lon_range * padding
    min_lat -= lat_range * padding
    max_lat += lat_range * padding

    # SVG viewport
    svg_width, svg_height = 800, 600

    def project(lon: float, lat: float) -> tuple[float, float]:
        """Equirectangular projection: lon/lat -> pixel coords."""
        x = (lon - min_lon) / (max_lon - min_lon) * svg_width
        y = (1 - (lat - min_lat) / (max_lat - min_lat)) * svg_height
        return x, y

    def coords_to_path(coords: list, geom_type: str) -> str:
        """Convert coordinates to SVG path data."""
        if geom_type == "Point":
            x, y = project(coords[0], coords[1])
            return f"M{x},{y}m-3,0a3,3 0 1,0 6,0a3,3 0 1,0 -6,0"
        elif geom_type == "LineString":
            path_data = "M" + " L".join(f"{project(c[0], c[1])[0]},{project(c[0], c[1])[1]}" for c in coords)
            return path_data
        elif geom_type == "Polygon":
            # First ring is outer boundary
            path_data = "M" + " L".join(f"{project(c[0], c[1])[0]},{project(c[0], c[1])[1]}" for c in coords[0])
            path_data += " Z"
            return path_data
        return ""

    # Build paths for each feature
    paths = []
    for feature in features:
        geom = feature.get("geometry", {})
        change_type = feature.get("properties", {}).get("_worlddiff_change", "unchanged")
        geom_type = geom.get("type", "")
        coords = geom.get("coordinates", [])

        if not coords or not geom_type:
            continue

        path_data = coords_to_path(coords, geom_type)
        if path_data:
            paths.append(f'  <path class="{change_type}" d="{path_data}"/>')

    # Build SVG
    svg = f'''<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 {svg_width} {svg_height}" width="{svg_width}" height="{svg_height}">
  <style>
    .added {{ fill: #22c55e; stroke: #000; stroke-width: 0.5; opacity: 0.8; }}
    .removed {{ fill: #ef4444; stroke: #000; stroke-width: 0.5; opacity: 0.8; }}
    .changed {{ fill: #f97316; stroke: #000; stroke-width: 0.5; opacity: 0.8; }}
    .unchanged {{ fill: #9ca3af; stroke: #000; stroke-width: 0.5; opacity: 0.3; }}
  </style>
{chr(10).join(paths)}
</svg>'''
    return svg


def run(old_path: str | Path, new_path: str | Path, key: str, out: str | Path, svg_path: str | Path | None = None) -> dict[str, Any]:
    old, new = load_features(old_path, key), load_features(new_path, key)
    summary = compare(old, new)
    output = Path(out)
    output.mkdir(parents=True, exist_ok=True)
    (output / "summary.json").write_text(json.dumps(summary, indent=2) + "\n", encoding="utf-8")
    layer = change_layer(old, new, summary)
    (output / "changes.geojson").write_text(json.dumps(layer, indent=2) + "\n", encoding="utf-8")
    if svg_path:
        svg_file = Path(svg_path)
        svg_file.parent.mkdir(parents=True, exist_ok=True)
        svg_file.write_text(svg_from_layer(layer), encoding="utf-8")
    return summary


def main() -> None:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("old")
    parser.add_argument("new")
    parser.add_argument("--key", default="id", help="top-level id or feature property name")
    parser.add_argument("--out", default="worlddiff-report")
    parser.add_argument("--svg", help="write change layer as SVG to this path")
    args = parser.parse_args()
    summary = run(args.old, args.new, args.key, args.out, svg_path=args.svg)
    print(
        f"{len(summary['added'])} added, {len(summary['removed'])} removed, "
        f"{len(summary['changed'])} changed"
    )


if __name__ == "__main__":
    main()