Skills › Scientific Agent Skills › geopandas › spatial_join_audit.py
spatial_join_audit.py
375 lines · skills/geopandas/scripts/spatial_join_audit.py at 5f2bc03
1#!/usr/bin/env python3
2"""Run a bounded local spatial join and emit only aggregate cardinality."""
4from __future__ import annotations
6import argparse
7import sys
8from collections import Counter
9from typing import Any
11from _common import (
12 ABSOLUTE_MAX_FEATURES,
13 DEFAULT_MAX_FEATURES,
14 DEFAULT_MAX_INPUT_BYTES,
15 CliError,
16 bounded_limit,
17 checked_input_file,
18 crs_summary,
19 duplicate_column_state,
20 emit_json,
21 fail_json,
22 finite_number,
23 geometry_state,
24 load_geodataframe,
25 native_versions,
26 package_versions,
27 positive_int,
28 sha256_file,
29)
31TOOL = "spatial_join_audit"
32DEFAULT_MAX_PAIRS = 1_000_000
33ABSOLUTE_MAX_PAIRS = 10_000_000
36def build_parser() -> argparse.ArgumentParser:
37 parser = argparse.ArgumentParser(
38 description=(
39 "Audit a bounded local binary-predicate or nearest spatial join. "
40 "Only aggregate multiplicity is emitted; no pairs, IDs, or coordinates."
41 ),
42 epilog=(
43 "Both inputs must be local allowlisted files under --root. Invalid or "
44 "CRS-incompatible inputs block execution. All operations are planar."
45 ),
46 )
47 parser.add_argument("left", help="Local left vector file")
48 parser.add_argument("right", help="Local right vector file")
49 parser.add_argument("--root", default=".", help="Existing local I/O root")
50 parser.add_argument("--left-layer", help="Exact left layer name")
51 parser.add_argument("--right-layer", help="Exact right layer name")
52 parser.add_argument("--left-id", help="Optional left stable ID column to audit")
53 parser.add_argument("--right-id", help="Optional right stable ID column to audit")
54 parser.add_argument(
55 "--mode",
56 choices=("predicate", "nearest"),
57 default="predicate",
58 help="Join mode (default: predicate)",
59 )
60 parser.add_argument(
61 "--predicate",
62 default="intersects",
63 help="Spatial-index predicate for predicate mode (default: intersects)",
64 )
65 parser.add_argument(
66 "--distance",
67 help="Positive CRS-unit scalar required for predicate=dwithin",
68 )
69 parser.add_argument(
70 "--max-distance",
71 help="Optional positive CRS-unit search limit for nearest mode",
72 )
73 parser.add_argument(
74 "--exclusive",
75 action="store_true",
76 help="Nearest mode: exclude geometrically equal candidates",
77 )
78 parser.add_argument(
79 "--on-attribute",
80 action="append",
81 default=[],
82 help="Additional equality column present on both sides; repeat as needed",
83 )
84 parser.add_argument(
85 "--max-input-bytes",
86 type=positive_int,
87 default=DEFAULT_MAX_INPUT_BYTES,
88 help="Maximum bytes for each primary input file",
89 )
90 parser.add_argument(
91 "--max-features",
92 type=positive_int,
93 default=DEFAULT_MAX_FEATURES,
94 help=f"Maximum features per input (default: {DEFAULT_MAX_FEATURES})",
95 )
96 parser.add_argument(
97 "--max-pairs",
98 type=positive_int,
99 default=DEFAULT_MAX_PAIRS,
100 help=f"Maximum joined pairs before failing closed (default: {DEFAULT_MAX_PAIRS})",
101 )
102 return parser
105def _work_frame(frame: Any, row_column: str, attributes: list[str]) -> Any:
106 active = frame.geometry.name
107 columns = [active, *attributes]
108 work = frame.loc[:, columns].reset_index(drop=True).copy()
109 if active != "_audit_geometry":
110 work = work.rename_geometry("_audit_geometry")
111 work[row_column] = range(len(work))
112 return work
115def _counter_summary(counter: Counter[int], total_features: int) -> dict[str, int]:
116 matched = len(counter)
117 return {
118 "matched_features": matched,
119 "unmatched_features": total_features - matched,
120 "features_with_one_match": sum(count == 1 for count in counter.values()),
121 "features_with_multiple_matches": sum(
122 count > 1 for count in counter.values()
123 ),
124 "maximum_matches_for_one_feature": max(counter.values(), default=0),
125 }
128def _run_chunked_join(
129 left: Any,
130 right: Any,
131 *,
132 mode: str,
133 predicate: str,
134 distance: float | None,
135 max_distance: float | None,
136 exclusive: bool,
137 attributes: list[str],
138 max_pairs: int,
139) -> dict[str, Any]:
140 try:
141 import geopandas as gpd
142 except ImportError as exc:
143 raise CliError("GeoPandas is required for join analysis") from exc
145 left_work = _work_frame(left, "_audit_left_row", attributes)
146 right_work = _work_frame(right, "_audit_right_row", attributes)
147 if len(right_work) > max_pairs:
148 raise CliError(
149 "right feature count exceeds max_pairs; one left feature could match all"
150 )
151 left_counts: Counter[int] = Counter()
152 right_counts: Counter[int] = Counter()
153 pair_count = 0
154 right_size = max(len(right_work), 1)
155 chunk_size = max(1, min(512, max_pairs // right_size or 1))
157 for start in range(0, len(left_work), chunk_size):
158 chunk = left_work.iloc[start : start + chunk_size]
159 if mode == "nearest":
160 joined = gpd.sjoin_nearest(
161 chunk,
162 right_work,
163 how="inner",
164 max_distance=max_distance,
165 exclusive=exclusive,
166 )
167 else:
168 joined = gpd.sjoin(
169 chunk,
170 right_work,
171 how="inner",
172 predicate=predicate,
173 distance=distance,
174 on_attribute=attributes or None,
175 )
176 pair_count += len(joined)
177 if pair_count > max_pairs:
178 raise CliError(
179 f"joined pair count exceeds the configured {max_pairs}-pair limit"
180 )
181 left_counts.update(int(value) for value in joined["_audit_left_row"])
182 right_counts.update(int(value) for value in joined["_audit_right_row"])
184 return {
185 "pair_count": pair_count,
186 "left": _counter_summary(left_counts, len(left)),
187 "right": _counter_summary(right_counts, len(right)),
188 "many_to_many_observed": bool(
189 any(count > 1 for count in left_counts.values())
190 and any(count > 1 for count in right_counts.values())
191 ),
192 "pairs_emitted": False,
193 "chunked": True,
194 "chunk_size": chunk_size,
195 }
198def audit(args: argparse.Namespace) -> dict[str, Any]:
199 max_features = bounded_limit(
200 args.max_features,
201 name="max_features",
202 maximum=ABSOLUTE_MAX_FEATURES,
203 )
204 max_pairs = bounded_limit(
205 args.max_pairs,
206 name="max_pairs",
207 maximum=ABSOLUTE_MAX_PAIRS,
208 )
209 left_path = checked_input_file(
210 args.left,
211 root=args.root,
212 max_bytes=args.max_input_bytes,
213 )
214 right_path = checked_input_file(
215 args.right,
216 root=args.root,
217 max_bytes=args.max_input_bytes,
218 )
219 if left_path == right_path and args.left_layer == args.right_layer:
220 raise CliError("left and right inputs resolve to the same layer")
222 left = load_geodataframe(
223 left_path,
224 layer=args.left_layer,
225 max_features=max_features,
226 )
227 right = load_geodataframe(
228 right_path,
229 layer=args.right_layer,
230 max_features=max_features,
231 )
232 left_state = geometry_state(left)
233 right_state = geometry_state(right)
234 left_id = duplicate_column_state(left, args.left_id)
235 right_id = duplicate_column_state(right, args.right_id)
237 try:
238 from pyproj import CRS
239 except ImportError as exc:
240 raise CliError("pyproj is required for CRS comparison") from exc
241 blockers: list[str] = []
242 same_crs = False
243 if left.crs is None or right.crs is None:
244 blockers.append("both inputs require explicit CRS metadata")
245 else:
246 same_crs = CRS.from_user_input(left.crs).equals(CRS.from_user_input(right.crs))
247 if not same_crs:
248 blockers.append("input CRS values are not equivalent")
249 if left_state["invalid"] or right_state["invalid"]:
250 blockers.append("invalid input geometries require separate repair review")
251 for attribute in args.on_attribute:
252 if attribute not in left.columns or attribute not in right.columns:
253 blockers.append("an on_attribute column is absent from one input")
254 break
255 if attribute in {
256 "_audit_geometry", "_audit_left_row", "_audit_right_row",
257 left.geometry.name, right.geometry.name,
258 }:
259 blockers.append("on_attribute must not use geometry or reserved audit columns")
260 if len(set(args.on_attribute)) != len(args.on_attribute):
261 blockers.append("on_attribute columns must not be repeated")
263 distance: float | None = None
264 max_distance: float | None = None
265 if args.mode == "predicate":
266 if args.exclusive or args.max_distance is not None:
267 blockers.append("--exclusive/--max-distance apply only to nearest mode")
268 if args.predicate == "dwithin":
269 if args.distance is None:
270 blockers.append("predicate dwithin requires --distance")
271 else:
272 distance = finite_number(args.distance, name="distance")
273 if distance <= 0:
274 blockers.append("distance must be greater than zero")
275 elif args.distance is not None:
276 blockers.append("--distance is valid only for predicate dwithin")
277 if not blockers:
278 valid = set(left.sindex.valid_query_predicates) & set(
279 right.sindex.valid_query_predicates
280 )
281 if args.predicate not in valid:
282 blockers.append("predicate is not supported by both spatial indexes")
283 else:
284 if args.distance is not None or args.on_attribute:
285 blockers.append("--distance/--on-attribute apply only to predicate mode")
286 if args.max_distance is not None:
287 max_distance = finite_number(args.max_distance, name="max_distance")
288 if max_distance <= 0:
289 blockers.append("max_distance must be greater than zero")
291 geographic = bool(left.crs is not None and left.crs.is_geographic)
292 if args.mode == "nearest" and geographic:
293 blockers.append("nearest joins are inaccurate in a geographic CRS")
294 if args.predicate == "dwithin" and geographic:
295 blockers.append("dwithin distance is angular in a geographic CRS")
297 pair_audit = None
298 if not blockers:
299 pair_audit = _run_chunked_join(
300 left,
301 right,
302 mode=args.mode,
303 predicate=args.predicate,
304 distance=distance,
305 max_distance=max_distance,
306 exclusive=args.exclusive,
307 attributes=args.on_attribute,
308 max_pairs=max_pairs,
309 )
311 warnings = [
312 "intersects includes boundary contact; contains/within and covers/covered_by differ at boundaries.",
313 "All GeoPandas joins are planar and ignore Z.",
314 "Declare expected cardinality separately and compare it with this aggregate audit.",
315 ]
316 if args.mode == "nearest":
317 warnings.append(
318 "Nearest returns every equidistant match and has no k parameter."
319 )
320 if left_state["missing"] or left_state["empty"]:
321 warnings.append("Missing/empty left geometries remain unmatched.")
322 if right_state["missing"] or right_state["empty"]:
323 warnings.append("Missing/empty right geometries remain unmatched.")
325 return {
326 "ok": not blockers,
327 "tool": TOOL,
328 "mode": args.mode,
329 "predicate": args.predicate if args.mode == "predicate" else None,
330 "distance_provided": distance is not None,
331 "max_distance_provided": max_distance is not None,
332 "exclusive": bool(args.exclusive),
333 "on_attribute_column_count": len(args.on_attribute),
334 "crs": {
335 "left": crs_summary(left.crs),
336 "right": crs_summary(right.crs),
337 "equivalent": same_crs,
338 },
339 "left": {
340 "source_sha256": sha256_file(left_path),
341 "geometry_state": left_state,
342 "stable_id_audit": left_id,
343 },
344 "right": {
345 "source_sha256": sha256_file(right_path),
346 "geometry_state": right_state,
347 "stable_id_audit": right_id,
348 },
349 "pair_audit": pair_audit,
350 "blockers": blockers,
351 "resource_limits": {
352 "max_input_bytes_each": args.max_input_bytes,
353 "max_features_each": max_features,
354 "max_pairs": max_pairs,
355 },
356 "network_accessed": False,
357 "coordinates_emitted": False,
358 "identifiers_emitted": False,
359 "stack": {"packages": package_versions(), "native": native_versions()},
360 "warnings": warnings,
361 }
364def main(argv: list[str] | None = None) -> int:
365 args = build_parser().parse_args(argv)
366 try:
367 report = audit(args)
368 emit_json(report)
369 return 0 if report["ok"] else 2
370 except Exception as exc: # noqa: BLE001 - errors are redacted at CLI boundary
371 return fail_json(TOOL, exc)
374if __name__ == "__main__":
375 sys.exit(main())