Skills › Scientific Agent Skills › geopandas › crs_reprojection_plan.py
crs_reprojection_plan.py
223 lines · skills/geopandas/scripts/crs_reprojection_plan.py at 5f2bc03
1#!/usr/bin/env python3
2"""Plan CRS and datum transformation semantics without transforming coordinates."""
4from __future__ import annotations
6import argparse
7import sys
8from typing import Any
10from _common import (
11 CliError,
12 crs_summary,
13 disable_proj_network,
14 emit_json,
15 fail_json,
16 finite_number,
17 native_versions,
18 package_versions,
19)
21TOOL = "crs_reprojection_plan"
22METRIC_OPERATIONS = {"area", "buffer", "distance", "nearest", "precision"}
25def build_parser() -> argparse.ArgumentParser:
26 parser = argparse.ArgumentParser(
27 description=(
28 "Inspect CRS axes, units, datum-operation candidates, grid availability, "
29 "and antimeridian risk without reading data or transforming coordinates."
30 ),
31 epilog=(
32 "PROJ network access is forcibly disabled. Optional bbox values are used "
33 "only as a protected area-of-interest hint and are never emitted."
34 ),
35 )
36 parser.add_argument("--source-crs", required=True, help="Authoritative source CRS")
37 parser.add_argument("--target-crs", required=True, help="Proposed target CRS")
38 parser.add_argument(
39 "--operation",
40 choices=(
41 "reproject",
42 "area",
43 "buffer",
44 "distance",
45 "nearest",
46 "precision",
47 "display",
48 ),
49 default="reproject",
50 help="Purpose used for unit-readiness checks (default: reproject)",
51 )
52 parser.add_argument(
53 "--bbox",
54 nargs=4,
55 metavar=("WEST_DEG", "SOUTH_DEG", "EAST_DEG", "NORTH_DEG"),
56 help=(
57 "Optional geographic source bbox in longitude/latitude degrees. "
58 "Projected sources are rejected; east < west marks an antimeridian "
59 "crossing. Values are not emitted."
60 ),
61 )
62 parser.add_argument(
63 "--allow-ballpark",
64 action="store_true",
65 help="Include ballpark candidates in this plan (not recommended for accuracy work)",
66 )
67 return parser
70def _dynamic(crs: Any) -> bool:
71 datum = getattr(crs, "datum", None)
72 type_name = str(getattr(datum, "type_name", "")).casefold()
73 return "dynamic" in type_name
76def _candidate_summary(transformer: Any) -> dict[str, Any]:
77 accuracy = float(transformer.accuracy)
78 area = transformer.area_of_use
79 return {
80 "description": str(transformer.description)[:300],
81 "accuracy_metres": accuracy if accuracy >= 0 else None,
82 "accuracy_known": accuracy >= 0,
83 "area_of_use_name": str(area.name)[:200] if area else None,
84 }
87def plan(args: argparse.Namespace) -> dict[str, Any]:
88 try:
89 from pyproj import CRS
90 from pyproj.aoi import AreaOfInterest
91 from pyproj.transformer import TransformerGroup
92 except ImportError as exc:
93 raise CliError("pyproj is required for CRS planning") from exc
95 disable_proj_network()
96 try:
97 source = CRS.from_user_input(args.source_crs)
98 target = CRS.from_user_input(args.target_crs)
99 except (TypeError, ValueError) as exc:
100 raise CliError("source or target CRS cannot be parsed") from exc
102 bbox_values: list[float] | None = None
103 crosses_antimeridian = False
104 area_of_interest = None
105 if args.bbox is not None:
106 if not source.is_geographic or not all(
107 str(axis.unit_name).casefold() == "degree"
108 for axis in source.axis_info[:2]
109 ):
110 raise CliError(
111 "--bbox requires a geographic source CRS with degree axes; "
112 "omit it for projected sources rather than treating metres/feet as degrees"
113 )
114 bbox_values = [
115 finite_number(item, name=f"bbox[{index}]")
116 for index, item in enumerate(args.bbox)
117 ]
118 west, south, east, north = bbox_values
119 if not (-180 <= west <= 180 and -180 <= east <= 180):
120 raise CliError("bbox longitude must be within [-180, 180] degrees")
121 if not (-90 <= south <= 90 and -90 <= north <= 90):
122 raise CliError("bbox latitude must be within [-90, 90] degrees")
123 if south > north:
124 raise CliError("bbox south must not exceed north")
125 crosses_antimeridian = bool(source.is_geographic and east < west)
126 if not crosses_antimeridian:
127 try:
128 area_of_interest = AreaOfInterest(west, south, east, north)
129 except (TypeError, ValueError) as exc:
130 raise CliError("bbox is not a valid area of interest") from exc
132 try:
133 group = TransformerGroup(
134 source,
135 target,
136 always_xy=True,
137 area_of_interest=area_of_interest,
138 allow_ballpark=args.allow_ballpark,
139 )
140 except (TypeError, ValueError) as exc:
141 raise CliError("PROJ could not construct transformation candidates") from exc
143 candidates = [_candidate_summary(item) for item in group.transformers[:5]]
144 source_info = crs_summary(source)
145 target_info = crs_summary(target)
146 linear_target = bool(
147 target.is_projected
148 and target.axis_info
149 and all(
150 "degree" not in str(axis.unit_name).casefold() for axis in target.axis_info[:2]
151 )
152 )
153 metric_ready = args.operation not in METRIC_OPERATIONS or linear_target
154 warnings: list[str] = []
155 if source.is_geographic:
156 warnings.append(
157 "Source coordinates are angular; GeoPandas segment transforms are vertex-wise, not geodesic."
158 )
159 if args.operation in METRIC_OPERATIONS and not linear_target:
160 warnings.append(
161 "Target CRS does not provide projected linear axes for the requested planar operation."
162 )
163 if crosses_antimeridian:
164 warnings.append(
165 "The bbox crosses the antimeridian; split/unwrap and densify as justified before to_crs."
166 )
167 if not group.best_available:
168 warnings.append(
169 "The best known operation is unavailable, commonly because a transformation grid is absent."
170 )
171 if not candidates:
172 warnings.append("No available transformation candidate satisfies this policy.")
173 if _dynamic(source) or _dynamic(target):
174 warnings.append(
175 "A dynamic CRS is involved; record and validate the coordinate epoch separately."
176 )
178 return {
179 "ok": bool(candidates) and metric_ready and bool(group.best_available),
180 "tool": TOOL,
181 "source_crs": source_info,
182 "target_crs": target_info,
183 "purpose": args.operation,
184 "coordinate_array_order": "x_y_via_always_xy",
185 "source_dynamic_crs": _dynamic(source),
186 "target_dynamic_crs": _dynamic(target),
187 "bbox": {
188 "provided": bbox_values is not None,
189 "values_emitted": False,
190 "used_as_area_of_interest": area_of_interest is not None,
191 "crosses_antimeridian": crosses_antimeridian,
192 },
193 "operation_policy": {
194 "allow_ballpark": bool(args.allow_ballpark),
195 "only_best_recommended_for_execution": True,
196 "proj_network_enabled": False,
197 "metric_operation_ready": metric_ready,
198 },
199 "candidate_count": len(group.transformers),
200 "candidate_summaries_first_five": candidates,
201 "best_available": bool(group.best_available),
202 "unavailable_operation_count": len(group.unavailable_operations),
203 "coordinates_transformed": False,
204 "network_accessed": False,
205 "coordinates_emitted": False,
206 "identifiers_emitted": False,
207 "stack": {"packages": package_versions(), "native": native_versions()},
208 "warnings": warnings,
209 }
212def main(argv: list[str] | None = None) -> int:
213 args = build_parser().parse_args(argv)
214 try:
215 report = plan(args)
216 emit_json(report)
217 return 0 if report["ok"] else 2
218 except Exception as exc: # noqa: BLE001 - errors are redacted at CLI boundary
219 return fail_json(TOOL, exc)
222if __name__ == "__main__":
223 sys.exit(main())