"""Precise CRS transformations: bundled IBGE grids and best-only transformations.
Converting coordinates between datums needs a datum shift. For Brazil's legacy
datums — SAD69, SAD69(96), Córrego Alegre 1961 and 1970-72 — the official shift
is IBGE's correction grid. pyproj wheels ship no grids, and without them PROJ
silently falls back to a formula that is metres off (SAD69 up to ~7.5 m) or, for
Córrego Alegre 1961, which has no formula, to a "ballpark" transformation that
ignores the datum shift altogether (~50 m off). Nothing signals the fallback, and
coordinates that fail to transform become ``inf``.
This module removes those silent paths:
* ``register_bundled_grids`` adds the four IBGE grids shipped with shapelib
(``shapelib/data/proj``) to PROJ's search path, so the official shift is used
offline;
* ``transform_geometries`` transforms with the best transformation PROJ knows
only (``only_best``), never a ballpark one, never a datum shift outside its
area of use (PROJ would otherwise apply, say, Mexico's NAD27 shift to a point in
São Paulo), and checks every coordinate; when that cannot be done it raises
``DatumTransformationError`` instead of degrading;
* ``datum_shift`` tells whether a CRS's datum needs a shift to reach WGS 84, and
how accurate that shift is, so data in a legacy datum can be flagged.
"""
from __future__ import annotations
import functools
import os
import re
import warnings
from dataclasses import dataclass
from pathlib import Path
from typing import Any
import numpy as np
import pyproj
import shapely
from pyproj import CRS, Transformer
from pyproj.aoi import AreaOfUse
from pyproj.exceptions import CRSError, ProjError
from pyproj.transformer import TransformerGroup
from shapelib.exceptions import DatumTransformationError
BUNDLED_GRIDS_DIR = Path(__file__).resolve().parents[1] / "data" / "proj"
WGS84 = CRS.from_epsg(4326)
# PROJ pipeline steps that move coordinates from one datum to another. A
# transformation without any of them only re-expresses the same datum (a
# projection, a unit or axis change, or a null shift between modern datums such
# as SIRGAS 2000 and WGS 84).
_DATUM_SHIFT_STEP = re.compile(r"\bproj=(?:hgridshift|gridshift|xyzgridshift|helmert|molodensky|deformation)\b")
[docs]
@dataclass(frozen=True)
class DatumShift:
"""A datum shift applied to reach WGS 84.
Attributes
----------
operation : str
The transformation used, as PROJ describes it.
accuracy : float | None
Its published accuracy in metres, or ``None`` when none is published
(e.g. ``TOWGS84`` parameters taken from a ``.prj``).
"""
operation: str
accuracy: float | None
[docs]
def register_bundled_grids() -> None:
"""
Add the IBGE grids bundled with shapelib to PROJ's search path.
Called when ``shapelib`` is imported, before any transformation. Idempotent.
Examples
--------
>>> register_bundled_grids()
>>> str(BUNDLED_GRIDS_DIR) in pyproj.datadir.get_data_dir()
True
"""
grids_dir = str(BUNDLED_GRIDS_DIR)
if grids_dir not in pyproj.datadir.get_data_dir().split(os.pathsep):
pyproj.datadir.append_data_dir(grids_dir)
[docs]
def datum_shift(crs: CRS, geometries: np.ndarray) -> DatumShift | None:
"""
The datum shift taking data in ``crs`` to WGS 84, or ``None`` if it needs none.
The transformation is validated where the data lies, as in
`transform_geometries`. PROJ may use different transformations across a region
(a grid covers one area and a formula another), so the shift is read at the
data's extreme vertices and in its middle, and the least accurate one found is
returned.
Parameters
----------
crs : CRS
The CRS of the data.
geometries : np.ndarray
The data's shapely geometries, in ``crs``.
Returns
-------
DatumShift | None
The shift, or ``None`` when reaching WGS 84 applies no datum shift
(WGS 84, SIRGAS 2000 and other modern datums).
Raises
------
DatumTransformationError
If the data cannot be transformed to WGS 84 with the best accuracy known.
Examples
--------
>>> datum_shift(CRS.from_epsg(29193), geometries) # SAD69 / UTM zone 23S
DatumShift(operation='... SAD69 to WGS 84 (15) ...', accuracy=2.0)
"""
if crs.is_exact_same(WGS84):
return None
transformer = _best_transformer(crs, WGS84)
shifts = []
for operation in _validate_where_the_data_lies(transformer, crs, WGS84, geometries):
if _DATUM_SHIFT_STEP.search(operation.definition):
accuracy = operation.accuracy if operation.accuracy >= 0 else None
shifts.append(DatumShift(operation.description, accuracy))
if not shifts:
return None
if any(shift.accuracy is None for shift in shifts):
return next(shift for shift in shifts if shift.accuracy is None)
return max(shifts, key=lambda shift: shift.accuracy or 0.0)
def _validate_where_the_data_lies(
transformer: Transformer, source: CRS, target: CRS, geometries: np.ndarray
) -> list[Transformer]:
"""
Check the transformation at the data's extreme vertices; return the operations used.
When no transformation between two datums is defined where a point lies, PROJ
still applies one defined elsewhere — e.g. Mexico's NAD27 shift to a point in
São Paulo — and the result looks valid while being meaningless. So at each
sample point, every datum shift PROJ used (each "Transformation" step; the
"Conversion" steps, such as a UTM projection, may legitimately run past their
nominal zone) must contain the point in its area of use.
The samples are the vertices at the data's west, east, south and north extremes
plus one in the middle: if all four extremes lie inside a shift's area of use
(a lon/lat rectangle), the data's whole bounding box does too.
"""
samples = _sample_points(geometries)
if not len(samples):
return []
to_lonlat = _to_geographic(source)
used = []
for x, y in samples:
try:
transformer.transform(x, y, errcheck=True)
except ProjError as error:
raise DatumTransformationError(_explain(source, target, error)) from error
operation = _last_used_operation(transformer)
lon, lat = to_lonlat.transform(x, y) if to_lonlat else (x, y)
for shift in _datum_shift_steps(operation):
area = shift.area_of_use
if area is not None and not _contains(area, lon, lat):
raise DatumTransformationError(
f"Cannot transform {_describe(source)} to {_describe(target)} where the data lies "
f"(lon {lon:.4f}, lat {lat:.4f}): the transformation between their datums, {shift.name}, is "
f"defined for {area.name.rstrip('.')}, not there, and applying it outside its area gives "
"meaningless coordinates. If the data really lies there, its CRS is probably declared wrong: "
"confirm it with whoever produced the data."
)
used.append(operation)
return used
def _datum_shift_steps(operation: Transformer) -> list[Any]:
"""The datum transformations within an operation, without its conversions (projections, axis swaps)."""
steps = tuple(operation.operations or ())
if steps:
return [step for step in steps if step.type_name == "Transformation"]
return [operation] if _DATUM_SHIFT_STEP.search(operation.definition) else []
def _contains(area: AreaOfUse, lon: float, lat: float) -> bool:
"""Whether an area of use contains a lon/lat point (areas may cross the antimeridian)."""
in_lon = area.west <= lon <= area.east if area.west <= area.east else lon >= area.west or lon <= area.east
return in_lon and area.south <= lat <= area.north
@functools.lru_cache(maxsize=128)
def _best_transformer(source: CRS, target: CRS) -> Transformer:
"""The transformer that only uses the best transformation and never a ballpark one."""
try:
return Transformer.from_crs(
_with_named_datum(source), _with_named_datum(target), always_xy=True, only_best=True, allow_ballpark=False
)
except ProjError as error:
raise DatumTransformationError(_explain(source, target, error)) from error
SIRGAS_2000 = CRS.from_epsg(4674)
def _named_by_epsg_code(crs: CRS) -> CRS | None:
"""The EPSG CRS ``crs`` is named after (``PROJCS["EPSG:31982",...]``) when its ellipsoid and projection agree.
QGIS writes the CRS name as its code but leaves the datum name empty, so PROJ cannot place the datum. The name is
only trusted when the EPSG CRS has the same ellipsoid, projection method and parameters; None otherwise.
"""
match = re.fullmatch(r"EPSG:(\d+)", crs.name)
if match is None or crs.ellipsoid is None:
return None
try:
named = CRS.from_epsg(int(match.group(1)))
except CRSError:
return None
if named.ellipsoid is None or named.is_projected != crs.is_projected:
return None
same_ellipsoid = named.ellipsoid.semi_major_metre == crs.ellipsoid.semi_major_metre and abs(
named.ellipsoid.inverse_flattening - crs.ellipsoid.inverse_flattening
) < 1e-6
if not same_ellipsoid:
return None
if named.is_projected:
# Parameters PROJ defaults when a WKT omits them (latitude of origin 0) read as 0.
parsed = {p.name: p.value for p in crs.coordinate_operation.params}
if named.coordinate_operation.method_name != crs.coordinate_operation.method_name or any(
abs(parsed.get(p.name, 0.0) - p.value) > 1e-9 for p in named.coordinate_operation.params
):
return None
return named
def _with_named_datum(crs: CRS) -> CRS:
"""``crs`` on SIRGAS 2000's datum when it is SIRGAS 2000 written with an unnamed datum.
Esri exports write SIRGAS 2000 as ``GEOGCS["GCS_SIRGAS_2000",DATUM["D_unnamed",SPHEROID["GRS_1980",...]]]``.
PROJ cannot place an unnamed datum, so the file would be refused like one that declares no datum at all. It is
not that case: the CRS names SIRGAS 2000 and carries its ellipsoid, GRS 1980, which no legacy Brazilian datum
uses (SAD69 is on GRS 1967, Córrego Alegre on International 1924), so reading it as SIRGAS 2000 is not a guess.
Every other CRS is returned unchanged.
"""
geodetic = crs.geodetic_crs
if geodetic is None or crs.datum is None or crs.ellipsoid is None:
return crs
names_sirgas_2000 = "sirgas" in geodetic.name.lower() and "2000" in geodetic.name
unnamed = crs.datum.name.lower().removeprefix("d_") == "unnamed"
grs_1980 = crs.ellipsoid.semi_major_metre == SIRGAS_2000.ellipsoid.semi_major_metre and abs(
crs.ellipsoid.inverse_flattening - SIRGAS_2000.ellipsoid.inverse_flattening
) < 1e-9
if not (names_sirgas_2000 and unnamed and grs_1980):
# GDAL reads a datum without a name as "unnamed", other readers as "" or "unknown".
nameless = crs.datum.name.lower().removeprefix("d_") in {"", "unnamed", "unknown"}
return (_named_by_epsg_code(crs) or crs) if nameless else crs
if crs.is_geographic:
return SIRGAS_2000
return pyproj.crs.ProjectedCRS(
conversion=crs.coordinate_operation, name=crs.name, cartesian_cs=crs.coordinate_system, geodetic_crs=SIRGAS_2000
)
@functools.lru_cache(maxsize=128)
def _to_geographic(crs: CRS) -> Transformer | None:
"""Conversion (no datum shift) from ``crs`` to lon/lat on its own datum; None if already lon/lat."""
if crs.is_geographic or crs.geodetic_crs is None:
return None
return Transformer.from_crs(crs, crs.geodetic_crs, always_xy=True)
def _last_used_operation(transformer: Transformer) -> Transformer:
"""The operation the last ``transform`` call used (the transformer itself if it has only one)."""
try:
return transformer.get_last_used_operation()
except ProjError:
return transformer
def _sample_points(geometries: np.ndarray) -> np.ndarray:
"""Up to five vertices: the data's west, east, south and north extremes and one in the middle."""
vertices: np.ndarray = shapely.get_coordinates(geometries)
if not len(vertices):
return vertices
picks = {
int(np.argmin(vertices[:, 0])),
int(np.argmax(vertices[:, 0])),
int(np.argmin(vertices[:, 1])),
int(np.argmax(vertices[:, 1])),
len(vertices) // 2,
}
return vertices[sorted(picks)]
def _declares_no_datum(crs: CRS) -> bool:
"""Whether ``crs`` names only an ellipsoid, without a datum (e.g. Esri's ``D_GRS_1980``, EPSG:4019)."""
datum = crs.datum.name.lower().removeprefix("d_") if crs.datum else ""
return datum.startswith(("not specified", "unknown"))
def _describe(crs: CRS) -> str:
code = crs.to_epsg()
return f"{crs.name} (EPSG:{code})" if code else crs.name
def _explain(source: CRS, target: CRS, error: ProjError) -> str:
"""Why ``source`` cannot be transformed to ``target`` precisely, in plain terms."""
with warnings.catch_warnings():
warnings.simplefilter("ignore") # TransformerGroup warns about the very grids listed below
operations = TransformerGroup(source, target, always_xy=True).unavailable_operations
missing = sorted({grid.short_name for operation in operations for grid in operation.grids if not grid.available})
prefix = f"Cannot transform {_describe(source)} to {_describe(target)} with the best accuracy known: "
if missing:
reason = f"the best transformation needs the correction grid(s) {', '.join(missing)}, which are not installed."
elif undeclared := next((crs for crs in (source, target) if _declares_no_datum(crs)), None):
ellipsoid = undeclared.ellipsoid.name if undeclared.ellipsoid else "an"
reason = (
f"{_describe(undeclared)} declares no datum, only the {ellipsoid} ellipsoid, so there is no defined way "
"to place it relative to another datum. Confirm the datum with whoever produced the data and declare it "
"in the file (e.g. a .prj for SIRGAS 2000). shapelib does not guess it: if the data were in a legacy "
"datum such as SAD69, a guess would put it tens of metres off."
)
else:
reason = (
"there is no known transformation between their datums — the file may name a datum PROJ does not "
"recognize, or carry only an ellipsoid. Falling back to a 'ballpark' transformation would ignore the "
"datum shift, which can be tens of metres."
)
return f"{prefix}{reason} shapelib refuses less accurate fallbacks. (PROJ: {error})"
# On import, before any transformation can run: shapelib's package __init__ imports
# this module (and every BaseShapes method that transforms depends on it).
register_bundled_grids()