Converting CADRG Maps to GeoJSON with Python for Emergency Response Workflows

import os
import json
import logging
import rasterio
from rasterio.features import shapes
from rasterio.transform import from_origin
from shapely.geometry import shape, mapping
from pyproj import Transformer
from pathlib import Path

logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")

def parse_cadrg_toc(toc_path: str) -> str:
    """Extract EPSG code from CADRG .toc header with strict fallback."""
    fallback_epsg = "EPSG:4326"
    if not Path(toc_path).exists():
        logging.warning("TOC missing. Falling back to WGS84. Verify spatial alignment manually.")
        return fallback_epsg

    try:
        with open(toc_path, "r", encoding="utf-8") as f:
            for line in f:
                if "ZONE" in line.upper() or "EPSG" in line.upper():
                    parts = line.replace("=", " ").split()
                    for p in parts:
                        if p.startswith("EPSG:"):
                            return p
        logging.warning("No explicit CRS in TOC. Defaulting to EPSG:4326.")
        return fallback_epsg
    except Exception as e:
        logging.error(f"TOC parse failure: {e}")
        return fallback_epsg

def transform_and_vectorize(raster_path: str, crs_src: str, threshold: int = 0) -> dict:
    """Convert CADRG raster to GeoJSON with memory-safe chunking and CRS transform."""
    try:
        with rasterio.open(raster_path) as src:
            # Validate source CRS
            if src.crs is None:
                src.crs = rasterio.crs.CRS.from_string(crs_src)

            transformer = Transformer.from_crs(src.crs, "EPSG:4326", always_xy=True)

            features = []
            # Process in bands to avoid OOM on large tactical sheets
            for band_idx in range(1, src.count + 1):
                band_data = src.read(band_idx)
                # Mask out background/noise
                mask = band_data > threshold

                for geom, val in shapes(band_data, mask=mask, transform=src.transform):
                    poly = shape(geom)
                    if not poly.is_valid:
                        poly = poly.buffer(0)  # Topology repair fallback

                    # Transform to WGS84
                    coords = list(poly.exterior.coords)
                    transformed_coords = [transformer.transform(x, y) for x, y in coords]
                    poly_wgs84 = shape({"type": "Polygon", "coordinates": [transformed_coords]})

                    features.append({
                        "type": "Feature",
                        "geometry": mapping(poly_wgs84),
                        "properties": {
                            "band": band_idx,
                            "pixel_value": int(val),
                            "source_crs": crs_src,
                            "processing_engine": "rasterio_pyproj"
                        }
                    })

            return {
                "type": "FeatureCollection",
                "features": features,
                "crs": {"type": "name", "properties": {"name": "urn:ogc:def:crs:OGC:1.3:CRS84"}}
            }
    except MemoryError:
        logging.critical("Memory threshold exceeded. Switching to GDAL CLI fallback.")
        raise RuntimeError("Use gdal_polygonize.py for datasets >2GB on constrained edge nodes.")
    except Exception as e:
        logging.error(f"Vectorization failed: {e}")
        raise

def export_geojson(geojson_data: dict, output_path: str, metadata: dict):
    """Attach ISO-compliant metadata and write to disk."""
    geojson_data["metadata"] = metadata
    with open(output_path, "w", encoding="utf-8") as f:
        json.dump(geojson_data, f, indent=2)
    logging.info(f"GeoJSON exported to {output_path}")

# Execution wrapper
if __name__ == "__main__":
    RASTER = "cadrg_sheet.tif"
    TOC = "cadrg_sheet.toc"
    OUT = "incident_zone.geojson"

    src_crs = parse_cadrg_toc(TOC)
    fc = transform_and_vectorize(RASTER, src_crs, threshold=1)
    export_geojson(fc, OUT, {
        "standard": "ISO 19115",
        "acquisition_utc": "2024-10-15T08:00:00Z",
        "source_edition": "CADRG v3.2",
        "crs_transform": f"{src_crs} -> EPSG:4326"
    })