From ae986526b1566b378eed1cff22e34cdbfea4128e Mon Sep 17 00:00:00 2001 From: amit-spatial Date: Tue, 7 Jul 2026 04:45:55 +0000 Subject: [PATCH 1/3] Add reproducible CS admin standard builder This script is the source-of-truth builder for the Core Stack village admin resource. It reads raw admin boundary inputs directly, records source trace and omission decisions, and emits the standardized GeoPackage, GeoJSON, CSV, and QA reports needed by downstream facilities, Antyodaya, and livestock assets. --- .../build_cs_admin_boundary_standard.py | 2546 +++++++++++++++++ 1 file changed, 2546 insertions(+) create mode 100644 utilities/scripts/admin_assets/build_cs_admin_boundary_standard.py diff --git a/utilities/scripts/admin_assets/build_cs_admin_boundary_standard.py b/utilities/scripts/admin_assets/build_cs_admin_boundary_standard.py new file mode 100644 index 00000000..64cff071 --- /dev/null +++ b/utilities/scripts/admin_assets/build_cs_admin_boundary_standard.py @@ -0,0 +1,2546 @@ +#!/usr/bin/env python3 +"""Build the Core Stack standardized village admin boundary asset. + +This script is intentionally separate from ``build_admin_boundary_assets.py``. +It reads the district GeoJSON files under ``data/admin-boundary/input`` and +keeps a source-row trace from those raw files into the final standard asset. + +Default outputs: + + data/admin-boundary/cs_admin_standard.gpkg + data/admin-boundary/cs_admin_standard.geojson + data/admin-boundary/cs_admin_standard.csv + data/admin-boundary/cs_admin_standard_reports/* + +The GeoPackage remains the primary output. The CSV is an attribute-only sibling +for review, joins, and quick QA. +""" + +from __future__ import annotations + +import argparse +from collections import Counter +import csv +import hashlib +import json +import math +from pathlib import Path +import resource +import sqlite3 +import tempfile +import time +from typing import Any, Iterable, Sequence + +try: + import geopandas as gpd + import pandas as pd + import pyogrio +except Exception as exc: # pragma: no cover - handled in ensure_dependencies + gpd = None + pd = None + pyogrio = None + GEOSTACK_IMPORT_ERROR = exc +else: + GEOSTACK_IMPORT_ERROR = None + +try: + from pyproj import Geod +except Exception as exc: # pragma: no cover - handled in ensure_dependencies + Geod = None + PYPROJ_IMPORT_ERROR = exc +else: + PYPROJ_IMPORT_ERROR = None + +try: + from shapely.geometry import MultiPolygon + from shapely.ops import unary_union + from shapely.wkb import loads as load_wkb +except Exception as exc: # pragma: no cover - handled in ensure_dependencies + MultiPolygon = None + unary_union = None + load_wkb = None + SHAPELY_IMPORT_ERROR = exc +else: + SHAPELY_IMPORT_ERROR = None + +try: + from shapely.validation import make_valid as shapely_make_valid +except Exception: # pragma: no cover - compatibility fallback + shapely_make_valid = None + + +ROOT_DIR = Path(__file__).resolve().parents[3] + +DEFAULT_ADMIN_INPUT_DIR = Path("data/admin-boundary/input") +DEFAULT_MULTIPART_ANALYSIS = Path("data/admin-boundary/multi_part_village_analysis.csv") +DEFAULT_OUTPUT_GPKG = Path("data/admin-boundary/cs_admin_standard.gpkg") +DEFAULT_OUTPUT_LAYER = "cs_admin_standard" +DEFAULT_OUTPUT_GEOJSON = Path("data/admin-boundary/cs_admin_standard.geojson") +DEFAULT_OUTPUT_CSV = Path("data/admin-boundary/cs_admin_standard.csv") +DEFAULT_REPORTS_DIR = Path("data/admin-boundary/cs_admin_standard_reports") +DEFAULT_CHUNK_SIZE = 50_000 +SQLITE_BATCH_SIZE = 5_000 +GEOMETRY_OVERLAP_AREA_EPSILON = 0.0 + +ADMIN_COLUMNS = [ + "state_name", + "district_name", + "TEHSIL", + "pc11_village_id", + "NAME", + "pc11_state_id", + "pc11_district_id", + "pc11_subdistrict_id", +] + +ADMIN_KEY_COLUMNS = ["state_name", "district_name", "TEHSIL", "pc11_village_id"] +ID_COLUMNS = { + "pc11_village_id", + "pc11_state_id", + "pc11_district_id", + "pc11_subdistrict_id", +} + +ADMIN_INT_COLUMNS = [ + "feature_part_index", + "feature_part_count", + "pc11_village_id", + "pc11_state_id", + "pc11_district_id", + "pc11_subdistrict_id", + "source_row_count", + "source_identity_count", + "source_geometry_count", + "source_file_count", + "standard_group_part_count", + "standard_name_count", +] + +STANDARD_COLUMNS = [ + "cs_feature_id", + "cs_admin_uid", + "core_admin_uid", + "source_group_key", + "standard_group_key", + "feature_part_index", + "feature_part_count", + "village_shape_hash", + "state_name", + "district_name", + "TEHSIL", + "pc11_village_id", + "village_id", + "NAME", + "pc11_state_id", + "pc11_district_id", + "pc11_subdistrict_id", + "source_row_count", + "source_identity_count", + "source_geometry_count", + "source_file_count", + "resolution_status", + "standard_action", + "standard_reason", + "standard_notes", + "standard_pattern_flags", + "standard_group_part_count", + "standard_part_area_km2", + "standard_part_area_share", + "standard_group_area_km2", + "standard_min_edge_distance_m", + "standard_min_centroid_distance_km", + "standard_max_centroid_distance_km", + "standard_name_count", +] + +GEOD = Geod(ellps="WGS84") if Geod is not None else None + + +def repo_path(path: Path | str) -> Path: + path = Path(path) + return path if path.is_absolute() else ROOT_DIR / path + + +def ensure_dependencies() -> None: + missing = [] + if gpd is None or pd is None or pyogrio is None: + missing.append(f"geopandas/pandas/pyogrio ({GEOSTACK_IMPORT_ERROR})") + if Geod is None: + missing.append(f"pyproj ({PYPROJ_IMPORT_ERROR})") + if MultiPolygon is None or unary_union is None or load_wkb is None: + missing.append(f"shapely ({SHAPELY_IMPORT_ERROR})") + if missing: + raise SystemExit( + "Missing geospatial dependencies: " + + ", ".join(missing) + + "\nRun with uv, for example: " + + "`uv run --with geopandas --with pyogrio --with shapely " + + "--with pyproj --with pandas python " + + "utilities/scripts/admin_assets/build_cs_admin_boundary_standard.py --overwrite`." + ) + + +def clean_text(value: Any) -> str | None: + if value is None: + return None + if isinstance(value, float) and math.isnan(value): + return None + text = str(value).strip() + if text.lower() in {"", "nan", "none", "null", ""}: + return None + return text + + +def clean_int(value: Any) -> int | None: + if value is None: + return None + if isinstance(value, bool): + return int(value) + if isinstance(value, int): + return value if value != 0 else None + if isinstance(value, float): + if math.isnan(value) or value == 0: + return None + if value.is_integer(): + return int(value) + text = str(value).strip() + if text.lower() in {"", "0", "0.0", "nan", "none", "null", ""}: + return None + if text.endswith(".0") and text[:-2].isdigit(): + text = text[:-2] + return int(text) if text.isdigit() else None + + +def normalize_label(value: Any) -> str: + text = clean_text(value) + return "" if text is None else " ".join(text.upper().split()) + + +def normalize_slug(value: Any) -> str: + text = clean_text(value) + if text is None: + return "unknown" + output = [] + last_was_sep = False + for char in text.lower(): + if char.isalnum(): + output.append(char) + last_was_sep = False + elif not last_was_sep: + output.append("_") + last_was_sep = True + slug = "".join(output).strip("_") + return slug or "unknown" + + +def quote_identifier(value: str) -> str: + return '"' + value.replace('"', '""') + '"' + + +def json_dumps(payload: Any, *, indent: bool = False) -> str: + return json.dumps( + payload, + indent=2 if indent else None, + ensure_ascii=False, + sort_keys=not indent, + default=str, + ) + + +def stable_shape_hash(geometry: Any) -> str | None: + if geometry is None or geometry.is_empty: + return None + normalised = geometry.normalize() if hasattr(geometry, "normalize") else geometry + return hashlib.sha256(normalised.wkb).hexdigest() + + +def normalize_property(column: str, value: Any) -> int | str | None: + if column in ID_COLUMNS: + return clean_int(value) + return clean_text(value) + + +def stable_identity_hash(properties: dict[str, Any], geometry_hash: str | None) -> str: + payload = { + "properties": {column: properties.get(column) for column in ADMIN_COLUMNS}, + "geometry_hash": geometry_hash, + } + return hashlib.sha256(json_dumps(payload).encode("utf-8")).hexdigest() + + +def source_group_key_for(properties: dict[str, Any], identity_hash: str) -> str: + village_id = properties.get("pc11_village_id") + if village_id is None: + return f"invalid:{identity_hash[:24]}" + parts = [normalize_slug(properties.get(column)) for column in ADMIN_KEY_COLUMNS[:-1]] + return "village:" + ":".join([*parts, str(village_id)]) + + +def repair_geometry(geometry: Any) -> Any | None: + if geometry is None or geometry.is_empty: + return None + if not geometry.is_valid: + geometry = shapely_make_valid(geometry) if shapely_make_valid is not None else geometry.buffer(0) + if geometry is None or geometry.is_empty: + return None + if geometry.geom_type == "Polygon": + return MultiPolygon([geometry]) + if geometry.geom_type == "MultiPolygon": + return geometry + if geometry.geom_type == "GeometryCollection": + polygons = [] + for part in geometry.geoms: + if part.geom_type == "Polygon": + polygons.append(part) + elif part.geom_type == "MultiPolygon": + polygons.extend(part.geoms) + return MultiPolygon(polygons) if polygons else None + return None + + +def geodesic_area_km2(geometry: Any) -> float: + if geometry is None or geometry.is_empty or GEOD is None: + return 0.0 + try: + area_m2, _ = GEOD.geometry_area_perimeter(geometry) + except Exception: + repaired = repair_geometry(geometry) + if repaired is None: + return 0.0 + area_m2, _ = GEOD.geometry_area_perimeter(repaired) + return abs(float(area_m2)) / 1_000_000.0 + + +def local_utm_epsg(geometries: Sequence[Any]) -> int: + non_empty = [geom for geom in geometries if geom is not None and not geom.is_empty] + if not non_empty: + return 32643 + centroid = unary_union(non_empty).centroid + zone = int((centroid.x + 180) / 6) + 1 + return 32600 + max(1, min(zone, 60)) + + +def group_distance_metrics(geometries: Sequence[Any]) -> dict[str, float | None]: + if len(geometries) < 2: + return { + "min_edge_distance_m": None, + "min_centroid_distance_km": None, + "max_centroid_distance_km": None, + } + repaired = [repair_geometry(geometry) for geometry in geometries] + repaired = [geometry for geometry in repaired if geometry is not None] + if len(repaired) < 2: + return { + "min_edge_distance_m": None, + "min_centroid_distance_km": None, + "max_centroid_distance_km": None, + } + projected = gpd.GeoSeries(repaired, crs="EPSG:4326").to_crs(epsg=local_utm_epsg(repaired)) + edge_distances = [] + centroid_distances = [] + for left_index, left_geometry in enumerate(projected): + for right_geometry in projected.iloc[left_index + 1 :]: + edge_distances.append(float(left_geometry.distance(right_geometry))) + centroid_distances.append(float(left_geometry.centroid.distance(right_geometry.centroid)) / 1000.0) + return { + "min_edge_distance_m": min(edge_distances) if edge_distances else None, + "min_centroid_distance_km": min(centroid_distances) if centroid_distances else None, + "max_centroid_distance_km": max(centroid_distances) if centroid_distances else None, + } + + +def ordered_distinct(values: Iterable[Any]) -> list[str]: + seen = set() + output = [] + for value in values: + text = clean_text(value) + if text is None or text in seen: + continue + seen.add(text) + output.append(text) + return output + + +def read_csv_with_fallback(path: Path) -> pd.DataFrame: + last_error: Exception | None = None + for encoding in ["utf-8-sig", "cp1252", "latin1"]: + try: + return pd.read_csv(path, dtype=str, encoding=encoding) + except UnicodeDecodeError as exc: + last_error = exc + raise SystemExit(f"Could not read {path} with UTF-8, cp1252, or latin1: {last_error}") + + +def load_multipart_notes(path: Path) -> dict[int, dict[str, Any]]: + if not path.exists(): + raise SystemExit(f"Multipart analysis CSV not found: {path}") + dataframe = read_csv_with_fallback(path) + notes: dict[int, dict[str, Any]] = {} + for row in dataframe.to_dict("records"): + village_id = clean_int(row.get("pc11_village_id")) + if village_id is None: + continue + notes[village_id] = row + return notes + + +def list_geojson_files( + admin_input_dir: Path, + *, + state: str | None = None, + district: str | None = None, + limit_files: int | None = None, +) -> list[Path]: + state_target = normalize_slug(state) if state else None + district_target = normalize_slug(district) if district else None + files: list[Path] = [] + for state_dir in sorted(path for path in admin_input_dir.iterdir() if path.is_dir()): + if state_target and normalize_slug(state_dir.name) != state_target: + continue + for path in sorted(state_dir.glob("*.geojson")): + if district_target and normalize_slug(path.stem) != district_target: + continue + files.append(path) + if limit_files is not None and len(files) >= limit_files: + return files + return files + + +def open_work_connection(sqlite_path: Path, *, overwrite: bool) -> sqlite3.Connection: + if sqlite_path.exists(): + if not overwrite: + raise FileExistsError(f"{sqlite_path} already exists. Pass --overwrite to rebuild it.") + sqlite_path.unlink() + sqlite_path.parent.mkdir(parents=True, exist_ok=True) + connection = sqlite3.connect(sqlite_path) + connection.row_factory = sqlite3.Row + connection.execute("PRAGMA journal_mode=WAL") + connection.execute("PRAGMA synchronous=NORMAL") + connection.execute("PRAGMA temp_store=MEMORY") + return connection + + +def create_work_schema(connection: sqlite3.Connection) -> None: + connection.executescript( + """ + CREATE TABLE source_rows ( + row_id INTEGER PRIMARY KEY, + source_file TEXT NOT NULL, + source_state_slug TEXT NOT NULL, + source_district_slug TEXT NOT NULL, + source_row_number INTEGER NOT NULL, + source_group_key TEXT NOT NULL, + identity_hash TEXT NOT NULL, + geometry_hash TEXT, + pc11_village_id INTEGER, + state_name TEXT, + district_name TEXT, + TEHSIL TEXT, + NAME TEXT, + pc11_state_id INTEGER, + pc11_district_id INTEGER, + pc11_subdistrict_id INTEGER + ); + + CREATE INDEX source_rows_identity_hash_idx ON source_rows(identity_hash); + CREATE INDEX source_rows_source_group_idx ON source_rows(source_group_key); + CREATE INDEX source_rows_pc11_village_id_idx ON source_rows(pc11_village_id); + + CREATE TABLE identities ( + identity_hash TEXT PRIMARY KEY, + source_group_key TEXT NOT NULL, + geometry_hash TEXT, + geom_wkb BLOB, + source_row_count INTEGER NOT NULL DEFAULT 0, + first_source_file TEXT NOT NULL, + first_source_row_number INTEGER NOT NULL, + state_name TEXT, + district_name TEXT, + TEHSIL TEXT, + pc11_village_id INTEGER, + NAME TEXT, + pc11_state_id INTEGER, + pc11_district_id INTEGER, + pc11_subdistrict_id INTEGER + ); + + CREATE INDEX identities_group_idx ON identities(source_group_key); + CREATE INDEX identities_pc11_village_id_idx ON identities(pc11_village_id); + + CREATE TABLE source_file_audit ( + source_file TEXT PRIMARY KEY, + source_state_slug TEXT NOT NULL, + source_district_slug TEXT NOT NULL, + read_status TEXT NOT NULL, + read_error TEXT, + feature_count INTEGER, + field_count INTEGER, + crs TEXT, + rows_read INTEGER NOT NULL, + valid_village_id_rows INTEGER NOT NULL, + invalid_village_id_rows INTEGER NOT NULL, + elapsed_seconds REAL NOT NULL + ); + + CREATE TABLE candidate_features ( + candidate_feature_id TEXT PRIMARY KEY, + source_group_key TEXT NOT NULL, + feature_part_index INTEGER NOT NULL, + feature_part_count INTEGER NOT NULL, + source_identity_hashes TEXT NOT NULL, + village_shape_hash TEXT, + pc11_village_id INTEGER, + state_name TEXT, + district_name TEXT, + TEHSIL TEXT, + NAME TEXT, + pc11_state_id INTEGER, + pc11_district_id INTEGER, + pc11_subdistrict_id INTEGER, + geom_wkb BLOB NOT NULL, + source_row_count INTEGER NOT NULL, + source_identity_count INTEGER NOT NULL, + source_geometry_count INTEGER NOT NULL, + source_file_count INTEGER NOT NULL, + resolution_status TEXT NOT NULL + ); + + CREATE INDEX candidate_features_pc11_village_id_idx ON candidate_features(pc11_village_id); + CREATE INDEX candidate_features_source_group_idx ON candidate_features(source_group_key); + + CREATE TABLE standard_trace_map ( + cs_feature_id TEXT, + trace_status TEXT NOT NULL, + candidate_feature_id TEXT, + identity_hash TEXT NOT NULL, + standard_reason TEXT + ); + + CREATE INDEX standard_trace_map_feature_idx ON standard_trace_map(cs_feature_id); + CREATE INDEX standard_trace_map_identity_idx ON standard_trace_map(identity_hash); + """ + ) + + +def insert_source_rows(connection: sqlite3.Connection, rows: list[tuple[Any, ...]]) -> None: + if rows: + connection.executemany( + """ + INSERT INTO source_rows ( + source_file, + source_state_slug, + source_district_slug, + source_row_number, + source_group_key, + identity_hash, + geometry_hash, + pc11_village_id, + state_name, + district_name, + TEHSIL, + NAME, + pc11_state_id, + pc11_district_id, + pc11_subdistrict_id + ) VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?) + """, + rows, + ) + + +def insert_identity_rows(connection: sqlite3.Connection, rows: list[tuple[Any, ...]]) -> None: + if rows: + connection.executemany( + """ + INSERT INTO identities ( + identity_hash, + source_group_key, + geometry_hash, + geom_wkb, + source_row_count, + first_source_file, + first_source_row_number, + state_name, + district_name, + TEHSIL, + pc11_village_id, + NAME, + pc11_state_id, + pc11_district_id, + pc11_subdistrict_id + ) VALUES (?, ?, ?, ?, 1, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?) + ON CONFLICT(identity_hash) DO UPDATE SET + source_row_count = source_row_count + 1 + """, + rows, + ) + + +def ingest_source_file(connection: sqlite3.Connection, path: Path) -> dict[str, Any]: + started = time.perf_counter() + source_file = path.relative_to(ROOT_DIR).as_posix() if path.is_relative_to(ROOT_DIR) else path.as_posix() + state_slug = normalize_slug(path.parent.name) + district_slug = normalize_slug(path.stem) + source_batch: list[tuple[Any, ...]] = [] + identity_batch: list[tuple[Any, ...]] = [] + rows_read = 0 + valid_rows = 0 + invalid_rows = 0 + read_status = "ok" + read_error = None + feature_count = None + field_count = None + crs = "" + try: + info = pyogrio.read_info(path) + feature_count = int(info.get("features") or 0) + fields = info.get("fields") + field_count = len(fields) if fields is not None else 0 + crs = str(info.get("crs") or "") + gdf = pyogrio.read_dataframe(path) + for row_number, row in enumerate(gdf.itertuples(index=False), start=1): + rows_read += 1 + row_data = row._asdict() + properties = { + column: normalize_property(column, row_data.get(column)) + for column in ADMIN_COLUMNS + } + geometry = repair_geometry(row_data.get(gdf.geometry.name)) + geometry_hash = stable_shape_hash(geometry) + identity_hash = stable_identity_hash(properties, geometry_hash) + source_group_key = source_group_key_for(properties, identity_hash) + if properties["pc11_village_id"] is None: + invalid_rows += 1 + else: + valid_rows += 1 + source_batch.append( + ( + source_file, + state_slug, + district_slug, + row_number, + source_group_key, + identity_hash, + geometry_hash, + properties["pc11_village_id"], + properties["state_name"], + properties["district_name"], + properties["TEHSIL"], + properties["NAME"], + properties["pc11_state_id"], + properties["pc11_district_id"], + properties["pc11_subdistrict_id"], + ) + ) + if geometry is not None: + identity_batch.append( + ( + identity_hash, + source_group_key, + geometry_hash, + sqlite3.Binary(geometry.wkb), + source_file, + row_number, + properties["state_name"], + properties["district_name"], + properties["TEHSIL"], + properties["pc11_village_id"], + properties["NAME"], + properties["pc11_state_id"], + properties["pc11_district_id"], + properties["pc11_subdistrict_id"], + ) + ) + if len(source_batch) >= SQLITE_BATCH_SIZE: + insert_source_rows(connection, source_batch) + insert_identity_rows(connection, identity_batch) + source_batch.clear() + identity_batch.clear() + insert_source_rows(connection, source_batch) + insert_identity_rows(connection, identity_batch) + connection.commit() + except Exception as exc: + read_status = "error" + read_error = str(exc) + audit = { + "source_file": source_file, + "source_state_slug": state_slug, + "source_district_slug": district_slug, + "read_status": read_status, + "read_error": read_error, + "feature_count": feature_count, + "field_count": field_count, + "crs": crs, + "rows_read": rows_read, + "valid_village_id_rows": valid_rows, + "invalid_village_id_rows": invalid_rows, + "elapsed_seconds": round(time.perf_counter() - started, 6), + } + connection.execute( + """ + INSERT INTO source_file_audit ( + source_file, + source_state_slug, + source_district_slug, + read_status, + read_error, + feature_count, + field_count, + crs, + rows_read, + valid_village_id_rows, + invalid_village_id_rows, + elapsed_seconds + ) VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?) + """, + ( + audit["source_file"], + audit["source_state_slug"], + audit["source_district_slug"], + audit["read_status"], + audit["read_error"], + audit["feature_count"], + audit["field_count"], + audit["crs"], + audit["rows_read"], + audit["valid_village_id_rows"], + audit["invalid_village_id_rows"], + audit["elapsed_seconds"], + ), + ) + connection.commit() + return audit + + +def load_identity_geometry(row: sqlite3.Row) -> Any | None: + if row["geom_wkb"] is None: + return None + return repair_geometry(load_wkb(bytes(row["geom_wkb"]))) + + +def dissolve_geometries(geometries: Sequence[Any]) -> Any | None: + repaired = [repair_geometry(geometry) for geometry in geometries] + repaired = [geometry for geometry in repaired if geometry is not None] + if not repaired: + return None + if len(repaired) == 1: + return repaired[0] + return repair_geometry(unary_union(repaired)) + + +def should_dissolve_geometry_pair(left_geometry: Any, right_geometry: Any) -> bool: + if left_geometry.equals(right_geometry): + return True + if not left_geometry.intersects(right_geometry): + return False + intersection = left_geometry.intersection(right_geometry) + if intersection.is_empty: + return False + return float(intersection.area) > GEOMETRY_OVERLAP_AREA_EPSILON + + +def geometry_components(identity_rows: Sequence[sqlite3.Row]) -> list[list[sqlite3.Row]]: + rows_with_geometry: list[tuple[sqlite3.Row, Any]] = [] + for row in identity_rows: + geometry = load_identity_geometry(row) + if geometry is not None: + rows_with_geometry.append((row, geometry)) + if not rows_with_geometry: + return [] + + parent = list(range(len(rows_with_geometry))) + + def find(index: int) -> int: + while parent[index] != index: + parent[index] = parent[parent[index]] + index = parent[index] + return index + + def union(left: int, right: int) -> None: + left_root = find(left) + right_root = find(right) + if left_root != right_root: + parent[right_root] = left_root + + for left_index, (_, left_geometry) in enumerate(rows_with_geometry): + for right_index in range(left_index + 1, len(rows_with_geometry)): + _, right_geometry = rows_with_geometry[right_index] + if should_dissolve_geometry_pair(left_geometry, right_geometry): + union(left_index, right_index) + + grouped: dict[int, list[sqlite3.Row]] = {} + order: list[int] = [] + for index, (row, _) in enumerate(rows_with_geometry): + root = find(index) + if root not in grouped: + grouped[root] = [] + order.append(root) + grouped[root].append(row) + return [grouped[root] for root in order] + + +def iter_identity_groups(connection: sqlite3.Connection) -> Iterable[tuple[str, list[sqlite3.Row]]]: + current_key = None + current_rows: list[sqlite3.Row] = [] + for row in connection.execute( + """ + SELECT * + FROM identities + ORDER BY source_group_key, first_source_file, first_source_row_number, identity_hash + """ + ): + source_group_key = row["source_group_key"] + if current_key is not None and source_group_key != current_key: + yield current_key, current_rows + current_rows = [] + current_key = source_group_key + current_rows.append(row) + if current_key is not None: + yield current_key, current_rows + + +def choose_property(values: Sequence[Any], *, numeric: bool) -> Any: + distinct = ordered_distinct(values) + if not distinct: + return None + if len(distinct) == 1: + return clean_int(distinct[0]) if numeric else distinct[0] + return clean_int(distinct[0]) if numeric else " | ".join(str(value) for value in distinct[:12]) + + +def candidate_feature_key(source_group_key: str, part_index: int, part_count: int) -> str: + if part_count == 1: + return source_group_key + return f"{source_group_key}:part:{part_index}" + + +def classify_candidate_resolution( + row_count: int, + identity_count: int, + geometry_count: int, + *, + part_count: int, +) -> str: + if part_count > 1: + if row_count == 1 and identity_count == 1 and geometry_count == 1: + return "split_distinct_geometry_part" + if geometry_count == 1: + return "split_collapsed_duplicate_geometry_part" + return "split_dissolved_overlapping_geometry_part" + if row_count == 1 and identity_count == 1 and geometry_count == 1: + return "single_source_row" + if identity_count == 1 and geometry_count == 1: + return "collapsed_exact_duplicate_rows" + if geometry_count == 1: + return "collapsed_attribute_variants_same_geometry" + return "dissolved_overlapping_geometries" + + +def flush_candidate_features(connection: sqlite3.Connection, rows: list[tuple[Any, ...]]) -> None: + if rows: + connection.executemany( + """ + INSERT INTO candidate_features ( + candidate_feature_id, + source_group_key, + feature_part_index, + feature_part_count, + source_identity_hashes, + village_shape_hash, + pc11_village_id, + state_name, + district_name, + TEHSIL, + NAME, + pc11_state_id, + pc11_district_id, + pc11_subdistrict_id, + geom_wkb, + source_row_count, + source_identity_count, + source_geometry_count, + source_file_count, + resolution_status + ) VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?) + """, + rows, + ) + + +def build_candidate_features(connection: sqlite3.Connection) -> dict[str, Any]: + started = time.perf_counter() + insert_rows: list[tuple[Any, ...]] = [] + resolution_counts: Counter[str] = Counter() + source_group_count = 0 + output_count = 0 + skipped_empty = 0 + split_source_groups = 0 + + for source_group_key, identity_rows in iter_identity_groups(connection): + source_group_count += 1 + components = geometry_components(identity_rows) + part_count = len(components) + if part_count == 0: + skipped_empty += 1 + continue + if part_count > 1: + split_source_groups += 1 + for part_index, component_rows in enumerate(components, start=1): + geometries = [load_identity_geometry(row) for row in component_rows] + geometry = dissolve_geometries(geometries) + if geometry is None: + skipped_empty += 1 + continue + source_row_count = sum(int(row["source_row_count"]) for row in component_rows) + identity_count = len(component_rows) + geometry_count = len({row["geometry_hash"] for row in component_rows if row["geometry_hash"]}) + source_file_count = len({row["first_source_file"] for row in component_rows}) + source_identity_hashes = "|".join(row["identity_hash"] for row in component_rows) + properties = { + column: choose_property( + [row[column] for row in component_rows], + numeric=column in ID_COLUMNS, + ) + for column in ADMIN_COLUMNS + } + resolution_status = classify_candidate_resolution( + source_row_count, + identity_count, + geometry_count, + part_count=part_count, + ) + resolution_counts[resolution_status] += 1 + insert_rows.append( + ( + candidate_feature_key(source_group_key, part_index, part_count), + source_group_key, + part_index, + part_count, + source_identity_hashes, + stable_shape_hash(geometry), + properties["pc11_village_id"], + properties["state_name"], + properties["district_name"], + properties["TEHSIL"], + properties["NAME"], + properties["pc11_state_id"], + properties["pc11_district_id"], + properties["pc11_subdistrict_id"], + sqlite3.Binary(geometry.wkb), + source_row_count, + identity_count, + geometry_count, + source_file_count, + resolution_status, + ) + ) + output_count += 1 + if len(insert_rows) >= SQLITE_BATCH_SIZE: + flush_candidate_features(connection, insert_rows) + insert_rows.clear() + flush_candidate_features(connection, insert_rows) + connection.commit() + return { + "candidate_source_groups": source_group_count, + "candidate_features": output_count, + "candidate_skipped_empty": skipped_empty, + "candidate_split_source_groups": split_source_groups, + "candidate_resolution_status_counts": dict(resolution_counts), + "candidate_build_seconds": round(time.perf_counter() - started, 6), + } + + +def candidate_rows_to_gdf(rows: Sequence[sqlite3.Row]) -> Any: + records = [] + for row in rows: + geometry = repair_geometry(load_wkb(bytes(row["geom_wkb"]))) if row["geom_wkb"] is not None else None + records.append( + { + "cs_feature_id": row["candidate_feature_id"], + "source_group_key": row["source_group_key"], + "feature_part_index": row["feature_part_index"], + "feature_part_count": row["feature_part_count"], + "source_identity_hashes": row["source_identity_hashes"], + "village_shape_hash": row["village_shape_hash"], + "state_name": row["state_name"], + "district_name": row["district_name"], + "TEHSIL": row["TEHSIL"], + "pc11_village_id": row["pc11_village_id"], + "NAME": row["NAME"], + "pc11_state_id": row["pc11_state_id"], + "pc11_district_id": row["pc11_district_id"], + "pc11_subdistrict_id": row["pc11_subdistrict_id"], + "source_row_count": row["source_row_count"], + "source_identity_count": row["source_identity_count"], + "source_geometry_count": row["source_geometry_count"], + "source_file_count": row["source_file_count"], + "resolution_status": row["resolution_status"], + "geometry": geometry, + } + ) + return gpd.GeoDataFrame(records, geometry="geometry", crs="EPSG:4326") + + +def standard_group_key(village_id: Any, fallback: str) -> str: + number = clean_int(village_id) + return f"village:{number}" if number is not None else fallback + + +def standard_uid(village_id: Any, *, part_index: int, part_count: int, fallback: str) -> str: + base = standard_group_key(village_id, fallback=fallback) + if clean_int(village_id) is None: + return fallback + return base if part_count == 1 else f"{base}:part:{part_index}" + + +def note_text(note_row: dict[str, Any] | None) -> str: + if not note_row: + return "" + return clean_text(note_row.get("notes")) or "" + + +def has_any(text: str, needles: Sequence[str]) -> bool: + lowered = text.lower() + return any(needle in lowered for needle in needles) + + +def classify_note_action(note: str) -> str | None: + lowered = note.lower() + if "discard outdated" in lowered and "telangana" in lowered: + return "prefer_telangana_duplicate" + explicit_keep = has_any( + lowered, + [ + "keep them both", + "keep them multipart", + "remain as multiparts", + "should remain", + "valid multipart", + "multipart village", + "local cluster", + "probable towns", + "townships", + "different village names", + "no probable fuzzy edge", + "pass a valid multipart", + ], + ) + negative_merge = has_any( + lowered, + [ + "not a good merge", + "should remain", + "keep them both", + "remain as multiparts", + "valid multipart", + ], + ) + explicit_merge = has_any( + lowered, + [ + "simply be merged", + "unary", + "candidate for merging", + "perfect candiates for merging", + "perfect candidates for merging", + "probably just be merged", + "one row with all parts together", + "complete surrounding", + "surrounding another, merge", + " split_distinct_geometry_part, merge", + ], + ) + if explicit_keep: + return "keep_multipart" + if explicit_merge and not negative_merge: + return "merge_parts" + return None + + +def all_geometries_equal(geometries: Sequence[Any]) -> bool: + repaired = [repair_geometry(geometry) for geometry in geometries] + repaired = [geometry for geometry in repaired if geometry is not None] + if len(repaired) < 2: + return False + first = repaired[0] + return all(first.equals(other) for other in repaired[1:]) + + +def choose_canonical_index(group: Any, areas: Sequence[float], *, prefer_state: str | None = None) -> Any: + ranking = group.copy() + ranking["_standard_area_km2"] = list(areas) + if prefer_state: + preferred = ranking[ranking["state_name"].map(normalize_label) == prefer_state] + if not preferred.empty: + ranking = preferred + ranking = ranking.sort_values( + by=["_standard_area_km2", "source_row_count", "source_identity_count"], + ascending=[False, False, False], + kind="mergesort", + ) + return ranking.index[0] + + +def geometry_set_signature(group: Any) -> str: + hashes = [] + for _, row in group.iterrows(): + shape_hash = clean_text(row.get("village_shape_hash")) + if shape_hash is None: + shape_hash = stable_shape_hash(repair_geometry(row.geometry)) + hashes.append(shape_hash or "") + payload = "\n".join(sorted(hashes)).encode("utf-8") + return hashlib.sha256(payload).hexdigest() + + +def is_cross_duplicate_instruction(note: str) -> bool: + lowered = note.lower() + return has_any( + lowered, + [ + "same identical group", + "identical same parts", + "repeated pattern", + "identical feature rows", + "same local urban", + "seem to be 3 identical rows", + ], + ) + + +def decision_for_group( + group: Any, + *, + note_row: dict[str, Any] | None, + dominant_share_threshold: float, + near_distance_m: float, + far_centroid_km: float, + many_part_threshold: int, +) -> dict[str, Any]: + geometries = [repair_geometry(geometry) for geometry in group.geometry.tolist()] + areas = [geodesic_area_km2(geometry) for geometry in geometries] + total_area = sum(areas) + largest_area = max(areas) if areas else 0.0 + largest_share = largest_area / total_area if total_area > 0 else 0.0 + distances = group_distance_metrics(geometries) + names = [normalize_label(value) for value in group["NAME"].tolist()] + names = [value for value in names if value] + name_count = len(set(names)) + states = {normalize_label(value) for value in group["state_name"].tolist() if normalize_label(value)} + districts = {normalize_label(value) for value in group["district_name"].tolist() if normalize_label(value)} + tehsils = {normalize_label(value) for value in group["TEHSIL"].tolist() if normalize_label(value)} + geometry_hashes = set(group["village_shape_hash"].dropna().astype(str)) + note = note_text(note_row) + note_action = classify_note_action(note) + identical_geometries = all_geometries_equal(geometries) or (len(geometry_hashes) == 1 and len(group) > 1) + + pattern_flags = [] + if identical_geometries: + pattern_flags.append("exact_duplicate_geometry") + if largest_share >= dominant_share_threshold: + pattern_flags.append("dominant_largest_part") + if name_count > 1: + pattern_flags.append("multiple_names") + if len(states) > 1 or len(districts) > 1 or len(tehsils) > 1: + pattern_flags.append("admin_variants") + if distances["min_edge_distance_m"] is not None and distances["min_edge_distance_m"] <= near_distance_m: + pattern_flags.append("near_or_touching_parts") + if distances["max_centroid_distance_km"] is not None and distances["max_centroid_distance_km"] >= far_centroid_km: + pattern_flags.append("far_apart_parts") + if len(group) >= many_part_threshold: + pattern_flags.append("many_part_group") + if note: + pattern_flags.append("manual_note") + + if note_action == "prefer_telangana_duplicate": + action = "keep_preferred_duplicate" + reason = "manual_note_prefers_telangana_over_outdated_duplicate" + canonical_index = choose_canonical_index(group, areas, prefer_state="TELANGANA") + elif identical_geometries and len(states) > 1 and "TELANGANA" in states: + action = "keep_preferred_duplicate" + reason = "identical_geometry_state_conflict_prefers_telangana" + canonical_index = choose_canonical_index(group, areas, prefer_state="TELANGANA") + elif identical_geometries: + action = "keep_largest_duplicate" + reason = "identical_geometry_duplicate" + canonical_index = choose_canonical_index(group, areas) + elif note_action == "keep_multipart": + action = "keep_multipart" + reason = "manual_note_keep_multipart" + canonical_index = choose_canonical_index(group, areas) + elif note_action == "merge_parts": + action = "merge_parts" + reason = "manual_note_merge" + canonical_index = choose_canonical_index(group, areas) + elif largest_share >= dominant_share_threshold: + action = "merge_parts" + reason = f"largest_part_share_ge_{dominant_share_threshold:g}" + canonical_index = choose_canonical_index(group, areas) + elif ( + distances["max_centroid_distance_km"] is not None + and distances["max_centroid_distance_km"] >= far_centroid_km + ): + action = "keep_multipart" + reason = f"far_apart_parts_ge_{far_centroid_km:g}km" + canonical_index = choose_canonical_index(group, areas) + elif name_count > 1 and len(states) <= 1 and len(districts) <= 1 and len(tehsils) <= 1: + action = "keep_multipart" + reason = "multiple_names_same_admin_context" + canonical_index = choose_canonical_index(group, areas) + elif len(group) >= many_part_threshold: + action = "keep_multipart" + reason = f"many_part_group_ge_{many_part_threshold}" + canonical_index = choose_canonical_index(group, areas) + elif ( + distances["min_edge_distance_m"] is not None + and distances["min_edge_distance_m"] <= near_distance_m + ): + action = "merge_parts" + reason = f"near_or_touching_parts_le_{near_distance_m:g}m" + canonical_index = choose_canonical_index(group, areas) + else: + action = "keep_multipart" + reason = "conservative_default_keep_multipart" + canonical_index = choose_canonical_index(group, areas) + + return { + "action": action, + "reason": reason, + "canonical_index": canonical_index, + "areas_km2": areas, + "total_area_km2": total_area, + "largest_area_km2": largest_area, + "largest_area_share": largest_share, + "min_edge_distance_m": distances["min_edge_distance_m"], + "min_centroid_distance_km": distances["min_centroid_distance_km"], + "max_centroid_distance_km": distances["max_centroid_distance_km"], + "name_count": name_count, + "state_count": len(states), + "district_count": len(districts), + "tehsil_count": len(tehsils), + "pattern_flags": pattern_flags, + "note": note, + } + + +def apply_forced_decision(decision: dict[str, Any], *, action: str, reason: str, flag: str) -> dict[str, Any]: + updated = dict(decision) + updated["action"] = action + updated["reason"] = reason + flags = list(updated.get("pattern_flags") or []) + if flag not in flags: + flags.append(flag) + updated["pattern_flags"] = flags + return updated + + +def aggregate_count(values: Any) -> int | None: + numbers = pd.to_numeric(values, errors="coerce").dropna() + return int(numbers.sum()) if not numbers.empty else None + + +def prepare_candidate_standard_columns(gdf: Any, *, action: str, reason: str) -> Any: + output = gdf.copy() + source_ids = output["cs_feature_id"].astype(str) + output["_candidate_feature_ids"] = source_ids + output["_source_identity_hashes"] = output.get("source_identity_hashes", "") + output["_trace_status"] = "output" + output["standard_group_key"] = [ + standard_group_key(village_id, fallback=source_id) + for village_id, source_id in zip(output["pc11_village_id"], source_ids) + ] + output["cs_admin_uid"] = [ + standard_uid(village_id, part_index=1, part_count=1, fallback=source_id) + for village_id, source_id in zip(output["pc11_village_id"], source_ids) + ] + output["cs_feature_id"] = output["cs_admin_uid"] + output["core_admin_uid"] = output["cs_admin_uid"] + output["village_id"] = pd.to_numeric(output["pc11_village_id"], errors="coerce").astype("Int64") + output["standard_action"] = action + output["standard_reason"] = reason + output["standard_notes"] = "" + output["standard_pattern_flags"] = "" + output["standard_group_part_count"] = 1 + output["standard_part_area_km2"] = float("nan") + output["standard_part_area_share"] = float("nan") + output["standard_group_area_km2"] = float("nan") + output["standard_min_edge_distance_m"] = float("nan") + output["standard_min_centroid_distance_km"] = float("nan") + output["standard_max_centroid_distance_km"] = float("nan") + output["standard_name_count"] = 1 + return coerce_and_order_output(output) + + +def coerce_and_order_output(gdf: Any) -> Any: + output = gdf.copy() + internal_columns = [ + column for column in ["_candidate_feature_ids", "_source_identity_hashes", "_trace_status"] + if column in output.columns + ] + for column in STANDARD_COLUMNS: + if column not in output.columns: + output[column] = None + for column in ADMIN_INT_COLUMNS: + if column in output.columns: + output[column] = pd.to_numeric(output[column], errors="coerce").astype("Int64") + for column in [ + "standard_part_area_km2", + "standard_part_area_share", + "standard_group_area_km2", + "standard_min_edge_distance_m", + "standard_min_centroid_distance_km", + "standard_max_centroid_distance_km", + ]: + output[column] = pd.to_numeric(output[column], errors="coerce").astype("Float64") + geometry_column = output.geometry.name + return output[[*STANDARD_COLUMNS, *internal_columns, geometry_column]] + + +def public_output_gdf(gdf: Any) -> Any: + geometry_column = gdf.geometry.name + columns = [column for column in STANDARD_COLUMNS if column in gdf.columns] + return gdf[[*columns, geometry_column]] + + +def make_output_row_from_series(row: Any, geometry: Any) -> dict[str, Any]: + record = row.to_dict() + record["geometry"] = geometry + return record + + +def transform_haryana_800372_group( + group: Any, + *, + note_row: dict[str, Any] | None, +) -> tuple[Any, list[dict[str, Any]], list[dict[str, Any]], dict[str, Any]]: + note = note_text(note_row) + village_id = 800372 + target_candidate_feature_id = "village:haryana:kurukshetra:thanesar:800372:part:1" + output_rows: list[dict[str, Any]] = [] + omission_rows: list[dict[str, Any]] = [] + part_rows: list[dict[str, Any]] = [] + areas = [geodesic_area_km2(repair_geometry(geometry)) for geometry in group.geometry.tolist()] + total_area = sum(areas) + distances = group_distance_metrics(group.geometry.tolist()) + canonical_matches = group[group["cs_feature_id"].astype(str) == target_candidate_feature_id] + canonical_index = canonical_matches.index[0] if not canonical_matches.empty else choose_canonical_index(group, areas) + + for position, (index, source_row) in enumerate(group.iterrows()): + area = float(areas[position]) + source_id = str(source_row["cs_feature_id"]) + geometry = repair_geometry(source_row.geometry) + row = source_row.copy() + row["_candidate_feature_ids"] = source_id + row["_source_identity_hashes"] = source_row.get("source_identity_hashes", "") + row["_trace_status"] = "output" + is_canonical = index == canonical_index + if is_canonical: + output_id = standard_uid(village_id, part_index=1, part_count=1, fallback=source_id) + row["pc11_village_id"] = village_id + row["village_id"] = village_id + row["cs_feature_id"] = output_id + row["cs_admin_uid"] = output_id + row["core_admin_uid"] = output_id + row["standard_group_key"] = f"village:{village_id}" + row["feature_part_index"] = 1 + row["feature_part_count"] = 1 + row["resolution_status"] = "standard_haryana_800372_kept_id_part" + row["standard_action"] = "keep_selected_part_reclassify_rest" + row["standard_reason"] = "manual_note_keep_only_kurukshetra_thanesar_part_1_with_village_id" + else: + shape_hash = stable_shape_hash(geometry) or hashlib.sha256(source_id.encode("utf-8")).hexdigest() + output_id = f"invalid:{shape_hash[:24]}" + row["pc11_village_id"] = pd.NA + row["village_id"] = pd.NA + row["cs_feature_id"] = output_id + row["cs_admin_uid"] = output_id + row["core_admin_uid"] = output_id + row["standard_group_key"] = output_id + row["feature_part_index"] = 1 + row["feature_part_count"] = 1 + row["resolution_status"] = "standard_reclassified_invalid_fragment" + row["standard_action"] = "reclassify_fragment_without_village_id" + row["standard_reason"] = "manual_note_remove_800372_from_nonselected_parts" + omission_rows.append( + { + "pc11_village_id": village_id, + "candidate_feature_id": source_id, + "omission_reason": "village_id_removed_but_geometry_retained_as_invalid_fragment", + "kept_candidate_feature_id": target_candidate_feature_id, + "standard_action": "reclassify_fragment_without_village_id", + "state_name": source_row["state_name"], + "district_name": source_row["district_name"], + "TEHSIL": source_row["TEHSIL"], + "NAME": source_row["NAME"], + } + ) + row["village_shape_hash"] = stable_shape_hash(geometry) + row["standard_notes"] = note + row["standard_pattern_flags"] = "manual_note|many_part_group|haryana_800372_reclassification" + row["standard_group_part_count"] = len(group) + row["standard_part_area_km2"] = round(area, 6) + row["standard_part_area_share"] = round(area / total_area, 8) if total_area > 0 else None + row["standard_group_area_km2"] = round(total_area, 6) + row["standard_min_edge_distance_m"] = distances["min_edge_distance_m"] + row["standard_min_centroid_distance_km"] = distances["min_centroid_distance_km"] + row["standard_max_centroid_distance_km"] = distances["max_centroid_distance_km"] + row["standard_name_count"] = len(set(normalize_label(value) for value in group["NAME"] if normalize_label(value))) + output_rows.append(make_output_row_from_series(row, geometry)) + part_rows.append( + { + "pc11_village_id": village_id, + "candidate_feature_id": source_id, + "source_group_key": source_row["source_group_key"], + "source_part_index": source_row["feature_part_index"], + "source_part_count": source_row["feature_part_count"], + "state_name": source_row["state_name"], + "district_name": source_row["district_name"], + "TEHSIL": source_row["TEHSIL"], + "NAME": source_row["NAME"], + "area_km2": round(area, 6), + "area_share": round(area / total_area, 8) if total_area > 0 else None, + "standard_action": row["standard_action"], + "standard_reason": row["standard_reason"], + } + ) + + output_gdf = gpd.GeoDataFrame(output_rows, geometry="geometry", crs=group.crs) + output_gdf = coerce_and_order_output(output_gdf) + decision_row = { + "pc11_village_id": village_id, + "input_rows": len(group), + "output_rows": len(output_gdf), + "omitted_rows": len(omission_rows), + "standard_action": "keep_selected_part_reclassify_rest", + "standard_reason": "manual_note_keep_only_kurukshetra_thanesar_part_1_with_village_id", + "pattern_flags": "manual_note|many_part_group|haryana_800372_reclassification", + "state_count": len(set(normalize_label(value) for value in group["state_name"] if normalize_label(value))), + "district_count": len(set(normalize_label(value) for value in group["district_name"] if normalize_label(value))), + "tehsil_count": len(set(normalize_label(value) for value in group["TEHSIL"] if normalize_label(value))), + "name_count": len(set(normalize_label(value) for value in group["NAME"] if normalize_label(value))), + "total_area_km2": round(total_area, 6), + "largest_area_km2": round(max(areas), 6) if areas else 0, + "largest_area_share": round(max(areas) / total_area, 8) if total_area > 0 and areas else 0, + "min_edge_distance_m": distances["min_edge_distance_m"], + "min_centroid_distance_km": distances["min_centroid_distance_km"], + "max_centroid_distance_km": distances["max_centroid_distance_km"], + "note": note, + "candidate_feature_ids": "|".join(str(value) for value in group["cs_feature_id"].tolist()[:20]) + + (f"|...(+{len(group) - 20})" if len(group) > 20 else ""), + } + return output_gdf, omission_rows, part_rows, decision_row + + +def omitted_cross_duplicate_group_reports( + group: Any, + *, + note_row: dict[str, Any] | None, + kept_village_id: int | None, + signature: str, + args: argparse.Namespace, +) -> tuple[list[dict[str, Any]], list[dict[str, Any]], dict[str, Any]]: + decision = decision_for_group( + group, + note_row=note_row, + dominant_share_threshold=args.dominant_share_threshold, + near_distance_m=args.near_distance_m, + far_centroid_km=args.far_centroid_km, + many_part_threshold=args.many_part_threshold, + ) + village_id = clean_int(group["pc11_village_id"].iloc[0]) + omission_rows = [] + part_rows = [] + areas = decision["areas_km2"] + total_area = float(decision["total_area_km2"] or 0.0) + for position, (_, row) in enumerate(group.iterrows()): + area = float(areas[position]) + input_id = str(row["cs_feature_id"]) + omission_rows.append( + { + "pc11_village_id": village_id, + "candidate_feature_id": input_id, + "omission_reason": "cross_village_duplicate_geometry_set", + "kept_pc11_village_id": kept_village_id, + "geometry_set_signature": signature, + "standard_action": "omit_cross_duplicate_group", + "state_name": row["state_name"], + "district_name": row["district_name"], + "TEHSIL": row["TEHSIL"], + "NAME": row["NAME"], + } + ) + part_rows.append( + { + "pc11_village_id": village_id, + "candidate_feature_id": input_id, + "source_group_key": row["source_group_key"], + "source_part_index": row["feature_part_index"], + "source_part_count": row["feature_part_count"], + "state_name": row["state_name"], + "district_name": row["district_name"], + "TEHSIL": row["TEHSIL"], + "NAME": row["NAME"], + "area_km2": round(area, 6), + "area_share": round(area / total_area, 8) if total_area > 0 else None, + "standard_action": "omit_cross_duplicate_group", + "standard_reason": "cross_village_duplicate_geometry_set", + } + ) + flags = list(decision.get("pattern_flags") or []) + for flag in ["cross_village_duplicate_geometry_set", "omitted_duplicate_group"]: + if flag not in flags: + flags.append(flag) + decision_row = { + "pc11_village_id": village_id, + "input_rows": len(group), + "output_rows": 0, + "omitted_rows": len(group), + "standard_action": "omit_cross_duplicate_group", + "standard_reason": "cross_village_duplicate_geometry_set", + "pattern_flags": "|".join(flags), + "state_count": decision["state_count"], + "district_count": decision["district_count"], + "tehsil_count": decision["tehsil_count"], + "name_count": decision["name_count"], + "total_area_km2": round(float(decision["total_area_km2"] or 0.0), 6), + "largest_area_km2": round(float(decision["largest_area_km2"] or 0.0), 6), + "largest_area_share": round(float(decision["largest_area_share"] or 0.0), 8), + "min_edge_distance_m": decision["min_edge_distance_m"], + "min_centroid_distance_km": decision["min_centroid_distance_km"], + "max_centroid_distance_km": decision["max_centroid_distance_km"], + "note": decision["note"], + "geometry_set_signature": signature, + "kept_pc11_village_id": kept_village_id, + "candidate_feature_ids": "|".join(str(value) for value in group["cs_feature_id"].tolist()[:20]) + + (f"|...(+{len(group) - 20})" if len(group) > 20 else ""), + } + return omission_rows, part_rows, decision_row + + +def transform_multipart_group( + group: Any, + *, + note_row: dict[str, Any] | None, + args: argparse.Namespace, + forced_action: str | None = None, + forced_reason: str | None = None, + forced_flag: str | None = None, +) -> tuple[Any, list[dict[str, Any]], list[dict[str, Any]], dict[str, Any], list[dict[str, Any]]]: + if clean_int(group["pc11_village_id"].iloc[0]) == 800372 or "rest should keep other values" in note_text(note_row).lower(): + output_gdf, omissions, parts, decision = transform_haryana_800372_group(group, note_row=note_row) + return output_gdf, omissions, parts, decision, [] + + decision = decision_for_group( + group, + note_row=note_row, + dominant_share_threshold=args.dominant_share_threshold, + near_distance_m=args.near_distance_m, + far_centroid_km=args.far_centroid_km, + many_part_threshold=args.many_part_threshold, + ) + if forced_action and forced_reason and forced_flag: + decision = apply_forced_decision( + decision, + action=forced_action, + reason=forced_reason, + flag=forced_flag, + ) + source_ids = [str(value) for value in group["cs_feature_id"].tolist()] + source_identity_hashes = [ + value + for value in group.get("source_identity_hashes", pd.Series(dtype=str)).dropna().astype(str).tolist() + if value + ] + village_id = clean_int(group["pc11_village_id"].iloc[0]) + fallback_group = str(group["source_group_key"].iloc[0]) + group_key = standard_group_key(village_id, fallback=fallback_group) + output_rows: list[dict[str, Any]] = [] + omission_rows: list[dict[str, Any]] = [] + part_rows: list[dict[str, Any]] = [] + areas = decision["areas_km2"] + total_area = float(decision["total_area_km2"] or 0.0) + + for position, (index, row) in enumerate(group.iterrows()): + area = float(areas[position]) + part_rows.append( + { + "pc11_village_id": village_id, + "candidate_feature_id": row["cs_feature_id"], + "source_group_key": row["source_group_key"], + "source_part_index": row["feature_part_index"], + "source_part_count": row["feature_part_count"], + "state_name": row["state_name"], + "district_name": row["district_name"], + "TEHSIL": row["TEHSIL"], + "NAME": row["NAME"], + "area_km2": round(area, 6), + "area_share": round(area / total_area, 8) if total_area > 0 else None, + "standard_action": decision["action"], + "standard_reason": decision["reason"], + } + ) + + canonical = group.loc[decision["canonical_index"]].copy() + canonical_geometry = repair_geometry(canonical.geometry) + pattern_flags = "|".join(decision["pattern_flags"]) + + if decision["action"] == "merge_parts": + geometries = [repair_geometry(geometry) for geometry in group.geometry.tolist()] + geometries = [geometry for geometry in geometries if geometry is not None] + merged_geometry = repair_geometry(unary_union(geometries)) if geometries else canonical_geometry + output_id = standard_uid(village_id, part_index=1, part_count=1, fallback=str(canonical["cs_feature_id"])) + canonical["cs_feature_id"] = output_id + canonical["_candidate_feature_ids"] = "|".join(source_ids) + canonical["_source_identity_hashes"] = "|".join(source_identity_hashes) + canonical["_trace_status"] = "output" + canonical["cs_admin_uid"] = output_id + canonical["core_admin_uid"] = output_id + canonical["standard_group_key"] = group_key + canonical["feature_part_index"] = 1 + canonical["feature_part_count"] = 1 + canonical["village_shape_hash"] = stable_shape_hash(merged_geometry) + canonical["village_id"] = village_id + canonical["source_row_count"] = aggregate_count(group["source_row_count"]) + canonical["source_identity_count"] = aggregate_count(group["source_identity_count"]) + canonical["source_geometry_count"] = len(set(group["village_shape_hash"].dropna().astype(str))) + canonical["source_file_count"] = aggregate_count(group["source_file_count"]) + canonical["resolution_status"] = "standard_merged_parts" + canonical["standard_action"] = decision["action"] + canonical["standard_reason"] = decision["reason"] + canonical["standard_notes"] = decision["note"] + canonical["standard_pattern_flags"] = pattern_flags + canonical["standard_group_part_count"] = len(group) + canonical["standard_part_area_km2"] = round(total_area, 6) + canonical["standard_part_area_share"] = 1.0 + canonical["standard_group_area_km2"] = round(total_area, 6) + canonical["standard_min_edge_distance_m"] = decision["min_edge_distance_m"] + canonical["standard_min_centroid_distance_km"] = decision["min_centroid_distance_km"] + canonical["standard_max_centroid_distance_km"] = decision["max_centroid_distance_km"] + canonical["standard_name_count"] = decision["name_count"] + output_rows.append(make_output_row_from_series(canonical, merged_geometry)) + elif decision["action"] in {"keep_preferred_duplicate", "keep_largest_duplicate"}: + output_id = standard_uid(village_id, part_index=1, part_count=1, fallback=str(canonical["cs_feature_id"])) + source_id = str(canonical["cs_feature_id"]) + canonical["cs_feature_id"] = output_id + canonical["_candidate_feature_ids"] = source_id + canonical["_source_identity_hashes"] = canonical.get("source_identity_hashes", "") + canonical["_trace_status"] = "output" + canonical["cs_admin_uid"] = output_id + canonical["core_admin_uid"] = output_id + canonical["standard_group_key"] = group_key + canonical["feature_part_index"] = 1 + canonical["feature_part_count"] = 1 + canonical["village_shape_hash"] = stable_shape_hash(canonical_geometry) + canonical["village_id"] = village_id + canonical["resolution_status"] = "standard_collapsed_duplicate" + canonical["standard_action"] = decision["action"] + canonical["standard_reason"] = decision["reason"] + canonical["standard_notes"] = decision["note"] + canonical["standard_pattern_flags"] = pattern_flags + canonical["standard_group_part_count"] = len(group) + canonical["standard_part_area_km2"] = round(max(areas), 6) if areas else None + canonical["standard_part_area_share"] = 1.0 + canonical["standard_group_area_km2"] = round(total_area, 6) + canonical["standard_min_edge_distance_m"] = decision["min_edge_distance_m"] + canonical["standard_min_centroid_distance_km"] = decision["min_centroid_distance_km"] + canonical["standard_max_centroid_distance_km"] = decision["max_centroid_distance_km"] + canonical["standard_name_count"] = decision["name_count"] + output_rows.append(make_output_row_from_series(canonical, canonical_geometry)) + for _, row in group.drop(index=decision["canonical_index"]).iterrows(): + omission_rows.append( + { + "pc11_village_id": village_id, + "candidate_feature_id": row["cs_feature_id"], + "omission_reason": decision["reason"], + "kept_candidate_feature_id": source_id, + "standard_action": decision["action"], + "state_name": row["state_name"], + "district_name": row["district_name"], + "TEHSIL": row["TEHSIL"], + "NAME": row["NAME"], + } + ) + else: + ranked = group.copy() + ranked["_standard_area_km2"] = areas + ranked = ranked.sort_values( + by=["_standard_area_km2", "cs_feature_id"], + ascending=[False, True], + kind="mergesort", + ) + part_count = len(ranked) + for part_index, (_, row) in enumerate(ranked.iterrows(), start=1): + geometry = repair_geometry(row.geometry) + source_id = str(row["cs_feature_id"]) + output_id = standard_uid(village_id, part_index=part_index, part_count=part_count, fallback=source_id) + row = row.copy() + row["cs_feature_id"] = output_id + row["_candidate_feature_ids"] = source_id + row["_source_identity_hashes"] = row.get("source_identity_hashes", "") + row["_trace_status"] = "output" + row["cs_admin_uid"] = output_id + row["core_admin_uid"] = output_id + row["standard_group_key"] = group_key + row["feature_part_index"] = part_index + row["feature_part_count"] = part_count + row["village_shape_hash"] = stable_shape_hash(geometry) + row["village_id"] = village_id + row["resolution_status"] = "standard_kept_multipart" + row["standard_action"] = "keep_multipart" + row["standard_reason"] = decision["reason"] + row["standard_notes"] = decision["note"] + row["standard_pattern_flags"] = pattern_flags + row["standard_group_part_count"] = part_count + area = float(row["_standard_area_km2"]) + row["standard_part_area_km2"] = round(area, 6) + row["standard_part_area_share"] = round(area / total_area, 8) if total_area > 0 else None + row["standard_group_area_km2"] = round(total_area, 6) + row["standard_min_edge_distance_m"] = decision["min_edge_distance_m"] + row["standard_min_centroid_distance_km"] = decision["min_centroid_distance_km"] + row["standard_max_centroid_distance_km"] = decision["max_centroid_distance_km"] + row["standard_name_count"] = decision["name_count"] + output_rows.append(make_output_row_from_series(row, geometry)) + + output_gdf = gpd.GeoDataFrame(output_rows, geometry="geometry", crs=group.crs) + output_gdf = coerce_and_order_output(output_gdf) + decision_row = { + "pc11_village_id": village_id, + "input_rows": len(group), + "output_rows": len(output_gdf), + "omitted_rows": len(omission_rows), + "standard_action": decision["action"], + "standard_reason": decision["reason"], + "pattern_flags": pattern_flags, + "state_count": decision["state_count"], + "district_count": decision["district_count"], + "tehsil_count": decision["tehsil_count"], + "name_count": decision["name_count"], + "total_area_km2": round(float(decision["total_area_km2"] or 0.0), 6), + "largest_area_km2": round(float(decision["largest_area_km2"] or 0.0), 6), + "largest_area_share": round(float(decision["largest_area_share"] or 0.0), 8), + "min_edge_distance_m": decision["min_edge_distance_m"], + "min_centroid_distance_km": decision["min_centroid_distance_km"], + "max_centroid_distance_km": decision["max_centroid_distance_km"], + "note": decision["note"], + "candidate_feature_ids": "|".join(source_ids[:20]) + (f"|...(+{len(source_ids) - 20})" if len(source_ids) > 20 else ""), + } + return output_gdf, omission_rows, part_rows, decision_row, [] + + +class ChunkedGeoJSONWriter: + """Stream GeoDataFrame chunks into a single GeoJSON FeatureCollection.""" + + def __init__(self, path: Path): + self.path = path + self.handle = None + self.first_feature = True + self.rows_written = 0 + + def __enter__(self): + self.path.parent.mkdir(parents=True, exist_ok=True) + self.handle = self.path.open("w", encoding="utf-8") + self.handle.write('{"type":"FeatureCollection","features":[\n') + return self + + def write(self, gdf: Any) -> int: + if gdf.empty: + return 0 + written = 0 + for feature in gdf.iterfeatures(na="null", drop_id=True): + if not self.first_feature: + self.handle.write(",\n") + self.handle.write(json_dumps(feature)) + self.first_feature = False + written += 1 + self.rows_written += written + return written + + def __exit__(self, exc_type, exc, traceback) -> None: + if self.handle is not None: + self.handle.write("\n]}\n") + self.handle.close() + + +class OutputWriters: + def __init__(self, args: argparse.Namespace, *, trace_connection: sqlite3.Connection): + self.args = args + self.trace_connection = trace_connection + self.first_gpkg_chunk = True + self.first_csv_chunk = True + self.rows_written = 0 + self.geojson_writer = None + + def __enter__(self): + if not self.args.skip_geojson: + self.geojson_writer = ChunkedGeoJSONWriter(repo_path(self.args.output_geojson)).__enter__() + return self + + def write(self, gdf: Any) -> None: + if gdf.empty: + return + gdf = coerce_and_order_output(gdf) + self.record_trace(gdf) + public_gdf = public_output_gdf(gdf) + output_gpkg = repo_path(self.args.output_gpkg) + pyogrio.write_dataframe( + public_gdf, + output_gpkg, + layer=self.args.output_layer, + driver="GPKG", + append=not self.first_gpkg_chunk, + promote_to_multi=True, + ) + self.first_gpkg_chunk = False + if self.geojson_writer is not None: + self.geojson_writer.write(public_gdf) + if not self.args.skip_csv: + csv_path = repo_path(self.args.output_csv) + csv_path.parent.mkdir(parents=True, exist_ok=True) + frame = pd.DataFrame(public_gdf.drop(columns=[public_gdf.geometry.name])) + frame.to_csv( + csv_path, + mode="w" if self.first_csv_chunk else "a", + header=self.first_csv_chunk, + index=False, + ) + self.first_csv_chunk = False + self.rows_written += len(public_gdf) + + def record_trace(self, gdf: Any) -> None: + rows = [] + for _, row in gdf.iterrows(): + identity_hashes = [ + value + for value in str(row.get("_source_identity_hashes") or "").split("|") + if value + ] + candidate_ids = clean_text(row.get("_candidate_feature_ids")) + trace_status = clean_text(row.get("_trace_status")) or "output" + reason = clean_text(row.get("standard_reason")) + for identity_hash in identity_hashes: + rows.append( + ( + row["cs_feature_id"], + trace_status, + candidate_ids, + identity_hash, + reason, + ) + ) + if rows: + self.trace_connection.executemany( + """ + INSERT INTO standard_trace_map ( + cs_feature_id, + trace_status, + candidate_feature_id, + identity_hash, + standard_reason + ) VALUES (?, ?, ?, ?, ?) + """, + rows, + ) + self.trace_connection.commit() + + def __exit__(self, exc_type, exc, traceback) -> None: + if self.geojson_writer is not None: + self.geojson_writer.__exit__(exc_type, exc, traceback) + + +def validate_output_path(path: Path, *, overwrite: bool, label: str) -> None: + if path.exists(): + if not overwrite: + raise SystemExit(f"{label} already exists: {path}. Pass --overwrite to rebuild it.") + path.unlink() + path.parent.mkdir(parents=True, exist_ok=True) + + +def prepare_outputs(args: argparse.Namespace) -> None: + validate_output_path(repo_path(args.output_gpkg), overwrite=args.overwrite, label="GeoPackage") + if not args.skip_geojson: + validate_output_path(repo_path(args.output_geojson), overwrite=args.overwrite, label="GeoJSON") + if not args.skip_csv: + validate_output_path(repo_path(args.output_csv), overwrite=args.overwrite, label="CSV") + + +def ensure_gpkg_indexes(gpkg_path: Path, layer: str) -> None: + if not gpkg_path.exists(): + return + index_columns = [ + ["cs_feature_id"], + ["cs_admin_uid"], + ["core_admin_uid"], + ["pc11_village_id"], + ["village_id"], + ["state_name", "district_name", "TEHSIL"], + ["pc11_state_id", "pc11_district_id", "pc11_subdistrict_id"], + ] + with sqlite3.connect(gpkg_path) as connection: + for columns in index_columns: + suffix = "_".join(column.lower() for column in columns) + index_name = f"idx_{layer}_{suffix}" + columns_sql = ", ".join(quote_identifier(column) for column in columns) + connection.execute( + f"CREATE INDEX IF NOT EXISTS {quote_identifier(index_name)} " + f"ON {quote_identifier(layer)} ({columns_sql})" + ) + connection.commit() + + +def write_gpkg_metadata(gpkg_path: Path, table_name: str, metadata: dict[str, Any], summary: dict[str, Any]) -> None: + with sqlite3.connect(gpkg_path) as connection: + connection.execute( + f""" + CREATE TABLE IF NOT EXISTS {quote_identifier(table_name)} ( + key TEXT PRIMARY KEY, + value TEXT NOT NULL + ) + """ + ) + connection.executemany( + f""" + INSERT INTO {quote_identifier(table_name)} (key, value) + VALUES (?, ?) + ON CONFLICT(key) DO UPDATE SET value = excluded.value + """, + [ + ("metadata", json_dumps(metadata)), + ("summary", json_dumps(summary)), + ], + ) + connection.commit() + + +def insert_review_ids(connection: sqlite3.Connection, review_ids: set[int]) -> None: + connection.execute("DROP TABLE IF EXISTS standard_review_village_ids") + connection.execute("CREATE TEMP TABLE standard_review_village_ids (village_id INTEGER PRIMARY KEY)") + connection.executemany( + "INSERT INTO standard_review_village_ids (village_id) VALUES (?)", + [(value,) for value in sorted(review_ids)], + ) + connection.commit() + + +def collect_review_village_ids(connection: sqlite3.Connection, notes: dict[int, dict[str, Any]]) -> set[int]: + review_ids = set(notes) + for row in connection.execute( + """ + SELECT pc11_village_id + FROM candidate_features + WHERE pc11_village_id IS NOT NULL + GROUP BY pc11_village_id + HAVING COUNT(*) > 1 + """ + ): + village_id = clean_int(row["pc11_village_id"]) + if village_id is not None: + review_ids.add(village_id) + return review_ids + + +def iter_singleton_candidate_chunks( + connection: sqlite3.Connection, + *, + chunk_size: int, +) -> Iterable[Any]: + offset = 0 + while True: + rows = connection.execute( + """ + SELECT * + FROM candidate_features + WHERE pc11_village_id IS NULL + OR pc11_village_id NOT IN ( + SELECT village_id FROM standard_review_village_ids + ) + ORDER BY + pc11_state_id, + pc11_district_id, + pc11_subdistrict_id, + pc11_village_id, + candidate_feature_id + LIMIT ? OFFSET ? + """, + (chunk_size, offset), + ).fetchall() + if not rows: + break + yield candidate_rows_to_gdf(rows) + offset += len(rows) + + +def load_review_candidate_groups( + connection: sqlite3.Connection, + notes: dict[int, dict[str, Any]], +) -> list[tuple[int | None, Any, dict[str, Any] | None, str]]: + rows = connection.execute( + """ + SELECT * + FROM candidate_features + WHERE pc11_village_id IN ( + SELECT village_id FROM standard_review_village_ids + ) + ORDER BY pc11_village_id, candidate_feature_id + """ + ).fetchall() + gdf = ( + candidate_rows_to_gdf(rows) + if rows + else gpd.GeoDataFrame({"geometry": []}, geometry="geometry", crs="EPSG:4326") + ) + groups = [] + if gdf.empty: + return groups + gdf["pc11_village_id"] = pd.to_numeric(gdf["pc11_village_id"], errors="coerce").astype("Int64") + for village_id, group in gdf.groupby("pc11_village_id", dropna=False, sort=True): + clean_village_id = clean_int(village_id) + note_row = notes.get(clean_village_id) if clean_village_id is not None else None + groups.append((clean_village_id, group.copy(), note_row, geometry_set_signature(group))) + return groups + + +def duplicate_signature_plan( + groups: list[tuple[int | None, Any, dict[str, Any] | None, str]] +) -> tuple[dict[str, int | None], set[str], set[str]]: + by_signature: dict[str, list[tuple[int | None, Any, dict[str, Any] | None, str]]] = {} + for item in groups: + by_signature.setdefault(item[3], []).append(item) + keep_by_signature: dict[str, int | None] = {} + duplicate_signatures: set[str] = set() + force_merge_signatures: set[str] = set() + for signature, items in by_signature.items(): + if len(items) <= 1: + continue + notes = [note_text(item[2]) for item in items] + if not any(is_cross_duplicate_instruction(note) for note in notes): + continue + duplicate_signatures.add(signature) + sorted_items = sorted(items, key=lambda item: (item[0] is None, item[0] or 0)) + keep_by_signature[signature] = sorted_items[0][0] + if any( + has_any( + note.lower(), + [ + "keep one row with all parts together", + "identical feature rows", + "seem to be 3 identical rows", + ], + ) + for note in notes + ): + force_merge_signatures.add(signature) + return keep_by_signature, duplicate_signatures, force_merge_signatures + + +def record_omitted_candidate_trace( + connection: sqlite3.Connection, + omission_rows: Sequence[dict[str, Any]], +) -> None: + rows = [] + for omission in omission_rows: + if omission.get("standard_action") == "reclassify_fragment_without_village_id": + continue + candidate_id = clean_text(omission.get("candidate_feature_id")) + if candidate_id is None: + continue + candidate = connection.execute( + """ + SELECT source_identity_hashes + FROM candidate_features + WHERE candidate_feature_id = ? + """, + (candidate_id,), + ).fetchone() + if candidate is None: + continue + for identity_hash in str(candidate["source_identity_hashes"] or "").split("|"): + if identity_hash: + rows.append( + ( + None, + "omitted", + candidate_id, + identity_hash, + clean_text(omission.get("omission_reason")), + ) + ) + if rows: + connection.executemany( + """ + INSERT INTO standard_trace_map ( + cs_feature_id, + trace_status, + candidate_feature_id, + identity_hash, + standard_reason + ) VALUES (?, ?, ?, ?, ?) + """, + rows, + ) + connection.commit() + + +def write_source_trace_table(output_gpkg: Path, work_connection: sqlite3.Connection) -> int: + trace_rows = work_connection.execute( + """ + SELECT + tm.cs_feature_id, + tm.trace_status, + tm.candidate_feature_id, + tm.standard_reason, + sr.source_file, + sr.source_row_number, + sr.source_state_slug, + sr.source_district_slug, + sr.identity_hash, + sr.geometry_hash, + sr.pc11_village_id, + sr.state_name, + sr.district_name, + sr.TEHSIL, + sr.NAME, + sr.pc11_state_id, + sr.pc11_district_id, + sr.pc11_subdistrict_id + FROM standard_trace_map tm + JOIN source_rows sr ON sr.identity_hash = tm.identity_hash + ORDER BY tm.cs_feature_id, tm.trace_status, sr.source_file, sr.source_row_number + """ + ).fetchall() + with sqlite3.connect(output_gpkg) as output_connection: + output_connection.execute("DROP TABLE IF EXISTS cs_admin_standard_source_trace") + output_connection.execute( + """ + CREATE TABLE cs_admin_standard_source_trace ( + cs_feature_id TEXT, + trace_status TEXT NOT NULL, + candidate_feature_id TEXT, + standard_reason TEXT, + source_file TEXT NOT NULL, + source_row_number INTEGER NOT NULL, + source_state_slug TEXT NOT NULL, + source_district_slug TEXT NOT NULL, + identity_hash TEXT NOT NULL, + geometry_hash TEXT, + pc11_village_id INTEGER, + state_name TEXT, + district_name TEXT, + TEHSIL TEXT, + NAME TEXT, + pc11_state_id INTEGER, + pc11_district_id INTEGER, + pc11_subdistrict_id INTEGER + ) + """ + ) + output_connection.executemany( + """ + INSERT INTO cs_admin_standard_source_trace ( + cs_feature_id, + trace_status, + candidate_feature_id, + standard_reason, + source_file, + source_row_number, + source_state_slug, + source_district_slug, + identity_hash, + geometry_hash, + pc11_village_id, + state_name, + district_name, + TEHSIL, + NAME, + pc11_state_id, + pc11_district_id, + pc11_subdistrict_id + ) VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?) + """, + [tuple(row) for row in trace_rows], + ) + output_connection.execute( + "CREATE INDEX IF NOT EXISTS idx_cs_admin_standard_source_trace_feature " + "ON cs_admin_standard_source_trace (cs_feature_id)" + ) + output_connection.execute( + "CREATE INDEX IF NOT EXISTS idx_cs_admin_standard_source_trace_source " + "ON cs_admin_standard_source_trace (source_file, source_row_number)" + ) + output_connection.commit() + return len(trace_rows) + + +def write_source_file_audit_report(reports_dir: Path, connection: sqlite3.Connection) -> None: + dataframe = pd.read_sql_query( + """ + SELECT * + FROM source_file_audit + ORDER BY source_file + """, + connection, + ) + dataframe.to_csv(reports_dir / "source_file_read_audit.csv", index=False) + + +def max_rss_mb() -> float: + # Linux reports ru_maxrss in KiB. + return resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1024.0 + + +def write_reports( + reports_dir: Path, + *, + metadata: dict[str, Any], + summary: dict[str, Any], + decision_rows: list[dict[str, Any]], + omission_rows: list[dict[str, Any]], + part_rows: list[dict[str, Any]], +) -> None: + reports_dir.mkdir(parents=True, exist_ok=True) + (reports_dir / "build_summary.json").write_text( + json_dumps({"metadata": metadata, "summary": summary}, indent=True), + encoding="utf-8", + ) + pd.DataFrame(decision_rows).to_csv(reports_dir / "standard_decision_log.csv", index=False) + pd.DataFrame(omission_rows).to_csv(reports_dir / "standard_omissions.csv", index=False) + pd.DataFrame(part_rows).to_csv(reports_dir / "standard_multipart_parts.csv", index=False) + + pattern_counts = Counter() + for row in decision_rows: + for flag in str(row.get("pattern_flags") or "").split("|"): + if flag: + pattern_counts[flag] += 1 + action_counts = Counter(row["standard_action"] for row in decision_rows) + pd.DataFrame( + [{"pattern_flag": key, "groups": value} for key, value in sorted(pattern_counts.items())] + ).to_csv(reports_dir / "standard_pattern_summary.csv", index=False) + + action_lines = [ + f"- `{key}`: `{value:,}`" + for key, value in sorted(summary.get("standard_action_counts", {}).items()) + ] or ["- None"] + pattern_lines = [ + f"- `{key}`: `{value:,}`" + for key, value in sorted(pattern_counts.items()) + ] or ["- None"] + report = f"""# CS Admin Standard Build Summary + +Generated by `utilities/scripts/admin_assets/build_cs_admin_boundary_standard.py`. + +## Outputs + +- GeoPackage: `{metadata['output_gpkg']}` +- Layer: `{metadata['output_layer']}` +- GeoJSON: `{metadata['output_geojson']}` +- CSV: `{metadata['output_csv']}` + +## Counts + +- Input features scanned: `{summary['input_features']:,}` +- Raw source identities: `{summary['source_identities']:,}` +- Candidate features after raw duplicate/overlap resolution: `{summary['candidate_features']:,}` +- Village ids held for note/multipart review: `{summary['review_village_ids']:,}` +- Unchanged singleton features written: `{summary['singleton_features_written']:,}` +- Standard output features written: `{summary['output_features']:,}` +- Multipart groups evaluated: `{summary['multipart_groups']:,}` +- Omission/reclassification report rows: `{summary['omitted_features']:,}` +- Raw source trace rows in GeoPackage: `{summary['trace_rows']:,}` +- Peak script RSS: `{summary['max_rss_mb']:.1f} MB` + +## Multipart Actions + +{chr(10).join(action_lines)} + +## Pattern Flags + +{chr(10).join(pattern_lines)} + +## Decision Rules + +- Raw boundary features are read from `{metadata['admin_input_dir']}`. +- Manual notes are read from `{metadata['multipart_analysis_csv']}` with an encoding fallback. +- Exact duplicate geometries collapse to one feature. +- Telangana is preferred when a manual note or identical-geometry state conflict identifies old Andhra/Telangana duplicates. +- Groups with a dominant largest part are dissolved unless manual notes explicitly preserve a valid multipart case. +- Lower-confidence groups remain multipart with part indexes sorted by descending area. +- Raw source trace is stored in `cs_admin_standard_source_trace` inside the GeoPackage. +""" + (reports_dir / "build_summary.md").write_text(report, encoding="utf-8") + + +def build_standard_asset(args: argparse.Namespace) -> dict[str, Any]: + ensure_dependencies() + started = time.perf_counter() + admin_input_dir = repo_path(args.admin_input_dir) + multipart_analysis_csv = repo_path(args.multipart_analysis_csv) + output_gpkg = repo_path(args.output_gpkg) + output_geojson = repo_path(args.output_geojson) + output_csv = repo_path(args.output_csv) + reports_dir = repo_path(args.reports_dir) + + if not admin_input_dir.exists(): + raise SystemExit(f"Admin input directory not found: {admin_input_dir}") + notes = load_multipart_notes(multipart_analysis_csv) + files = list_geojson_files( + admin_input_dir, + state=args.state, + district=args.district, + limit_files=args.limit_files, + ) + if not files: + raise SystemExit(f"No district GeoJSON files found under {admin_input_dir}") + + metadata = { + "admin_input_dir": admin_input_dir.as_posix(), + "multipart_analysis_csv": multipart_analysis_csv.as_posix(), + "output_gpkg": output_gpkg.as_posix(), + "output_layer": args.output_layer, + "output_geojson": None if args.skip_geojson else output_geojson.as_posix(), + "output_csv": None if args.skip_csv else output_csv.as_posix(), + "reports_dir": reports_dir.as_posix(), + "source_file_count": len(files), + "state_filter": args.state, + "district_filter": args.district, + "limit_files": args.limit_files, + "chunk_size": args.chunk_size, + "dominant_share_threshold": args.dominant_share_threshold, + "near_distance_m": args.near_distance_m, + "far_centroid_km": args.far_centroid_km, + "many_part_threshold": args.many_part_threshold, + "work_db": repo_path(args.work_db).as_posix() if args.work_db else None, + } + + if args.validate_only: + validation = { + **metadata, + "multipart_note_rows": len(notes), + "first_source_files": [ + path.relative_to(ROOT_DIR).as_posix() if path.is_relative_to(ROOT_DIR) else path.as_posix() + for path in files[:5] + ], + } + print(json_dumps(validation, indent=True), flush=True) + return validation + + prepare_outputs(args) + reports_dir.mkdir(parents=True, exist_ok=True) + + temp_dir: tempfile.TemporaryDirectory[str] | None = None + if args.work_db: + work_db = repo_path(args.work_db) + elif args.keep_work_db: + work_db = reports_dir / "cs_admin_standard_work.sqlite3" + else: + temp_dir = tempfile.TemporaryDirectory(prefix="cs_admin_standard_") + work_db = Path(temp_dir.name) / "cs_admin_standard_work.sqlite3" + + connection = open_work_connection(work_db, overwrite=True) + try: + create_work_schema(connection) + print(f"[standard] ingesting {len(files):,} raw admin GeoJSON files", flush=True) + audits = [] + for index, path in enumerate(files, start=1): + audit = ingest_source_file(connection, path) + audits.append(audit) + if index == 1 or index % 25 == 0 or index == len(files): + print( + f"[standard] ingested {index:,}/{len(files):,}; " + f"latest={path.parent.name}/{path.name}; " + f"rows={audit['rows_read']:,}; status={audit['read_status']}; " + f"rss={max_rss_mb():.1f} MB", + flush=True, + ) + + print("[standard] resolving raw identities into candidate features", flush=True) + candidate_summary = build_candidate_features(connection) + print( + f"[standard] candidate features: {candidate_summary['candidate_features']:,}; " + f"rss={max_rss_mb():.1f} MB", + flush=True, + ) + + review_ids = collect_review_village_ids(connection, notes) + insert_review_ids(connection, review_ids) + review_groups = load_review_candidate_groups(connection, notes) + keep_by_signature, duplicate_signatures, force_merge_signatures = duplicate_signature_plan(review_groups) + + decision_rows: list[dict[str, Any]] = [] + omission_rows: list[dict[str, Any]] = [] + part_rows: list[dict[str, Any]] = [] + action_counts: Counter[str] = Counter() + singleton_features_written = 0 + transformed_rows = 0 + omitted_duplicate_groups = 0 + print( + f"[standard] writing singleton candidate features; " + f"{len(review_ids):,} village ids held for note/multipart review", + flush=True, + ) + + with OutputWriters(args, trace_connection=connection) as writers: + for chunk_index, singleton in enumerate( + iter_singleton_candidate_chunks(connection, chunk_size=args.chunk_size), + start=1, + ): + singleton = prepare_candidate_standard_columns( + singleton, + action="kept_single", + reason="single_candidate_from_raw_inputs", + ) + writers.write(singleton) + singleton_features_written += len(singleton) + if chunk_index == 1 or singleton_features_written % (args.chunk_size * 2) == 0: + print( + f"[standard] singleton features written={singleton_features_written:,}; " + f"rss={max_rss_mb():.1f} MB", + flush=True, + ) + + print(f"[standard] applying note/multipart rules to {len(review_groups):,} groups", flush=True) + for group_index, (clean_village_id, group, note_row, signature) in enumerate(review_groups, start=1): + forced_action = None + forced_reason = None + forced_flag = None + if signature in duplicate_signatures: + keep_village_id = keep_by_signature.get(signature) + if clean_village_id != keep_village_id: + omissions, parts, decision = omitted_cross_duplicate_group_reports( + group.copy(), + note_row=note_row, + kept_village_id=keep_village_id, + signature=signature, + args=args, + ) + decision_rows.append(decision) + omission_rows.extend(omissions) + part_rows.extend(parts) + record_omitted_candidate_trace(connection, omissions) + action_counts[decision["standard_action"]] += 1 + omitted_duplicate_groups += 1 + if group_index == 1 or group_index % 500 == 0 or group_index == len(review_groups): + print( + f"[standard] review groups {group_index:,}/{len(review_groups):,}; " + f"transformed_rows={transformed_rows:,}; omitted_duplicate_groups={omitted_duplicate_groups:,}; " + f"rss={max_rss_mb():.1f} MB", + flush=True, + ) + continue + if signature in force_merge_signatures: + forced_action = "merge_parts" + forced_reason = "cross_village_duplicate_geometry_set_keep_one_row_with_all_parts" + forced_flag = "cross_village_duplicate_canonical" + + output_gdf, omissions, parts, decision, _ = transform_multipart_group( + group.copy(), + note_row=note_row, + args=args, + forced_action=forced_action, + forced_reason=forced_reason, + forced_flag=forced_flag, + ) + writers.write(output_gdf) + transformed_rows += len(output_gdf) + decision_rows.append(decision) + omission_rows.extend(omissions) + part_rows.extend(parts) + record_omitted_candidate_trace(connection, omissions) + action_counts[decision["standard_action"]] += 1 + if group_index == 1 or group_index % 500 == 0 or group_index == len(review_groups): + print( + f"[standard] review groups {group_index:,}/{len(review_groups):,}; " + f"transformed_rows={transformed_rows:,}; omitted_duplicate_groups={omitted_duplicate_groups:,}; " + f"rss={max_rss_mb():.1f} MB", + flush=True, + ) + + print("[standard] creating source trace table inside output GeoPackage", flush=True) + trace_rows = write_source_trace_table(output_gpkg, connection) + ensure_gpkg_indexes(output_gpkg, args.output_layer) + write_source_file_audit_report(reports_dir, connection) + source_counts = { + "input_features": connection.execute("SELECT COUNT(*) FROM source_rows").fetchone()[0], + "source_identities": connection.execute("SELECT COUNT(*) FROM identities").fetchone()[0], + "candidate_features": connection.execute("SELECT COUNT(*) FROM candidate_features").fetchone()[0], + "trace_rows": trace_rows, + } + summary = { + **source_counts, + **candidate_summary, + "source_file_count": len(files), + "read_error_files": sum(1 for audit in audits if audit["read_status"] != "ok"), + "review_village_ids": len(review_ids), + "multipart_groups": len(review_groups), + "singleton_features_written": singleton_features_written, + "multipart_output_features": transformed_rows, + "omitted_features": len(omission_rows), + "omitted_duplicate_groups": omitted_duplicate_groups, + "output_features": singleton_features_written + transformed_rows, + "standard_action_counts": dict(action_counts), + "duplicate_geometry_set_signatures": len(duplicate_signatures), + "output_gpkg_bytes": output_gpkg.stat().st_size if output_gpkg.exists() else None, + "output_geojson_bytes": output_geojson.stat().st_size if output_geojson.exists() else None, + "output_csv_bytes": output_csv.stat().st_size if output_csv.exists() else None, + "max_rss_mb": max_rss_mb(), + "total_seconds": round(time.perf_counter() - started, 6), + } + if args.keep_work_db or args.work_db: + metadata["work_db"] = work_db.as_posix() + write_reports( + reports_dir, + metadata=metadata, + summary=summary, + decision_rows=decision_rows, + omission_rows=omission_rows, + part_rows=part_rows, + ) + write_gpkg_metadata(output_gpkg, "cs_admin_standard_metadata", metadata, summary) + finally: + connection.close() + if temp_dir is not None: + temp_dir.cleanup() + + print( + f"[standard] complete: {summary['output_features']:,} features written " + f"from {summary['input_features']:,} raw source rows in {summary['total_seconds']:.1f}s; " + f"rss={summary['max_rss_mb']:.1f} MB", + flush=True, + ) + return {"metadata": metadata, "summary": summary} + + +def parse_args(argv: Sequence[str] | None = None) -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--admin-input-dir", type=Path, default=DEFAULT_ADMIN_INPUT_DIR) + parser.add_argument("--multipart-analysis-csv", type=Path, default=DEFAULT_MULTIPART_ANALYSIS) + parser.add_argument("--output-gpkg", type=Path, default=DEFAULT_OUTPUT_GPKG) + parser.add_argument("--output-layer", default=DEFAULT_OUTPUT_LAYER) + parser.add_argument("--output-geojson", type=Path, default=DEFAULT_OUTPUT_GEOJSON) + parser.add_argument("--output-csv", type=Path, default=DEFAULT_OUTPUT_CSV) + parser.add_argument("--reports-dir", type=Path, default=DEFAULT_REPORTS_DIR) + parser.add_argument("--state", help="Optional state slug/name filter for smoke tests") + parser.add_argument("--district", help="Optional district slug/name filter for smoke tests") + parser.add_argument("--limit-files", type=int, help="Optional raw file limit for smoke tests") + parser.add_argument("--chunk-size", type=int, default=DEFAULT_CHUNK_SIZE) + parser.add_argument("--dominant-share-threshold", type=float, default=0.9) + parser.add_argument("--near-distance-m", type=float, default=10.0) + parser.add_argument("--far-centroid-km", type=float, default=10.0) + parser.add_argument("--many-part-threshold", type=int, default=20) + parser.add_argument("--work-db", type=Path, help="Optional SQLite work database path") + parser.add_argument("--keep-work-db", action="store_true", help="Keep the temporary work database under reports_dir") + parser.add_argument("--overwrite", action="store_true") + parser.add_argument("--skip-geojson", action="store_true") + parser.add_argument("--skip-csv", action="store_true") + parser.add_argument("--validate-only", action="store_true") + return parser.parse_args(argv) + + +def main(argv: Sequence[str] | None = None) -> None: + args = parse_args(argv) + if args.chunk_size <= 0: + raise SystemExit("--chunk-size must be positive") + if not 0 < args.dominant_share_threshold <= 1: + raise SystemExit("--dominant-share-threshold must be in (0, 1]") + build_standard_asset(args) + + +if __name__ == "__main__": + main() From 4a7dd8a0019a8d838f551e9493a8411ff9147fed Mon Sep 17 00:00:00 2001 From: amit-spatial Date: Tue, 7 Jul 2026 04:45:56 +0000 Subject: [PATCH 2/3] Document CS admin standard build decisions The workplan captures the data-quality rules, multipart-village handling, and lineage expectations behind cs_admin_standard. Keeping it beside the builder makes the asset reproducible and gives reviewers the rationale for non-trivial geometry standardisation decisions. --- .../cs_admin_boundary_standard_workplan.md | 134 ++++++++++++++++++ 1 file changed, 134 insertions(+) create mode 100644 utilities/scripts/admin_assets/cs_admin_boundary_standard_workplan.md diff --git a/utilities/scripts/admin_assets/cs_admin_boundary_standard_workplan.md b/utilities/scripts/admin_assets/cs_admin_boundary_standard_workplan.md new file mode 100644 index 00000000..f1680a05 --- /dev/null +++ b/utilities/scripts/admin_assets/cs_admin_boundary_standard_workplan.md @@ -0,0 +1,134 @@ +# CS Admin Boundary Standard Workplan + +Date: 2026-07-05 +Working branch at start: `feat/facilities_v2` + +## Operating Constraints + +- Keep RAM usage under control during geospatial builds. Avoid running two large full-India GeoPandas/Fiona writes at the same time. +- Prefer `uv run --with ...` for Python/data runs when dependencies are needed. +- Preserve existing pipeline scripts unless the task explicitly calls for a replacement script. +- Add new scripts for the new standard path: + - `utilities/scripts/admin_assets/build_cs_admin_boundary_standard.py` + - `utilities/scripts/admin_assets/build_cs_admin_boundary_assets.py` +- Keep `data/admin-boundary/cs_admin_standard.gpkg` as the primary standard asset. +- Also export: + - `data/admin-boundary/cs_admin_standard.geojson` + - `data/admin-boundary/cs_admin_standard.csv` +- Maintain consistent village identifiers and dtypes across admin, facilities, livestock, and Antyodaya outputs. + +## Admin Boundary Standard + +Primary inputs to inspect before implementation: + +- `data/admin-boundary/input/` +- `data/admin-boundary/multi_part_village_analysis.csv` +- `data/admin-boundary/multi_part_village_analysis_file_data.csv` +- `data/admin-boundary/cs_admin_boundary_village_shapes_analysis.py` +- `data/admin-boundary/cs_admin_boundary_village_shapes_analysis.geojson` +- `utilities/scripts/admin_assets/build_admin_boundary_assets.py` + +`data/admin-boundary/cs_admin_sanitised.gpkg` is useful as a previous-build reference only. The new standard asset must maintain traceability from raw `data/admin-boundary/input/` rows to `data/admin-boundary/cs_admin_standard.gpkg`, not from `cs_admin_sanitised.gpkg`. + +Expected behavior: + +- Produce a reproducible standardized village asset with traceable decisions. +- Log omissions, duplicate handling, multipart handling, geometry repair, dtype coercion, and raw source rows. +- Store raw source trace in the output GeoPackage through a non-spatial table. +- Keep multipart villages when the analysis notes or geometry context indicate they are valid distinct/urban/town-like cases. +- Merge or dissolve only when rules and notes make it safe. +- Keep part indexes sorted by descending area share when multipart geometry remains. + +## Facilities + +Branch-aware work: + +- Facilities rebuild work spans `feat/facilities_v2` and `feat/facilities_data_processing`. +- Rebuild proximity using: + - `data/facilities/outputs/pan_india_facilities.gpkg` + - `data/admin-boundary/cs_admin_standard.gpkg` +- Target proximity output: + - `data/facilities/outputs/cs_village_facility_proximity.gpkg` + - plus a useful CSV sibling for local inspection/API parity. +- Update `computing/misc/facilities_proximity.py` only in the branch where that runtime API change belongs. + +## Joined Admin Assets + +Create a revised builder from the useful patterns in `build_admin_boundary_assets.py`, but do not depend on the old script at runtime. + +Targets: + +- With `utilities/scripts/admin_assets/asset_configs/livestock.json`, rebuild: + - `data/livestock/cs_village_livestock_census_20.gpkg` + - `data/livestock/cs_village_livestock_census_20.geojson` +- With `utilities/scripts/admin_assets/asset_configs/antyodaya.json`, rebuild: + - `data/antyodaya/output/cs_antyodaya_2020.gpkg` + - `data/antyodaya/output/cs_antyodaya_2020.geojson` + +The new builder resolves legacy config references to +`data/admin-boundary/cs_admin_standard.gpkg` / layer `cs_admin_standard` by +default. It prepends the standard identity columns (`cs_feature_id`, +`cs_admin_uid`, `core_admin_uid`, `pc11_village_id`, `village_id`) before the +configured joined columns and writes unmatched-admin and unmatched-source CSV +reports beside each asset. + +## Monitoring + +For long-running steps, run one heavy geospatial build at a time unless a companion task is clearly lightweight. Check memory/process status periodically with commands such as: + +```bash +ps -eo pid,ppid,pmem,rss,etime,cmd --sort=-rss | head -20 +free -h +``` + +Use concise monitoring notes and avoid dumping full logs unless diagnosing a failure. + +Current full admin-standard build is running comfortably below the original +RAM ceiling. If a future rebuild is bottlenecked mostly by conservative +chunking, it is acceptable to tune chunk sizes for a managed target around +4 GB RAM, while still avoiding simultaneous full-India geometry writes. + +External data-preparation scripts may use `pyogrio`/GeoPandas for speed. The +Django/API facilities runtime in `computing/misc/facilities_proximity.py` +should not introduce a `pyogrio` dependency. + +## 2026-07-05 Build Checkpoint + +Completed and validated: + +- `data/admin-boundary/cs_admin_standard.gpkg` + - Layer rows: `630,898` + - Raw source-trace rows: `983,551` + - Trace status counts: `971,657` output, `11,894` omitted + - Source files read OK: `660/660` + - `pc11_village_id = 800372`: exactly `1` retained ID row; `96` + geometries retained with the village id removed as invalid fragments + - SQLite `quick_check`: `ok` + - Peak logged RSS: `3,876.5 MB` +- `data/livestock/cs_village_livestock_census_20.gpkg` + - Layer rows: `630,898` + - Matched admin rows: `474,240` + - SQLite `quick_check`: `ok` + - GeoJSON sibling written +- `data/antyodaya/output/cs_antyodaya_2020.gpkg` + - Layer rows: `630,898` + - Matched admin rows: `545,874` + - SQLite `quick_check`: `ok` + - GeoJSON sibling written + +Facilities correction: + +- `data/facilities/outputs/cs_village_facility_proximity.gpkg` + - The final village feature layer is `village_shapes`, not `village_points`. + - `village_shapes` geometry type: `MULTIPOLYGON`, CRS `EPSG:4326` + - `village_shapes` rows: `630,898` + - `proximity_l3` rows: `15,772,450` + - `proximity_l2_materialized` rows: `6,939,878` + - `proximity_l1_materialized` rows: `3,154,490` + - `proximity_nearest_facilities` rows: `1,711,631` + - SQLite `quick_check`: `ok` +- `data/facilities/outputs/cs_village_facility_proximity_summary.csv` + - 40 rows: table counts, 25 L3 class summaries, metadata +- Smoke build from `utilities/scripts/facilities_utils/facility_proximity_finder.py` + writes `village_shapes` as `MULTIPOLYGON` and does not create a + `village_points` layer. From 8269523e4c6835f4932dda0f58da9383a5f5840f Mon Sep 17 00:00:00 2001 From: amit-spatial Date: Mon, 20 Jul 2026 11:57:50 +0000 Subject: [PATCH 3/3] Refresh changed CS admin standard inputs --- .../build_cs_admin_boundary_standard.py | 2 - .../refresh_cs_admin_boundary_standard.py | 400 ++++++++++++++++++ 2 files changed, 400 insertions(+), 2 deletions(-) create mode 100644 utilities/scripts/admin_assets/refresh_cs_admin_boundary_standard.py diff --git a/utilities/scripts/admin_assets/build_cs_admin_boundary_standard.py b/utilities/scripts/admin_assets/build_cs_admin_boundary_standard.py index 64cff071..84601d99 100644 --- a/utilities/scripts/admin_assets/build_cs_admin_boundary_standard.py +++ b/utilities/scripts/admin_assets/build_cs_admin_boundary_standard.py @@ -20,7 +20,6 @@ import argparse from collections import Counter -import csv import hashlib import json import math @@ -2209,7 +2208,6 @@ def write_reports( for flag in str(row.get("pattern_flags") or "").split("|"): if flag: pattern_counts[flag] += 1 - action_counts = Counter(row["standard_action"] for row in decision_rows) pd.DataFrame( [{"pattern_flag": key, "groups": value} for key, value in sorted(pattern_counts.items())] ).to_csv(reports_dir / "standard_pattern_summary.csv", index=False) diff --git a/utilities/scripts/admin_assets/refresh_cs_admin_boundary_standard.py b/utilities/scripts/admin_assets/refresh_cs_admin_boundary_standard.py new file mode 100644 index 00000000..ce548e79 --- /dev/null +++ b/utilities/scripts/admin_assets/refresh_cs_admin_boundary_standard.py @@ -0,0 +1,400 @@ +#!/usr/bin/env python3 +"""Refresh changed district rows in the standard admin GeoPackage. + +The first run treats input files newer than the existing GeoPackage as changed. +Successful runs persist an input manifest inside the GeoPackage, so later runs +also detect new, removed, or replaced files. Use ``--force-rebuild`` when a +global standardisation rule or the multipart review file has changed. +""" + +from __future__ import annotations + +import argparse +import json +import os +from pathlib import Path +import shutil +import sqlite3 +import tempfile +import time +from typing import Any, Sequence + +import pyogrio + +from build_cs_admin_boundary_standard import ( + DEFAULT_ADMIN_INPUT_DIR, + DEFAULT_MULTIPART_ANALYSIS, + DEFAULT_OUTPUT_GPKG, + DEFAULT_OUTPUT_LAYER, + DEFAULT_REPORTS_DIR, + ROOT_DIR, + build_standard_asset, + normalize_slug, + parse_args as parse_builder_args, + repo_path, +) + + +DEFAULT_BASE_RESOURCE = Path("data/base_resources/cs_admin_standard.gpkg") +MANIFEST_TABLE = "cs_admin_standard_input_manifest" +REFRESH_TABLE = "cs_admin_standard_refresh_log" +TRACE_TABLE = "cs_admin_standard_source_trace" + + +def source_file_name(path: Path) -> str: + return path.relative_to(ROOT_DIR).as_posix() if path.is_relative_to(ROOT_DIR) else path.as_posix() + + +def current_inputs(input_dir: Path) -> dict[str, dict[str, Any]]: + inputs: dict[str, dict[str, Any]] = {} + for path in sorted(input_dir.glob("*/*.geojson")): + stat = path.stat() + inputs[source_file_name(path)] = { + "source_file": source_file_name(path), + "source_state_slug": normalize_slug(path.parent.name), + "source_district_slug": normalize_slug(path.stem), + "size_bytes": stat.st_size, + "mtime_ns": stat.st_mtime_ns, + } + return inputs + + +def load_manifest(gpkg_path: Path) -> dict[str, dict[str, Any]]: + if not gpkg_path.exists(): + return {} + with sqlite3.connect(gpkg_path) as connection: + exists = connection.execute( + "SELECT 1 FROM sqlite_master WHERE type = 'table' AND name = ?", + (MANIFEST_TABLE,), + ).fetchone() + if not exists: + return {} + columns = [row[1] for row in connection.execute(f'PRAGMA table_info("{MANIFEST_TABLE}")')] + return { + row[0]: dict(zip(columns, row)) + for row in connection.execute(f'SELECT * FROM "{MANIFEST_TABLE}"') + } + + +def changed_inputs( + current: dict[str, dict[str, Any]], + previous: dict[str, dict[str, Any]], + gpkg_path: Path, +) -> tuple[list[dict[str, Any]], list[dict[str, Any]]]: + if previous: + changed = [ + item + for name, item in current.items() + if name not in previous + or int(previous[name]["size_bytes"]) != item["size_bytes"] + or int(previous[name]["mtime_ns"]) != item["mtime_ns"] + ] + removed = [item for name, item in previous.items() if name not in current] + return changed, removed + + output_mtime = gpkg_path.stat().st_mtime_ns if gpkg_path.exists() else 0 + return [item for item in current.values() if item["mtime_ns"] > output_mtime], [] + + +def write_manifest(connection: sqlite3.Connection, inputs: dict[str, dict[str, Any]]) -> None: + connection.execute(f'DROP TABLE IF EXISTS "{MANIFEST_TABLE}"') + connection.execute( + f'''CREATE TABLE "{MANIFEST_TABLE}" ( + source_file TEXT PRIMARY KEY, + source_state_slug TEXT NOT NULL, + source_district_slug TEXT NOT NULL, + size_bytes INTEGER NOT NULL, + mtime_ns INTEGER NOT NULL + )''' + ) + connection.executemany( + f'INSERT INTO "{MANIFEST_TABLE}" VALUES (?, ?, ?, ?, ?)', + [ + ( + item["source_file"], + item["source_state_slug"], + item["source_district_slug"], + item["size_bytes"], + item["mtime_ns"], + ) + for item in inputs.values() + ], + ) + + +def record_refresh(connection: sqlite3.Connection, details: dict[str, Any]) -> None: + connection.execute( + f'''CREATE TABLE IF NOT EXISTS "{REFRESH_TABLE}" ( + refreshed_at_utc TEXT NOT NULL, + mode TEXT NOT NULL, + changed_files INTEGER NOT NULL, + removed_files INTEGER NOT NULL, + details_json TEXT NOT NULL + )''' + ) + connection.execute( + f'INSERT INTO "{REFRESH_TABLE}" VALUES (datetime("now"), ?, ?, ?, ?)', + ( + details["mode"], + len(details.get("changed_files", [])), + len(details.get("removed_files", [])), + json.dumps(details, sort_keys=True), + ), + ) + + +def copy_atomically(source: Path, destination: Path) -> None: + destination.parent.mkdir(parents=True, exist_ok=True) + with tempfile.NamedTemporaryFile( + prefix=f".{destination.name}.", suffix=".tmp", dir=destination.parent, delete=False + ) as handle: + temporary = Path(handle.name) + try: + shutil.copy2(source, temporary) + os.replace(temporary, destination) + finally: + temporary.unlink(missing_ok=True) + + +def wait_for_file(path: Path, *, timeout_seconds: float = 30.0) -> None: + """Wait for a replaced file to become visible on mounted/network filesystems.""" + + deadline = time.monotonic() + timeout_seconds + while time.monotonic() < deadline: + if path.exists(): + return + time.sleep(0.25) + raise FileNotFoundError(f"Replaced file did not become visible within {timeout_seconds}s: {path}") + + +def layer_labels_for_scope( + connection: sqlite3.Connection, + state_slug: str, + district_slug: str, +) -> set[tuple[str, str]]: + return { + (state_name, district_name) + for state_name, district_name in connection.execute( + f'''SELECT DISTINCT state_name, district_name + FROM "{TRACE_TABLE}" + WHERE source_state_slug = ? AND source_district_slug = ? + AND state_name IS NOT NULL AND district_name IS NOT NULL''', + (state_slug, district_slug), + ) + } + + +def replace_scope( + target: Path, + replacement: Path | None, + *, + layer: str, + state_slug: str, + district_slug: str, +) -> tuple[int, int]: + with sqlite3.connect(target) as connection: + labels = layer_labels_for_scope(connection, state_slug, district_slug) + + replacement_frame = None + replacement_trace_rows: list[tuple[Any, ...]] = [] + trace_columns: list[str] = [] + if replacement is not None: + replacement_frame = pyogrio.read_dataframe(replacement, layer=layer) + labels.update( + (str(state_name), str(district_name)) + for state_name, district_name in replacement_frame[["state_name", "district_name"]] + .drop_duplicates() + .itertuples(index=False, name=None) + ) + with sqlite3.connect(replacement) as source: + trace_columns = [row[1] for row in source.execute(f'PRAGMA table_info("{TRACE_TABLE}")')] + replacement_trace_rows = source.execute(f'SELECT * FROM "{TRACE_TABLE}"').fetchall() + + with sqlite3.connect(target) as connection: + deleted_rows = 0 + for state_name, district_name in labels: + cursor = connection.execute( + f'DELETE FROM "{layer}" WHERE state_name = ? AND district_name = ?', + (state_name, district_name), + ) + deleted_rows += cursor.rowcount + connection.execute( + f'DELETE FROM "{TRACE_TABLE}" WHERE source_state_slug = ? AND source_district_slug = ?', + (state_slug, district_slug), + ) + connection.commit() + + inserted_rows = 0 + if replacement_frame is not None and not replacement_frame.empty: + pyogrio.write_dataframe( + replacement_frame, + target, + layer=layer, + driver="GPKG", + append=True, + ) + inserted_rows = len(replacement_frame) + if replacement_trace_rows: + quoted = ", ".join(f'"{column}"' for column in trace_columns) + placeholders = ", ".join("?" for _ in trace_columns) + with sqlite3.connect(target) as connection: + connection.executemany( + f'INSERT INTO "{TRACE_TABLE}" ({quoted}) VALUES ({placeholders})', + replacement_trace_rows, + ) + connection.commit() + return deleted_rows, inserted_rows + + +def build_scope( + input_dir: Path, + multipart_csv: Path, + state_slug: str, + district_slug: str, + output_dir: Path, + layer: str, +) -> Path: + output = output_dir / f"{state_slug.lower()}_{district_slug.lower()}.gpkg" + args = parse_builder_args( + [ + "--admin-input-dir", str(input_dir), + "--multipart-analysis-csv", str(multipart_csv), + "--output-gpkg", str(output), + "--output-layer", layer, + "--reports-dir", str(output_dir / "reports" / state_slug.lower() / district_slug.lower()), + "--state", state_slug, + "--district", district_slug, + "--skip-geojson", + "--skip-csv", + "--overwrite", + ] + ) + build_standard_asset(args) + return output + + +def force_rebuild(args: argparse.Namespace, inputs: dict[str, dict[str, Any]]) -> dict[str, Any]: + target = repo_path(args.output_gpkg) + target.parent.mkdir(parents=True, exist_ok=True) + with tempfile.TemporaryDirectory(prefix="cs_admin_full_rebuild_", dir=target.parent) as directory: + staging = Path(directory) / target.name + builder_args = parse_builder_args( + [ + "--admin-input-dir", str(args.admin_input_dir), + "--multipart-analysis-csv", str(args.multipart_analysis_csv), + "--output-gpkg", str(staging), + "--output-layer", args.output_layer, + "--reports-dir", str(args.reports_dir), + "--skip-geojson", + "--skip-csv", + "--overwrite", + ] + ) + result = build_standard_asset(builder_args) + details = {"mode": "force_rebuild", "changed_files": sorted(inputs), "removed_files": []} + with sqlite3.connect(staging) as connection: + write_manifest(connection, inputs) + record_refresh(connection, details) + connection.commit() + os.replace(staging, target) + wait_for_file(target) + copy_atomically(target, repo_path(args.base_resource)) + return {**details, "summary": result["summary"]} + + +def refresh(args: argparse.Namespace) -> dict[str, Any]: + input_dir = repo_path(args.admin_input_dir) + multipart_csv = repo_path(args.multipart_analysis_csv) + target = repo_path(args.output_gpkg) + current = current_inputs(input_dir) + previous = load_manifest(target) + changed, removed = changed_inputs(current, previous, target) + + if args.force_rebuild or not target.exists(): + return force_rebuild(args, current) + + scopes = sorted( + { + (item["source_state_slug"], item["source_district_slug"]) + for item in [*changed, *removed] + } + ) + details: dict[str, Any] = { + "mode": "incremental" if scopes else "unchanged", + "changed_files": sorted(item["source_file"] for item in changed), + "removed_files": sorted(item["source_file"] for item in removed), + "scopes": [], + } + if args.dry_run: + details["scopes"] = [f"{state}/{district}" for state, district in scopes] + return details + + if scopes: + with tempfile.TemporaryDirectory(prefix="cs_admin_refresh_", dir=target.parent) as directory: + temporary_dir = Path(directory) + staging = temporary_dir / target.name + shutil.copy2(target, staging) + for state_slug, district_slug in scopes: + scope_exists = any( + item["source_state_slug"] == state_slug + and item["source_district_slug"] == district_slug + for item in current.values() + ) + replacement = ( + build_scope(input_dir, multipart_csv, state_slug, district_slug, temporary_dir, args.output_layer) + if scope_exists + else None + ) + deleted, inserted = replace_scope( + staging, + replacement, + layer=args.output_layer, + state_slug=state_slug, + district_slug=district_slug, + ) + details["scopes"].append( + { + "state": state_slug, + "district": district_slug, + "deleted_rows": deleted, + "inserted_rows": inserted, + } + ) + with sqlite3.connect(staging) as connection: + write_manifest(connection, current) + record_refresh(connection, details) + connection.commit() + os.replace(staging, target) + wait_for_file(target) + elif not previous: + with sqlite3.connect(target) as connection: + write_manifest(connection, current) + record_refresh(connection, details) + connection.commit() + + copy_atomically(target, repo_path(args.base_resource)) + return details + + +def parse_args(argv: Sequence[str] | None = None) -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--admin-input-dir", type=Path, default=DEFAULT_ADMIN_INPUT_DIR) + parser.add_argument("--multipart-analysis-csv", type=Path, default=DEFAULT_MULTIPART_ANALYSIS) + parser.add_argument("--output-gpkg", type=Path, default=DEFAULT_OUTPUT_GPKG) + parser.add_argument("--output-layer", default=DEFAULT_OUTPUT_LAYER) + parser.add_argument("--reports-dir", type=Path, default=DEFAULT_REPORTS_DIR) + parser.add_argument("--base-resource", type=Path, default=DEFAULT_BASE_RESOURCE) + parser.add_argument("--force-rebuild", action="store_true", help="Rebuild the full national GeoPackage") + parser.add_argument("--dry-run", action="store_true", help="Report detected changes without writing outputs") + return parser.parse_args(argv) + + +def main(argv: Sequence[str] | None = None) -> None: + started = time.perf_counter() + result = refresh(parse_args(argv)) + result["elapsed_seconds"] = round(time.perf_counter() - started, 3) + print(json.dumps(result, indent=2, sort_keys=True)) + + +if __name__ == "__main__": + main()