"""Resolving and applying a coordinate transformation, with the rules enforced.
Four things are checked here that PROJ will not check for you, because PROJ's
job is to compute and this package's job is to refuse to compute something
untrustworthy:
1. If an operation was requested, the transformer that gets built is inspected
to confirm it really contains that operation. PROJ will otherwise build a
perfectly functional transformer using a different one.
2. A ballpark path is refused outright rather than returned as an approximate
number with no usable accuracy statement.
3. A grid the operation depends on but that is not installed is an error, named
specifically, rather than a quiet fall back to a grid-free operation.
4. A time-dependent operation without a coordinate epoch is an error. A dynamic
reference frame alone does not require an epoch if the operation ignores time.
Coordinate **values** are in ``xy`` order, in and out: easting or longitude
first wherever a CRS has an easting and a northing, whatever order it declares
them in. A CRS with no such pair (Krovak's southing and westing, a geocentric
X/Y/Z, a plant grid declaring north and west) keeps its declared order, as
PROJ's ``always_xy`` keeps it; ``value_axis_order`` states the order axis by
axis. The CRSs' declared axis order is reported separately and is not changed
by this; see :mod:`geodetic_engine.geodesy.crs`. PROJ does not honour this
contract at an engineering CRS on its own; see the workaround block above
``_correct_engineering_axes``.
"""
from __future__ import annotations
import json
import math
import os
import warnings
from collections.abc import Iterable, Iterator, Sequence
from contextlib import contextmanager
from dataclasses import dataclass, replace
from functools import lru_cache
from pathlib import Path
from typing import Any, cast
from pyproj import CRS, Transformer, datadir
from pyproj.crs import CoordinateOperation
from pyproj.enums import TransformDirection
from pyproj.exceptions import CRSError, ProjError
from pyproj.network import is_network_enabled
from pyproj.transformer import TransformerGroup
from geodetic_engine.geodesy.crs import AxisSpec, CoordinateReferenceSystem
from geodetic_engine.geodesy.database import (
DatabaseIdentity,
database_fingerprints,
database_identity,
)
from geodetic_engine.geodesy.errors import (
AmbiguousOperationError,
BallparkTransformationError,
CoordinateOutOfRangeError,
MissingCoordinateEpochError,
MissingGridError,
OperationNotAvailableError,
TransformationFailedError,
)
from geodetic_engine.geodesy.operation import (
_VISUALIZATION_SUFFIX,
AppliedOperation,
AreaOfUse,
GridUsage,
OperationCandidate,
OperationReference,
OperationRequest,
OperationRoute,
OperationStep,
StatedOperation,
base_authority,
datum_operation_count,
fully_requested,
grid_usages,
has_inverted_step,
is_ballpark,
operation_names,
operation_references,
parse_operations,
requires_epoch,
)
from geodetic_engine.geodesy.result import Coordinates, TransformationResult
_VERTICAL_DIRECTIONS = frozenset({"up", "down"})
# Axis directions that carry latitude and longitude in a geographic CRS.
_NORTHINGS = frozenset({"north", "south"})
_EASTINGS = frozenset({"east", "west"})
# PROJ transforms at most x, y, z: a fourth spatial component has no meaning to
# it, so one extra value beyond what a CRS declares is tolerated (a height
# alongside a 2D horizontal CRS, carried through unchanged) but no more than
# that.
_MAX_COORDINATE_VALUES = 3
# PROJ's flag for running a step backwards. A bare token, never a parameter
# with a value, so it can be added or removed by name.
_INVERSE_FLAG = "inv"
# How PROJ words a grid it cannot open while building a pipeline step. The
# message names the step, never the file.
_MISSING_FILE_MARKER = "File not found or invalid"
# Methods that restate axes rather than move coordinates. PROJ inserts these
# when it normalises axis order, and they are not the method a caller means.
_BOOKKEEPING_METHODS = frozenset(
{
"axis order reversal (2d)",
"axis order reversal (geographic3d horizontal)",
"change of vertical unit",
"geographic3d to geographic2d conversion",
}
)
# Vertical operation methods whose shift is the same everywhere. Their result
# does not read the horizontal position, unlike a grid interpolation or
# "Vertical Offset and Slope", so a position supplied alongside the height is
# carried rather than interpreted as a latitude and longitude.
_POSITION_FREE_VERTICAL_METHODS = frozenset(
{
"vertical offset",
"change of vertical unit",
"height depth reversal",
}
)
@dataclass(frozen=True, slots=True)
class _Pipeline:
"""One or more PROJ transformers applied in sequence.
A single operation is a pipeline of one. Chaining exists for the case where
the requested operation is defined between different CRSs than the caller's
pair, for example a datum shift published between geographic CRSs that a
caller wants applied between two projected ones.
"""
steps: tuple[tuple[Transformer, TransformDirection], ...]
core: Transformer
route: OperationRoute
identified_by: tuple[OperationRequest, ...] = ()
"""Operations that name this pipeline when the caller did not name one.
Set for a bound CRS, whose definition states the operation itself: one
entry per bound end. Kept apart from the caller's request so that a result
reports what was applied without claiming it was asked for.
"""
corrected: tuple[tuple[Transformer, TransformDirection], ...] | None = None
"""The steps actually run, when a PROJ workaround had to rewrite one.
``steps`` keeps the transformers PROJ built, which carry the provenance
(definition, accuracy, grids). This holds the same sequence with the
rewritten core in place of the original, and is what executes and what is
reported as the pipeline. None when nothing was rewritten. See
``_correct_engineering_axes``.
"""
@property
def executed(self) -> tuple[tuple[Transformer, TransformDirection], ...]:
"""The steps that run, corrected where a PROJ workaround applies."""
return self.steps if self.corrected is None else self.corrected
def run(
self, columns: Sequence[Sequence[float]], epoch: float | None
) -> list[list[float]]:
"""Apply every step in order, carrying all coordinate components through."""
values: list[list[float]] = [list(column) for column in columns]
for transformer, direction in self.executed:
values = _apply(transformer, values, epoch, direction)
return values
@property
def core_direction(self) -> TransformDirection:
"""Direction in which the raw core definition is executed."""
return next(
direction
for transformer, direction in self.steps
if transformer is self.core
)
@property
def definition(self) -> dict[str, Any]:
"""PROJJSON of the resolved operation."""
try:
return dict(self.core.to_json_dict())
except (TypeError, ProjError):
return {}
@property
def accuracy(self) -> float | None:
"""Stated accuracy in metres, or None when PROJ reports none."""
value = self.core.accuracy
return None if value is None or value < 0 else float(value)
@property
def text(self) -> str | None:
"""The whole chain as one PROJ pipeline definition, ready to be rebuilt.
Every step is flattened into a single ``proj=pipeline``, and a step
that was applied backwards is written out with PROJ's ``inv`` flag and
its own steps reversed, so what comes back is what ran rather than a
listing of the parts. Rebuild it with
:meth:`pyproj.Transformer.from_pipeline`.
It reads and writes PROJ's own components in PROJ's own order, which
at a vertical end is not the caller's ``xy`` order; see
``_order_horizontal_for_pipeline``.
Returns:
The definition, or None when the steps cannot be written as one
pipeline without changing what they do. Provenance that cannot be
replayed is worse than none, so nothing is reported rather than a
string that only looks executable.
"""
return _compose_pipeline(self.executed)
@property
def reads_declared_horizontal(self) -> bool:
"""Whether the first step consumes its horizontal pair northing-first.
``always_xy`` normalises the axis order of each end that has horizontal
axes. A vertical CRS has none, so the ``axisswap`` PROJ inserts for the
operation's own geographic order survives normalisation there and the
accompanying position is read latitude-first. Reading it off the
pipeline is the only statement of that order there is: the vertical
CRS does not declare one, and the operation's own geographic end is not
necessarily either CRS the caller named.
PROJ workaround, not a permanent design: see the note above
``_entry_step`` below for the upstream bug and how to retire this once
it is fixed.
"""
if not self.executed:
return False
transformer, direction = self.executed[0]
return _swaps_horizontal(_entry_step(transformer.definition, direction))
def _columns(
crs: CoordinateReferenceSystem, points: Iterable[Iterable[float]] | Iterable[float]
) -> tuple[tuple[float, ...], ...]:
"""Reshape points into one tuple of values per axis, the shape PROJ wants.
Args:
crs: The CRS the points are expressed in, whose declared axis count
bounds how many values each point may carry.
points: A single point's values, given flat -- ``(lon, lat)`` -- or a
batch: an iterable of coordinate iterables, each holding one
point's values in ``xy`` order. A 2D numpy array of shape
``(n_points, n_axes)`` works, one row per point. A row may carry
one value more than ``crs`` declares -- a height alongside a 2D
horizontal CRS -- which is carried through unchanged rather than
consumed, matching how :meth:`pyproj.Transformer.transform` accepts
an optional ``zz`` regardless of what the CRS pair declares.
Returns:
One tuple of values per axis, so a whole batch crosses into PROJ in a
single call instead of once per point.
Raises:
ValueError: If points disagree on how many values they carry, or that
count is not ``crs``'s declared dimension, or one more.
"""
# Materialised up front: points may be a one-shot iterable or a numpy array
# (not a Sequence), and each axis is read once below. Whether an element is
# a point or one value of a single flat point is only known at runtime.
materialized: list[Any] = list(points)
if materialized and not isinstance(materialized[0], Iterable):
# A lone point given flat, e.g. (lon, lat), rather than [(lon, lat)].
# Unambiguous whenever a point has more than one value: only a single
# flat point looks like a list of bare numbers rather than of rows.
materialized = [materialized]
rows = [tuple(float(value) for value in point) for point in materialized]
widths = {len(row) for row in rows}
if len(widths) > 1:
raise ValueError(f"points have differing numbers of values: {sorted(widths)}")
width = widths.pop() if widths else crs.dimension
_require_width(crs, width)
return tuple(tuple(row[axis] for row in rows) for axis in range(width))
def _columns_from_axes(
crs: CoordinateReferenceSystem,
x: Iterable[float] | float,
y: Iterable[float] | float,
z: Iterable[float] | float | None,
) -> tuple[tuple[float, ...], ...]:
"""Reshape separate per-axis values into columns, broadcasting a lone scalar.
Mirrors :meth:`pyproj.Transformer.transform`'s ``xx, yy, zz`` convention:
each axis is either the whole batch's values, or a lone scalar applied to
every point -- a fixed height for many horizontal points, for example,
given once rather than repeated per point.
Args:
crs: The CRS the points are expressed in, whose declared axis count
bounds how many axes may be given.
x: First axis's values: a scalar for one point, or a sequence (a list
or a 1D numpy array) for a batch.
y: Second axis's values, in the same shape as ``x``.
z: Third axis's values, if any, in the same shape as ``x`` and ``y``.
Returns:
One tuple of values per axis, so a whole batch crosses into PROJ in a
single call instead of once per point.
Raises:
ValueError: If the sequence axes disagree on how many points they
hold, or the number of axes given is not ``crs``'s declared
dimension, or one more.
"""
axes = (x, y) if z is None else (x, y, z)
_require_width(crs, len(axes))
values: list[tuple[float, ...] | float] = [
tuple(float(value) for value in axis)
if isinstance(axis, Iterable)
else float(axis)
for axis in axes
]
lengths = {len(column) for column in values if isinstance(column, tuple)}
if len(lengths) > 1:
raise ValueError(f"axes have differing batch sizes: {sorted(lengths)}")
count = lengths.pop() if lengths else 1
return tuple(
column if isinstance(column, tuple) else (column,) * count for column in values
)
def _require_width(crs: CoordinateReferenceSystem, width: int) -> None:
"""Check how many values a point carries against what ``crs`` allows.
Raises:
ValueError: If ``width`` is not ``crs``'s declared dimension or one
more (a height alongside a 2D horizontal CRS, carried through
unchanged); or, for a vertical CRS, anything other than a full
``(lon, lat, h)`` triple.
"""
if crs.dimension == 1 and _is_vertical(crs):
# A height shifted by a geoid or vertical datum grid is only defined
# where the grid is read, so the horizontal position that locates it
# has to be given even though the CRS declares no axis for it. The
# values go to PROJ in x, y, z order, so a lone height would be read
# as a longitude rather than as a height.
if width != _MAX_COORDINATE_VALUES:
raise ValueError(
f"{width} values were given per point but {crs!r} is a vertical "
f"CRS; exactly {_MAX_COORDINATE_VALUES} are accepted, "
"(lon, lat, h), because the height is only defined at the "
"horizontal position its grid is read at"
)
return
allowed = {crs.dimension}
if crs.dimension < _MAX_COORDINATE_VALUES:
allowed.add(crs.dimension + 1)
if width not in allowed:
raise ValueError(
f"{width} values were given per point but {crs!r} declares "
f"{crs.dimension} axes; {sorted(allowed)} values are accepted "
"(the extra being a height PROJ carries through unchanged)"
)
[docs]
def available_operations(
source_crs: Any,
target_crs: Any,
*,
authority: str | None = "any",
accuracy: float | None = None,
allow_superseded: bool = True,
allow_ballpark: bool = True,
) -> tuple[OperationCandidate, ...]:
"""List every coordinate operation PROJ offers between two CRSs.
This package's equivalent of inspecting a
:class:`pyproj.transformer.TransformerGroup` directly: every candidate is
described, including a ballpark fallback or one whose grid is not
installed, so an ``operation=`` argument for :class:`Transformation` can be
chosen with full information instead of by trial and error. Nothing here
is applied to coordinates or checked against a request.
A deprecated EPSG operation is never among the candidates: PROJ's own
operation search excludes deprecated operations unconditionally, with no
option to include them, so there is no ``allow_deprecated`` filter here to
match -- one would silently do nothing.
Args:
source_crs: CRS the input coordinates would be in.
target_crs: CRS to produce coordinates in.
authority: Restrict candidates to those published by this authority,
for example ``"EPSG"``. ``"any"`` searches every authority without
the preference PROJ otherwise gives the source/target CRS's own
authority. Defaults to "any"; pass None to use PROJ's preference.
accuracy: Discard candidates stated as less accurate than this, in
metres. Omitted by default, so every accuracy is considered.
allow_superseded: Whether to include an operation EPSG has marked as
superseded by a newer one. True by default, since a superseded
operation is still valid, just no longer preferred.
allow_ballpark: Whether to include a ballpark approximation among the
candidates. True by default, so its presence and its lack of a
usable accuracy are visible here rather than only discovered when
:class:`Transformation` refuses it.
Returns:
One candidate per operation PROJ offers, in the order PROJ ranks them,
which weighs area of use and other criteria alongside accuracy: for
ED50 to WGS 84 the first candidate is EPSG:1133 (10 m) ahead of
EPSG:1612 (1 m). Filter on
:attr:`~geodetic_engine.geodesy.operation.OperationCandidate.accuracy`
or pass ``accuracy=`` rather than taking the first entry as the best.
Pass any entry's
:attr:`~geodetic_engine.geodesy.operation.OperationCandidate.authority_code`
as ``Transformation``'s ``operation=`` argument.
Example:
>>> candidates = available_operations("EPSG:4230", "EPSG:4326")
>>> candidates[0].authority_code
'EPSG:1133'
>>> most_accurate = available_operations(
... "EPSG:4230", "EPSG:4326", accuracy=1.0, allow_ballpark=False
... )
>>> all(c.accuracy is not None and c.accuracy <= 1.0 for c in most_accurate)
True
"""
source = CoordinateReferenceSystem.from_user_input(source_crs)
target = CoordinateReferenceSystem.from_user_input(target_crs)
with _proj_construction(source, target):
group = TransformerGroup(
source.crs,
target.crs,
always_xy=True,
authority=authority,
accuracy=accuracy,
allow_ballpark=allow_ballpark,
allow_superseded=allow_superseded,
crs_extent_use="none",
grid_check="none",
)
return tuple(_describe_candidate(transformer) for transformer in group.transformers)
def _describe_candidate(transformer: Transformer) -> OperationCandidate:
"""Describe one candidate transformer without applying or requesting it."""
definition = transformer.to_json_dict()
operations = tuple(transformer.operations or ())
if not operations:
try:
operations = (CoordinateOperation.from_json_dict(definition),)
except CRSError:
operations = ()
steps = _substantive_steps(operations)
# The whole candidate's own id, which a registered concatenated operation
# such as EPSG:8047 carries even though it applies two Helmerts. Absent
# only when PROJ assembled the chain itself, and then no code may be
# reported: borrowing a step's would understate the rest of the pipeline.
identifier = _identifier(definition) or (
(steps[0].auth_name, steps[0].code)
if len(steps) == 1 and steps[0].auth_name and steps[0].code
else None
)
area = transformer.area_of_use
ballpark = is_ballpark(definition)
grids = grid_usages(operations)
return OperationCandidate(
auth_name=None if identifier is None else identifier[0],
code=None if identifier is None else identifier[1],
name=_candidate_name(definition, steps, transformer),
method_name=steps[0].method_name if len(steps) == 1 else None,
accuracy=transformer.accuracy if transformer.accuracy >= 0 else None,
area_of_use=(
None
if area is None
else AreaOfUse(
west=area.west,
south=area.south,
east=area.east,
north=area.north,
name=area.name,
)
),
ballpark=ballpark,
requires_epoch=requires_epoch(definition, operations),
grids=grids,
steps=steps,
projjson=json.dumps(definition),
usable=not ballpark and all(grid.available for grid in grids),
)
def _candidate_name(
definition: dict[str, Any],
steps: tuple[OperationStep, ...],
transformer: Transformer,
) -> str:
"""The candidate's name, without PROJ's axis-normalisation annotation.
Normalising axis order makes PROJ append "(with axis order normalized for
visualization)" to the name of the operation it wraps, which is wording
the caller never asked to have on a result. A single step's own name is
already clean; otherwise the annotation is stripped from the whole.
"""
if len(steps) == 1:
return steps[0].name
name = str(definition.get("name") or transformer.description)
return name.removesuffix(_VISUALIZATION_SUFFIX) or " + ".join(
step.name for step in steps
)
def _substantive_steps(
operations: tuple[CoordinateOperation, ...],
) -> tuple[OperationStep, ...]:
"""Every operation that does geodetic work, not the axis bookkeeping."""
return tuple(
OperationStep(
auth_name=None if identifier is None else identifier[0],
code=None if identifier is None else identifier[1],
name=str(node.get("name") or operation.name),
method_name=(str(operation.method_name) if operation.method_name else None),
)
for operation in operations
if not _is_bookkeeping(operation)
for node in (operation.to_json_dict(),)
for identifier in (_identifier(node),)
)
def _is_bookkeeping(operation: CoordinateOperation) -> bool:
"""Whether a step restates axes or units rather than moving coordinates."""
name = operation.method_name
return bool(name) and name.strip().lower() in _BOOKKEEPING_METHODS
def _cache_key(crs: Any) -> str:
"""Render a CRS input as a hashable definition string."""
return CoordinateReferenceSystem.from_user_input(crs).definition
@lru_cache(maxsize=128)
def _cached_transformation(
source: str,
target: str,
operation: (
str
| int
| OperationReference
| tuple[str | int | OperationReference, ...]
| None
),
allow_any_operation: bool,
identity: DatabaseIdentity,
) -> Transformation:
"""Resolve and cache a transformation by its textual inputs."""
return Transformation(
source, target, operation, allow_any_operation=allow_any_operation
)
@contextmanager
def _proj_construction(
source: CoordinateReferenceSystem, target: CoordinateReferenceSystem
) -> Iterator[None]:
"""Report PROJ's refusal to build anything for a CRS pair as a package error.
PROJ raises for pairs it has no notion of a path between at all, such as a
vertical CRS to a geographic one. That is a real failure and is not hidden,
but it reaches the caller as :class:`OperationNotAvailableError` naming both
CRSs rather than as a bare ``ProjError`` from inside pyproj.
pyproj's warning about the best ranked candidate is dropped here. Candidates
are only being enumerated, most are discarded, and this package never takes
PROJ's ranking silently: it applies the operation the caller named, or
refuses the pair as ambiguous. What the chosen operation needs is reported
by :attr:`Transformation.grids`, and a grid that is genuinely missing raises
:class:`MissingGridError` naming it.
"""
try:
with warnings.catch_warnings():
warnings.filterwarnings(
"ignore",
message="Best transformation is not available",
category=UserWarning,
)
yield
except (ProjError, CRSError, IndexError) as error:
raise OperationNotAvailableError(
f"PROJ cannot build a transformation from {_label(source)} to "
f"{_label(target)}: {error}"
) from error
def _resolve(
source: CoordinateReferenceSystem,
target: CoordinateReferenceSystem,
requests: tuple[OperationRequest, ...],
*,
allow_any_operation: bool,
) -> _Pipeline:
"""Build the pipeline that will be applied, recording how it was found."""
if not requests:
return _resolve_without_request(
source, target, allow_any_operation=allow_any_operation
)
if any(request.definition is not None for request in requests):
if len(requests) > 1:
raise OperationNotAvailableError(
f"{', '.join(str(r) for r in requests)} were requested for "
f"{_label(source)} to {_label(target)}, but an operation stated "
"outright is applied on its own; it cannot be combined with "
"others, which are selected from what PROJ offers for the pair"
)
# Never offered to the candidate search: what PROJ publishes under the
# same code is a different object from the one stated here.
return _from_operation(source, target, requests[0])
found = _from_transformer_group(source, target, requests)
if found is not None:
return _Pipeline(
steps=((found, TransformDirection.FORWARD),),
core=found,
route=OperationRoute.TRANSFORMER_GROUP,
)
if len(requests) == 1:
return _from_operation(source, target, requests[0])
raise OperationNotAvailableError(
f"{', '.join(str(r) for r in requests)} were requested for "
f"{_label(source)} to {_label(target)}, but no candidate PROJ offers "
"for this CRS pair applies all of them together; naming more than one "
"operation is only supported among PROJ's own candidates, not chained "
"by hand"
)
def _resolve_without_request(
source: CoordinateReferenceSystem,
target: CoordinateReferenceSystem,
*,
allow_any_operation: bool,
) -> _Pipeline:
"""Resolve conversions or bound operations, refusing unrequested datum shifts.
WKT1 and ensemble representations can give equivalent frames different
names. An empty request set admits a candidate only when every step is
a conversion (or explicitly authorized by a bound CRS).
"""
if _datum_names(source.crs) != _datum_names(target.crs):
bound = _from_bound_crs(source, target)
if bound is not None:
return bound
conversion = _from_transformer_group(source, target, ())
if conversion is not None:
return _Pipeline(
steps=((conversion, TransformDirection.FORWARD),),
core=conversion,
route=OperationRoute.TRANSFORMER_GROUP,
)
raise AmbiguousOperationError(
f"{_label(source)} to {_label(target)} involves a datum change; "
"every datum operation must be named explicitly. "
"allow_any_operation no longer "
"permits automatic selection or ballpark results"
)
with _proj_construction(source, target):
transformer = Transformer.from_crs(
source.crs, target.crs, always_xy=True, allow_ballpark=False
)
return _Pipeline(
steps=((transformer, TransformDirection.FORWARD),),
core=transformer,
route=OperationRoute.PROJ_DEFAULT,
)
def _from_bound_crs(
source: CoordinateReferenceSystem, target: CoordinateReferenceSystem
) -> _Pipeline | None:
"""Use the transformation a bound CRS states as its own definition.
A bound CRS names exactly one transformation to its hub, so there is no
choice left for PROJ to make and nothing for the caller to disambiguate.
That is early binding, and it is the one datum change this package will
apply without being told which operation to use: the operation was declared
by whoever defined the CRS, not guessed here.
The operation is read out of the bound CRS and then resolved through the
same transformer group as a named one, rather than letting the bound CRS
build the transformer by itself. A bound CRS with a projected base needs
the map projection applied around the datum shift, and going through the
group is what supplies those steps and keeps the applied operation
identifiable.
Returns:
The pipeline, or None if neither CRS is bound, in which case the datum
change really is ambiguous.
"""
requests = tuple(
request
for request in (_bound_operation(source), _bound_operation(target))
if request is not None
)
if not requests:
return None
found = _from_transformer_group(source, target, requests)
if found is None:
found = _bound_transformer(source, target, requests)
if found is None:
return None
return _Pipeline(
steps=((found, TransformDirection.FORWARD),),
core=found,
route=OperationRoute.BOUND,
identified_by=requests,
)
def _bound_transformer(
source: CoordinateReferenceSystem,
target: CoordinateReferenceSystem,
requests: tuple[OperationRequest, ...],
) -> Transformer | None:
"""Chain the bound CRSs' own transformations without the candidate search.
``TransformerGroup`` withholds a candidate PROJ reports as not
instantiable, and it judges that on the grid name the authority published
rather than on the one it would actually read: the NADCON pair binding
EPSG:1188 to EPSG:15851 is withheld over ``conus.las`` while PROJ builds it
perfectly well from the installed ``us_noaa_conus.tif``. Early binding is
then left with nothing to apply and the caller is told the datum change is
ambiguous, when both CRSs had in fact declared their operation.
Building the pair outright goes through the same early binding without that
filter. Ballpark remains refused, and the result is accepted only once
every declared operation is confirmed present in what PROJ built, so this
cannot quietly substitute a different one.
"""
with _proj_construction(source, target):
transformer = Transformer.from_crs(
source.crs, target.crs, always_xy=True, allow_ballpark=False
)
definition = transformer.to_json_dict()
if all(
request.is_satisfied_by(definition) for request in requests
) and fully_requested(definition, requests):
return transformer
return None
def _bound_operation(
crs: CoordinateReferenceSystem,
) -> OperationRequest | None:
"""The operation a bound CRS embeds, by authority code or else by name.
A collapsed concatenated operation carries no identifier of its own, having
been synthesised rather than published, so the name is the only handle on
it. See :mod:`geodetic_engine.geodesy.utils.helmert`.
"""
if not crs.crs.is_bound:
return None
node = crs.crs.to_json_dict().get("transformation")
if not isinstance(node, dict):
return None
identifier = node.get("id")
if isinstance(identifier, dict):
authority, code = identifier.get("authority"), identifier.get("code")
if authority is not None and code is not None:
return OperationRequest.parse(f"{base_authority(str(authority))}:{code}")
name = node.get("name")
return OperationRequest.parse(str(name)) if name else None
def _from_transformer_group(
source: CoordinateReferenceSystem,
target: CoordinateReferenceSystem,
requests: tuple[OperationRequest, ...],
) -> Transformer | None:
"""Find a candidate PROJ offers that satisfies every one of the requests.
Extent filtering and grid filtering are both turned off, so that an
operation is not hidden merely because its grid is missing. A missing grid
is then reported as a missing grid rather than as a missing operation.
"""
with _proj_construction(source, target):
group = TransformerGroup(
source.crs,
target.crs,
always_xy=True,
allow_ballpark=False,
allow_superseded=True,
crs_extent_use="none",
grid_check="none",
)
for transformer in group.transformers:
definition = transformer.to_json_dict()
authorized = (
*requests,
*(
request
for crs in (source, target)
if (request := _bound_operation(crs)) is not None
),
)
if all(
request.is_satisfied_by(definition) for request in requests
) and fully_requested(definition, authorized):
return transformer
return None
def _grids_of(request: OperationRequest) -> tuple[GridUsage, ...]:
"""Grids the requested operation declares, read from the registry.
The registry answers even when the operation cannot be built, which is
what makes it usable for explaining that failure.
Args:
request: The operation the caller asked for.
Returns:
Its grids, or empty when the request names no authority code or the
registry does not know it.
"""
if request.auth_name is None or request.code is None:
return ()
try:
operation = CoordinateOperation.from_authority(request.auth_name, request.code)
except (CRSError, ProjError):
return ()
return grid_usages((operation,))
def _from_operation(
source: CoordinateReferenceSystem,
target: CoordinateReferenceSystem,
request: OperationRequest,
) -> _Pipeline:
"""Build the named operation itself, wrapping it in same-datum conversions.
Reached when the operation is not among the candidates for this CRS pair,
which typically means it is published between geographic CRSs while the
caller is working in projected ones, or when the operation was stated
outright and so was never a candidate to begin with.
"""
if request.definition is not None:
core = _stated_transformer(request, source, target)
else:
try:
core = Transformer.from_pipeline(
request.urn or request.text, always_xy=True
)
except (ProjError, CRSError) as error:
# PROJ reports an absent grid here as a malformed pipeline step,
# which says nothing about the grid. Its availability is only
# consulted now that building has already failed: PROJ resolves
# legacy grid names through proj.db's alternatives, so an operation
# whose grid reads as unavailable often still transforms, and
# checking earlier would refuse work that succeeds.
_require_grids(_grids_of(request), source, target)
raise OperationNotAvailableError(
f"{request} could not be built as a coordinate operation, and is "
f"not among the operations PROJ offers for {_label(source)} to "
f"{_label(target)}: {error}"
) from error
ends = _operation_ends(core)
if ends is None:
raise OperationNotAvailableError(
f"{request} does not declare the CRSs it operates between, so it "
f"cannot be chained into {_label(source)} to {_label(target)}"
)
op_source, op_target = ends
if _datum_names(source.crs) & _datum_names(op_source):
direction = TransformDirection.FORWARD
entry, exit_ = op_source, op_target
elif _datum_names(source.crs) & _datum_names(op_target):
direction = TransformDirection.INVERSE
entry, exit_ = op_target, op_source
else:
raise OperationNotAvailableError(
f"{request} operates between {op_source.name!r} and "
f"{op_target.name!r}, neither of which shares a datum with "
f"{_label(source)}; it cannot be applied here"
)
steps: list[tuple[Transformer, TransformDirection]] = []
if _datum_names(source.crs) != _datum_names(entry) or source.crs != entry:
steps.append(
(_conversion(source.crs, entry, request), TransformDirection.FORWARD)
)
steps.append((core, direction))
if _datum_names(target.crs) != _datum_names(exit_) or target.crs != exit_:
steps.append(
(_conversion(exit_, target.crs, request), TransformDirection.FORWARD)
)
return _Pipeline(steps=tuple(steps), core=core, route=OperationRoute.CHAINED)
def _stated_transformer(
request: OperationRequest,
source: CoordinateReferenceSystem,
target: CoordinateReferenceSystem,
) -> Transformer:
"""Run the operation a caller stated, rather than one PROJ looked up.
The definition is handed over whole, so a chain stays a chain: nothing has
to be collapsed into a single step the way a bound CRS would require.
Raises:
MissingGridError: If PROJ cannot build the operation because a grid it
reads is not installed. PROJ reports that as a malformed pipeline
step, which would otherwise send the caller looking for another
operation instead of for the grid.
OperationNotAvailableError: If PROJ cannot build it for any other
reason.
"""
definition = request.definition
if definition is None:
raise ValueError(f"{request} states no operation to run")
if has_inverted_step(definition.to_json_dict()):
# PROJJSON writes an inverted step with its forward parameters.
raise OperationNotAvailableError(
f"{request} states an operation with a step applied inverted, which "
"PROJ would run forwards once handed over; state that step in the "
"direction it is applied"
)
try:
return Transformer.from_pipeline(definition.to_json(), always_xy=True)
except (ProjError, CRSError) as error:
if _MISSING_FILE_MARKER in str(error):
_require_grids(grid_usages((definition,)), source, target)
raise OperationNotAvailableError(
f"{request} states an operation PROJ cannot run: {error}"
) from error
def _conversion(source: CRS, target: CRS, request: OperationRequest) -> Transformer:
"""Build a step that changes representation without changing datum.
Guards the chain: if this step would move between datums it would apply a
datum shift the caller never asked for, on top of the one they did.
"""
if _datum_names(source) != _datum_names(target):
raise OperationNotAvailableError(
f"applying {request} between {source.name!r} and {target.name!r} "
"would require an additional, unrequested datum change; name the "
"full operation instead"
)
return Transformer.from_crs(source, target, always_xy=True, allow_ballpark=False)
def _operation_ends(transformer: Transformer) -> tuple[CRS, CRS] | None:
"""The CRSs an operation is defined between, as full CRS objects."""
source = transformer.source_crs
target = transformer.target_crs
if source is None or target is None:
return None
return CRS.from_wkt(source.to_wkt()), CRS.from_wkt(target.to_wkt())
def _datum_names(crs: CRS) -> frozenset[str]:
"""Names of every datum the CRS is built on, including compound components.
WKT1 cannot represent ensembles. A plain geodetic datum is normalized to
a registered ensemble only when its geographic CRS is equivalent to that
registry definition. The registry lookup only proposes a candidate; its
confidence score does not establish equivalence. A name suffix is not evidence
that two frames are interchangeable.
"""
parts = crs.sub_crs_list or [crs]
names = set()
for part in parts:
datum = part.datum
if datum is not None:
name = str(datum.name)
geographic = part.geodetic_crs
if (
datum.type_name == "Geodetic Reference Frame"
and geographic is not None
and geographic.is_geographic
):
horizontal = geographic.to_2d()
if identified := horizontal.to_authority(
auth_name="EPSG", min_confidence=0
):
registered = CRS.from_authority(*identified).to_2d()
ensemble = registered.datum
if (
ensemble is not None
and ensemble.type_name == "Datum Ensemble"
and horizontal.equals(registered, ignore_axis_order=True)
):
name = str(ensemble.name)
names.add(name)
return frozenset(names)
def _mandatory_pipeline_grids(pipeline: _Pipeline) -> set[str]:
"""Grid files the compiled pipeline must read, under PROJ's own filenames.
A name prefixed with ``@`` is one PROJ will run without, so it proves
nothing about what is installed and is not returned here.
"""
found: set[str] = set()
for token in (pipeline.text or "").split():
if not token.startswith("grids="):
continue
found.update(
name
for name in token.removeprefix("grids=").split(",")
if name and not name.startswith("@")
)
return found
def _confirm_installed(
grids: tuple[GridUsage, ...], pipeline: _Pipeline
) -> tuple[GridUsage, ...]:
"""Correct the registry's account of what is installed with the pipeline's.
The registry names the grid the authority published; PROJ substitutes its
own distribution of it through proj.db's ``grid_alternatives``. EPSG:15851
cites ``conus.las`` and ``conus.los``, neither of which PROJ ships any
more: it reads the installed ``us_noaa_conus.tif`` instead, while the
registry still reports the published names as missing. Believing the
registry there refuses a transformation that demonstrably runs, which is
the same failure in the opposite direction to the one
:func:`_require_grids` exists to prevent.
The compiled pipeline is the complete statement of what will be read, so
it is the pipeline's files that are checked, and checked on disk: PROJ
keeps opened grids in memory, and a pipeline compiled while a file was
present still compiles once the file is gone, only to fail at the first
coordinate. A registry grid reported missing is taken as satisfied only
when every file the pipeline must read is installed; a file the pipeline
must read and that is not installed is reported missing under PROJ's own
name, whatever the registry said. With PROJ's network access enabled a
file can be fetched on demand, so its absence proves nothing and the
pipeline is trusted as compiled.
Args:
grids: What the registry says the applied operations depend on.
pipeline: The pipeline that was compiled from them.
Returns:
The same grids, with availability corrected from the pipeline and the
disk, plus one entry per file the pipeline reads that the registry did
not name and that is not installed.
"""
reads = _mandatory_pipeline_grids(pipeline)
if not reads:
return grids
if is_network_enabled():
# A grid PROJ substituted may be fetched on demand, so its absence
# proves nothing; a name the registry calls missing and the pipeline
# reads unchanged still says nothing new.
return tuple(
grid
if grid.available or grid.name in reads
else replace(grid, available=True)
for grid in grids
)
absent = sorted(name for name in reads if not _installed_grid(name))
if not absent:
return tuple(replace(grid, available=True) for grid in grids)
# The pipeline is the complete list of what will be read: a registry name
# it does not read is satisfied, and only its own absent files are missing.
named = {grid.name for grid in grids}
return (
*(replace(grid, available=grid.name not in absent) for grid in grids),
*(
GridUsage(
name=name,
full_name="",
package_name="",
url="",
available=False,
open_license=False,
direct_download=False,
)
for name in absent
if name not in named
),
)
def _installed_grid(name: str) -> bool:
"""Whether PROJ will find a grid file on disk under the name it uses.
Searched where PROJ searches: the data directories on its search path, in
order, and the user-writable directory that ``projsync`` and PROJ's own
network cache write to.
"""
path = Path(name)
if path.is_absolute():
return path.is_file()
return any((Path(directory) / name).is_file() for directory in _grid_directories())
def _grid_directories() -> tuple[str, ...]:
"""PROJ's grid search path, first match winning, without duplicates."""
found: list[str] = []
for entry in (
*datadir.get_data_dir().split(os.pathsep),
*os.environ.get("PROJ_DATA", "").split(os.pathsep),
*os.environ.get("PROJ_LIB", "").split(os.pathsep),
datadir.get_user_data_dir(),
):
if entry and entry not in found:
found.append(entry)
return tuple(found)
def _constituent_operations(
pipeline: _Pipeline, definition: dict[str, Any]
) -> tuple[CoordinateOperation, ...]:
"""Every operation being applied, so their grids can be inspected.
Read from what PROJ built, never from the EPSG registry. The registry names
the grid the authority published, while PROJ substitutes its own
distribution of it: EPSG:3858 cites
``Und_min2.5x2.5_egm2008_isw=82_WGS84_TideFree``, which is not installed,
where PROJ actually reads ``us_nga_egm08_25.tif``, which is. Asking the
registry would report a missing grid for a transformation that works.
"""
found: list[CoordinateOperation] = []
for transformer, _ in pipeline.steps:
found.extend(transformer.operations or ())
try:
found.append(CoordinateOperation.from_json_dict(transformer.to_json_dict()))
except (CRSError, TypeError, ProjError):
continue
return tuple(found)
def _require_grids(
grids: tuple[GridUsage, ...],
source: CoordinateReferenceSystem,
target: CoordinateReferenceSystem,
) -> None:
"""Refuse to transform when a grid the operation depends on is absent."""
missing = [grid for grid in grids if not grid.available]
if not missing:
return
described = ", ".join(
f"{grid.name}"
+ (f" (from {grid.package_name})" if grid.package_name else "")
+ (f" at {grid.url}" if grid.url else "")
for grid in missing
)
raise MissingGridError(
f"transforming {_label(source)} to {_label(target)} needs "
f"{len(missing)} grid file(s) that are not installed: {described}"
)
def _describe(
requests: tuple[OperationRequest, ...],
pipeline: _Pipeline,
definition: dict[str, Any],
operations: tuple[CoordinateOperation, ...],
*,
ballpark: bool,
requires_epoch: bool = False,
) -> AppliedOperation:
"""Record which operation was applied, against what was asked for.
Provenance comes from the requested operation's own node in the tree when
exactly one was named, or from the operation a bound CRS names for
itself. When more than one was named, no single node represents "the"
operation -- the whole applied step does, since that is what actually
fused them -- so the top-level definition is reported as-is. Otherwise it
comes from the substantive step: normalising axis order renames the
top-level operation, appending "(with axis order normalized for
visualization)" to its name, so using it directly would leak that wording
into a result the caller never asked to have annotated.
Narrowing to one node is refused outright when the caller named nothing,
the pipeline applies more than one datum transformation, and no authority
publishes the whole of it. That is what a bound CRS on each side of the
pair produces: each names its own shift to the hub, so the identification
one of them supplies would understate the result by the other. Both are
then reported as :attr:`AppliedOperation.bound_operations` instead, since
a chain of two identified operations has no single code but is not
unidentified either. A registered concatenated operation such as
EPSG:8047 is unaffected, since the identifier is then on the top-level
node, and so is an operation the caller named, which is reported as asked
for.
"""
node = definition
unnamed_chain = (
not requests
and datum_operation_count(definition) > 1
and _identifier(definition) is None
)
identifier_of = (
requests[0]
if len(requests) == 1
else pipeline.identified_by[0]
if len(pipeline.identified_by) == 1
else None
)
if not unnamed_chain:
if identifier_of is not None:
matched = identifier_of.find_in(definition)
if matched is not None:
node = matched
elif not requests:
substantive = _substantive_operation(operations)
if substantive is not None:
node = substantive.to_json_dict()
identifier = _identifier(node)
steps = tuple(
name for name in operation_names(definition) if name != definition.get("name")
)
method = node.get("method")
name = str(node.get("name") or pipeline.core.description).removesuffix(
_VISUALIZATION_SUFFIX
)
if pipeline.core_direction is TransformDirection.INVERSE:
name = f"Inverse of {name}"
return AppliedOperation(
requested=None if not requests else " + ".join(r.text for r in requests),
auth_name=None if identifier is None else identifier[0],
code=None if identifier is None else identifier[1],
name=name,
method_name=(
None
if unnamed_chain
else str(method["name"])
if isinstance(method, dict) and "name" in method
else _substantive_method(operations)
),
accuracy=pipeline.accuracy,
route=pipeline.route,
ballpark=ballpark,
requires_epoch=requires_epoch,
steps=tuple(sorted(steps)),
projjson=json.dumps(node),
execution_direction=pipeline.core_direction,
bound_operations=tuple(r.text for r in pipeline.identified_by),
axis_order_corrected=pipeline.corrected is not None,
)
def _substantive_operation(
operations: tuple[CoordinateOperation, ...],
) -> CoordinateOperation | None:
"""The operation that does the geodetic work, not the axis bookkeeping.
Normalising axis order makes PROJ prepend an axis-reversal conversion, so
the first step of a concatenated operation is often bookkeeping rather than
the transformation the caller cares about.
"""
for operation in operations:
if not _is_bookkeeping(operation):
return operation
return None
def _substantive_method(operations: tuple[CoordinateOperation, ...]) -> str | None:
"""Name of the method applied, or None when more than one operation is chained.
A chain has no single method, and naming the first one would understate
what the other steps do.
"""
substantive = [
operation for operation in operations if not _is_bookkeeping(operation)
]
if len(substantive) != 1 or not substantive[0].method_name:
return None
return str(substantive[0].method_name)
def _identifier(definition: dict[str, Any]) -> tuple[str, str] | None:
"""The authority code of the operation as a whole, if it has one."""
identity = definition.get("id")
if isinstance(identity, dict):
authority = identity.get("authority")
code = identity.get("code")
if authority is not None and code is not None:
return base_authority(str(authority)), str(code)
return None
def _apply(
transformer: Transformer,
values: list[list[float]],
epoch: float | None,
direction: TransformDirection,
) -> list[list[float]]:
"""Run one transformer over every coordinate component at once."""
count = len(values[0]) if values else 0
arguments: dict[str, Any] = {
"xx": values[0],
"yy": values[1] if len(values) > 1 else [0.0] * count,
}
if len(values) > 2:
arguments["zz"] = values[2]
if epoch is not None:
arguments["tt"] = [epoch] * count
produced = transformer.transform(**arguments, direction=direction, errcheck=True)
return [list(component) for component in produced[: len(values)]]
def _compose_pipeline(
steps: Sequence[tuple[Transformer, TransformDirection]],
) -> str | None:
"""Write a sequence of transformers out as one runnable PROJ pipeline.
Running several pipelines back to back is the same as running one pipeline
holding all of their steps in order, so the chain can be stated as a single
definition rather than as a list of parts a caller would have to reassemble
to reproduce a result.
Args:
steps: The transformers and the direction each one is applied in.
Returns:
A definition :meth:`pyproj.Transformer.from_pipeline` accepts, or None
when a step carries pipeline-level parameters, which PROJ applies to
every step of their own pipeline and which would silently reach the
other steps once merged.
"""
composed: list[str] = []
for transformer, direction in steps:
flattened = _pipeline_steps(transformer.definition)
if flattened is None:
return None
if direction is TransformDirection.INVERSE:
flattened = [_inverted_step(step) for step in reversed(flattened)]
composed.extend(flattened)
if not composed:
return "proj=noop"
return " ".join(["proj=pipeline", *(f"step {step}" for step in composed)])
def _pipeline_steps(definition: str) -> list[str] | None:
"""Split one PROJ definition into its steps, without the ``step`` keywords.
Returns:
One entry per step, a definition that is not a pipeline being a single
step in itself, or None when the pipeline carries parameters outside
any step.
"""
tokens = [token.removeprefix("+") for token in definition.split()]
if not tokens:
return None
if tokens[0] != "proj=pipeline":
return [" ".join(tokens)]
steps: list[list[str]] = []
for token in tokens[1:]:
if token == "step":
steps.append([])
elif not steps:
return None
else:
steps[-1].append(token)
return [" ".join(step) for step in steps if step]
def _inverted_step(step: str) -> str:
"""The same PROJ step run backwards, by adding or removing its ``inv`` flag."""
tokens = step.split()
if _INVERSE_FLAG in tokens:
tokens.remove(_INVERSE_FLAG)
return " ".join(tokens)
return f"{_INVERSE_FLAG} {step}"
def _output_indices(
target: CoordinateReferenceSystem, produced: int
) -> tuple[int, ...]:
"""Map PROJ's output components onto the axes the target CRS declares.
PROJ returns as many components as it was given. The target CRS decides how
many of them are coordinates in that CRS: a vertical CRS declares one axis
and its value is the height component, not the first one. A source height
supplied alongside a horizontal-only pair is one component more than the
target declares; it is passed through rather than dropped, since the
caller gave it deliberately and PROJ already carried it through unchanged.
"""
if target.dimension == 1 and _is_vertical(target):
if produced < 3:
raise TransformationFailedError(
f"{_label(target)} declares a height axis, but only {produced} "
"coordinate components were supplied; give the source height too"
)
return (2,)
if target.dimension > produced:
raise TransformationFailedError(
f"{_label(target)} declares {target.dimension} axes but only "
f"{produced} coordinate components were produced; supply "
f"{target.dimension} values per point"
)
if produced - target.dimension > 1:
raise TransformationFailedError(
f"{_label(target)} declares {target.dimension} axes but {produced} "
"coordinate components were produced; at most one extra "
"(a pass-through height) is carried through"
)
return tuple(range(produced))
def _is_vertical(crs: CoordinateReferenceSystem) -> bool:
"""Whether the CRS's single axis is a height or a depth."""
return crs.axes[0].direction.lower() in _VERTICAL_DIRECTIONS
# ---------------------------------------------------------------------------
# PROJ workaround, not permanent design. ``always_xy`` is supposed to
# guarantee lon/lat, E/N in and out; it does not for a pipeline with a vertical
# (or otherwise horizontal-axis-less) end, because there is no declared order
# there for PROJ to normalise against, and the operation's own residual
# ``axisswap`` survives. ``_entry_step``, ``_swaps_horizontal``,
# ``reads_declared_horizontal`` and ``_order_horizontal_for_pipeline`` exist
# only to detect and undo that residual swap on the way in.
#
# Once PROJ strips the residual axisswap itself,
# ``Transformation.transform`` no longer needs to call
# ``_order_horizontal_for_pipeline`` at all, and this whole block --
# including ``reads_declared_horizontal`` and its call site -- can be deleted
# outright rather than adapted: nothing else in this module depends on it.
# ---------------------------------------------------------------------------
def _entry_step(definition: str, direction: TransformDirection) -> str:
"""The step a PROJ pipeline applies first when run in ``direction``.
An inverted pipeline runs its steps back to front, so its entry step is the
last one written. ``axisswap`` is its own inverse, so the step reads the
same either way.
"""
steps = definition.split(" step ")
if len(steps) == 1:
return definition
return steps[1] if direction is TransformDirection.FORWARD else steps[-1]
def _swaps_horizontal(step: str) -> bool:
"""Whether a PROJ step exchanges the first two coordinate components."""
if "proj=axisswap" not in step:
return False
for token in step.split():
if token.startswith("order="):
order = token.removeprefix("order=").split(",")
return order[:2] == ["2", "1"]
return False
def _order_horizontal_for_pipeline(
source: CoordinateReferenceSystem,
pipeline: _Pipeline,
columns: tuple[tuple[float, ...], ...],
) -> tuple[tuple[tuple[float, ...], ...], bool]:
"""Hand a vertical CRS's position to PROJ in the order the pipeline reads it.
Coordinate values reach this package in ``xy`` order, and PROJ's
``always_xy`` normally guarantees PROJ reads them that way. It cannot for a
vertical source: with no horizontal axes to normalise, the pipeline keeps
the ``axisswap`` belonging to the operation's own geographic end and reads
the accompanying position latitude-first. Left uncorrected the grid is
interpolated at the transposed point, which is a wrong height wherever the
transposed point is still inside the grid.
Args:
source: Source CRS.
pipeline: The pipeline the values are about to be run through.
columns: One tuple of values per axis, in ``xy`` order.
Returns:
The columns with the horizontal pair transposed when the pipeline reads
it northing-first, and unchanged otherwise, and whether it transposed
them. The caller reports the transposition as a step of the pipeline,
so it needs to be told rather than left to compare tuples.
"""
if len(columns) != _MAX_COORDINATE_VALUES:
return columns, False
if not (source.dimension == 1 and _is_vertical(source)):
return columns, False
if not pipeline.reads_declared_horizontal:
return columns, False
return (columns[1], columns[0], columns[2]), True
def _report_horizontal_swap(pipeline: str | None) -> str | None:
"""Write the transposition above into the pipeline that gets reported.
The reordering happens in Python, before PROJ is handed the values, so a
pipeline rebuilt from the reported text alone would read the caller's
``xy`` values transposed and interpolate the grid at the wrong point.
Stating it as the ``axisswap`` step it is keeps the reported pipeline
reproducing the result from the values the caller actually gave.
"""
if pipeline is None:
return None
body = pipeline.removeprefix("proj=pipeline ")
if body == pipeline:
body = f"step {pipeline}"
return f"proj=pipeline step proj=axisswap order=2,1 {body}"
# ---------------------------------------------------------------------------
# PROJ/EPSG workaround, not permanent design. Two things go wrong at an
# engineering CRS, and both transpose coordinates without any error:
#
# 1. ``always_xy`` never normalises an engineering CRS. PROJ's
# ``mustAxisOrderBeSwitchedForVisualization`` considers geographic,
# projected and derived projected CRSs only, so a plant grid declared
# (northing, easting) is read and written in that order while every other
# end of the pipeline is east-first. This package promises ``xy`` wherever
# a CRS has an easting and a northing, so the ``axisswap`` PROJ omits is
# added here. Cartesian Grid Offsets (EPSG:9656) is the exception: PROJ
# renders it as an east-first ``affine`` and adapts projected ends around
# it but leaves engineering ends alone, so at such an end the values PROJ
# consumes already are ``xy``. A grid with no east/north pair at all
# (EPSG:5800 declares north and west) is left in declared order, which is
# what ``value_axis_order`` reports for it and what PROJ does with the
# corresponding projected case (Krovak's south and west).
#
# 2. PROJ reads the "Ordinate 1/2 of evaluation point in target CRS" of a
# Similarity transformation (EPSG:9621), and A0/B0 of an Affine parametric
# transformation (EPSG:9624), as the target CRS's first and second declared
# axes, then appends an ``axisswap`` to normalise a northing-first target.
# EPSG's own data does not follow that reading: EPSG:1035 states the Astra
# Minas grid origin as (2610200.48, 4905282.73) in EPSG:22192, which is
# declared (northing, easting) -- read that way the origin lies in
# Antarctica, read as (easting, northing) it lies at Comodoro Rivadavia,
# the operation's area of use. Registers of plant grids are authored both
# ways. So the convention is established from the data: the evaluation
# point is unprojected under both readings and the one that falls inside
# the area of use wins. When the ordinates turn out east-first, PROJ's
# appended ``axisswap`` transposes a result that was already ``xy`` and is
# removed. When neither or both readings are plausible the operation is
# refused rather than guessed at. An evaluation point stated in an
# engineering CRS cannot be unprojected, so PROJ's declared-order reading
# is kept there.
#
# Retire part 1 when PROJ normalises engineering CRSs under always_xy; the
# canary is tests/geodesy/test_engineering_axes.py::
# test_proj_still_leaves_engineering_axes_alone, which fails the day it does.
# Retire part 2 when EPSG and PROJ agree on which axis each ordinate of
# methods 9621/9624 refers to. The corrected core is rebuilt from its text
# and reported as TransformationResult.pipeline, so replaying that pipeline
# reproduces the result from the caller's own ``xy`` values.
# ---------------------------------------------------------------------------
_AXIS_SWAP_STEP = "proj=axisswap order=2,1"
_AFFINE_STEP = "proj=affine"
_INVERSE_METHOD_PREFIX = "inverse of "
_EVALUATION_POINT_METHODS = frozenset(
{"similarity transformation", "affine parametric transformation"}
)
_GRID_OFFSETS_METHOD = "cartesian grid offsets"
# Slack around an area of use when placing an evaluation point, in degrees.
# Extents are quoted coarsely; the two readings differ by thousands of
# kilometres, so this only has to absorb a sloppy bounding box.
_AREA_MARGIN_DEGREES = 1.0
def _correct_engineering_axes(pipeline: _Pipeline) -> _Pipeline:
"""Rewrite the core so that its engineering ends honour the ``xy`` contract.
See the block comment above for what is corrected and why.
Returns:
The pipeline unchanged when no engineering CRS and no northing-first
evaluation point is involved, otherwise with ``corrected`` set to the
steps to execute.
Raises:
OperationNotAvailableError: If the operation is a Similarity or Affine
parametric transformation whose evaluation point is stated in a
northing-first projected CRS and cannot be placed inside the area
of use under exactly one reading of its ordinates, or if PROJ
built it in a shape this correction does not recognise.
"""
nodes = _substantive_nodes(pipeline.definition)
steps = _pipeline_steps(pipeline.core.definition)
if not nodes or steps is None:
return pipeline
rewritten = list(steps)
changed = _drop_spurious_ordinate_swap(rewritten, nodes, pipeline.core)
changed = _add_engineering_swaps(rewritten, nodes) or changed
if not changed:
return pipeline
corrected = Transformer.from_pipeline(
" ".join(["proj=pipeline", *(f"step {step}" for step in rewritten)])
)
return replace(
pipeline,
corrected=tuple(
(corrected if transformer is pipeline.core else transformer, direction)
for transformer, direction in pipeline.steps
),
)
def _substantive_nodes(definition: dict[str, Any]) -> list[dict[str, Any]]:
"""The operations in a PROJJSON tree that move coordinates, in order.
``always_xy`` wraps the operation in a concatenated operation together
with the axis-order bookkeeping it inserts; that bookkeeping is skipped.
"""
if definition.get("type") == "ConcatenatedOperation":
return [
node
for node in definition.get("steps", ())
if isinstance(node, dict)
and _method_of(node)[0] not in _BOOKKEEPING_METHODS
]
return [definition] if isinstance(definition.get("method"), dict) else []
def _method_of(node: dict[str, Any]) -> tuple[str, bool]:
"""A node's method name, lower-cased, and whether PROJ inverted the node.
PROJ exports an inverted step with its forward parameters, the ends
swapped, and ``Inverse of`` prefixed to the method name.
"""
method = node.get("method")
name = str(method.get("name", "")).lower() if isinstance(method, dict) else ""
if name.startswith(_INVERSE_METHOD_PREFIX):
return name[len(_INVERSE_METHOD_PREFIX) :], True
return name, False
def _node_crs(node: dict[str, Any], end: str) -> CRS | None:
"""The CRS a PROJJSON operation node names at ``end``, as declared."""
value = node.get(end)
if not isinstance(value, dict):
return None
try:
return CRS.from_json_dict(value)
except CRSError:
return None
def _north_first(crs: CRS) -> bool:
"""Whether the CRS declares a northing before an easting.
Read through :attr:`CoordinateReferenceSystem.value_axis_order` so that
the correction and the order this package reports cannot disagree.
"""
order = CoordinateReferenceSystem(crs, crs.name).value_axis_order
return order[:2] == (1, 0)
def _add_engineering_swaps(steps: list[str], nodes: list[dict[str, Any]]) -> bool:
"""Part 1: swap the values at a northing-first engineering end.
Returns:
Whether ``steps`` was changed.
"""
changed = False
entry = _node_crs(nodes[0], "source_crs")
if (
entry is not None
and entry.is_engineering
and _north_first(entry)
and _method_of(nodes[0])[0] != _GRID_OFFSETS_METHOD
):
steps.insert(0, _AXIS_SWAP_STEP)
changed = True
exit_ = _node_crs(nodes[-1], "target_crs")
if (
exit_ is not None
and exit_.is_engineering
and _north_first(exit_)
and _method_of(nodes[-1])[0] != _GRID_OFFSETS_METHOD
):
steps.append(_AXIS_SWAP_STEP)
changed = True
return changed
def _drop_spurious_ordinate_swap(
steps: list[str], nodes: list[dict[str, Any]], core: Transformer
) -> bool:
"""Part 2: remove PROJ's normalising swap when the ordinates are east-first.
Returns:
Whether ``steps`` was changed.
Raises:
OperationNotAvailableError: See :func:`_correct_engineering_axes`.
"""
evaluated = []
for node in nodes:
method, inverted = _method_of(node)
if method not in _EVALUATION_POINT_METHODS:
continue
# The evaluation point is stated in the operation's own target CRS,
# which is the node's source once PROJ has inverted it.
crs = _node_crs(node, "source_crs" if inverted else "target_crs")
if crs is not None and crs.is_projected and _north_first(crs):
evaluated.append((node, crs))
if not evaluated:
return False
if len(nodes) > 1:
raise OperationNotAvailableError(
f"{_node_label(evaluated[0][0])} states its evaluation point in "
f"{evaluated[0][1].name!r}, which declares its northing first, and "
"PROJ chained it with other steps; whether the ordinates follow "
"that declared order or are east-first cannot be established for "
"a step inside a chain, so the operation is refused rather than "
"risk a transposed result. Apply the operation between its own "
"CRSs instead"
)
node, crs = evaluated[0]
affine = [index for index, step in enumerate(steps) if _AFFINE_STEP in step.split()]
if len(affine) != 1:
raise OperationNotAvailableError(
f"PROJ built {_node_label(node)} as {len(affine)} affine steps "
"rather than one, a shape this package does not recognise; "
"refusing rather than risk a transposed result"
)
index = affine[0]
inverted = _INVERSE_FLAG in steps[index].split()
first, second = _affine_offsets(steps[index])
convention = _ordinate_convention(
crs, first, second, _area_bounds(core.area_of_use, crs, node)
)
if convention is None:
raise OperationNotAvailableError(
f"{_node_label(node)} states its evaluation point as ({first:g}, "
f"{second:g}) in {crs.name!r}, which declares its northing first, "
"and the point does not fall inside the area of use under exactly "
"one reading of those ordinates (declared order, or easting "
"first); which axis each ordinate refers to cannot be established, "
"so the operation is refused rather than risk a transposed result"
)
if convention == "declared":
return False
swap = index - 1 if inverted else index + 1
if not (0 <= swap < len(steps) and _swaps_horizontal(steps[swap])):
raise OperationNotAvailableError(
f"PROJ built {_node_label(node)} without the axis swap this "
"package expects beside its affine step; refusing rather than "
"risk a transposed result"
)
del steps[swap]
return True
def _affine_offsets(step: str) -> tuple[float, float]:
"""The translation of a PROJ ``affine`` step, in the CRS's own units."""
offsets = {"xoff": 0.0, "yoff": 0.0}
for token in step.split():
key, separator, value = token.partition("=")
if separator and key in offsets:
offsets[key] = float(value)
return offsets["xoff"], offsets["yoff"]
def _area_bounds(
area: Any, crs: CRS, node: dict[str, Any]
) -> tuple[float, float, float, float] | None:
"""The bounding box to place an evaluation point against.
The operation's own area of use is the most specific. Failing that, the
CRS's -- read from the registry by the authority code the operation names
for it when the embedded definition dropped its usage, as PROJ's export of
an operation's ends does.
"""
bounds = getattr(area, "bounds", None)
if bounds is not None:
return cast("tuple[float, float, float, float]", bounds)
if crs.area_of_use is not None:
return crs.area_of_use.bounds
for end in ("target_crs", "source_crs"):
value = node.get(end)
if not isinstance(value, dict) or value.get("name") != crs.name:
continue
identifier = _identifier(value)
if identifier is None:
return None
try:
registered = CRS.from_authority(*identifier)
except CRSError:
return None
if registered.area_of_use is not None and registered.equals(
crs, ignore_axis_order=False
):
return registered.area_of_use.bounds
return None
def _ordinate_convention(
crs: CRS,
first: float,
second: float,
bounds: tuple[float, float, float, float] | None,
) -> str | None:
"""Which axes an evaluation point's ordinates refer to, from where it lands.
Args:
crs: The northing-first projected CRS the point is stated in.
first: The first ordinate as stated.
second: The second ordinate as stated.
bounds: The area of use to place the point against, or None when
nothing states one.
Returns:
``"declared"`` when only the declared reading (first ordinate a
northing) lands inside the area of use, ``"east-first"`` when only
the other does, None when neither or both do or there is no area.
"""
geodetic = crs.geodetic_crs
if bounds is None or geodetic is None:
return None
to_geographic = Transformer.from_crs(crs, geodetic, always_xy=True)
# An area of use is stated from Greenwich; a base CRS on another prime
# meridian (NGO 1948 counts from Oslo) reports longitudes from there.
meridian = geodetic.prime_meridian
offset = (
0.0
if meridian is None
else math.degrees(
float(meridian.longitude) * float(meridian.unit_conversion_factor)
)
)
def lands(easting: float, northing: float) -> bool:
longitude, latitude = to_geographic.transform(easting, northing)
return _inside(bounds, longitude + offset, latitude)
declared = lands(second, first)
east_first = lands(first, second)
if declared and not east_first:
return "declared"
if east_first and not declared:
return "east-first"
return None
def _inside(
bounds: tuple[float, float, float, float], longitude: float, latitude: float
) -> bool:
"""Whether a position lies within a bounding box, allowing for slack."""
if not (math.isfinite(longitude) and math.isfinite(latitude)):
return False
longitude = (longitude + 180.0) % 360.0 - 180.0
west, south, east, north = bounds
margin = _AREA_MARGIN_DEGREES
if not south - margin <= latitude <= north + margin:
return False
if west <= east:
return west - margin <= longitude <= east + margin
# A box across the antimeridian.
return longitude >= west - margin or longitude <= east + margin
def _node_label(node: dict[str, Any]) -> str:
"""Short identification of an operation node for error messages."""
identifier = _identifier(node)
name = node.get("name") or "the operation"
return f"{name!r}" if identifier is None else f"{identifier[0]}:{identifier[1]}"
def _carry_unread_horizontal(
source: CoordinateReferenceSystem,
target: CoordinateReferenceSystem,
applied: AppliedOperation,
grids: tuple[GridUsage, ...],
columns: tuple[tuple[float, ...], ...],
) -> tuple[tuple[float, ...], ...]:
"""Withhold a horizontal position that the applied shift never reads.
A height travelling between two vertical CRSs is usually accompanied by
the horizontal position its grid is read at. Where the shift is instead
the same everywhere -- a plain vertical offset, a unit change, a
height/depth reversal -- that position is not part of the calculation,
and it need not even be geographic: it is commonly an engineering or
projected coordinate, which PROJ's ``geogoffset`` step would reject as an
impossible latitude. Both CRSs being vertical, the position cannot reach
the result either way, so it is held back rather than offered to PROJ as
something it is not.
Args:
source: Source CRS.
target: Target CRS.
applied: The operation being applied.
grids: Grids that operation depends on.
columns: One tuple of values per axis.
Returns:
The columns, with the horizontal pair blanked when it cannot be read,
and unchanged otherwise.
"""
if len(columns) != _MAX_COORDINATE_VALUES or grids:
return columns
if not (source.dimension == 1 and _is_vertical(source)):
return columns
if not (target.dimension == 1 and _is_vertical(target)):
return columns
if (applied.method_name or "").strip().lower() not in (
_POSITION_FREE_VERTICAL_METHODS
):
return columns
blank = (0.0,) * len(columns[0])
return (blank, blank, columns[2])
def _require_in_range(
source: CoordinateReferenceSystem, columns: tuple[tuple[float, ...], ...]
) -> None:
"""Refuse a latitude or longitude a geographic CRS's axis unit cannot represent.
PROJ rejects an impossible latitude too, but only as "Invalid latitude",
naming neither the CRS nor the units nor which value it read as latitude.
Since the usual cause is projected coordinates in metres handed to a
geographic CRS, or latitude passed first, the message has to name all
three to be actionable. An impossible longitude PROJ does not reject at
all: it wraps it, so a longitude of 400 degrees quietly becomes 40.
The limits are derived from each axis's own unit rather than assumed to
be 90 and 180, so a CRS declaring grads (``EPSG:4807``) is held to 100
rather than wrongly refused. Longitude is allowed a full turn either way,
since a dataset counted from 0 to 360 is a convention, not a mistake.
"""
if not source.is_geographic:
return
order = source.value_axis_order
checks: list[tuple[int, AxisSpec, float]] = []
for value_index, declared_index in enumerate(order[: len(columns)]):
axis = source.axes[declared_index]
direction = axis.direction.lower()
if direction in _NORTHINGS:
limit = (math.pi / 2) / axis.unit_conversion_factor
checks.insert(0, (value_index, axis, limit))
elif direction in _EASTINGS:
limit = (2 * math.pi) / axis.unit_conversion_factor
checks.append((value_index, axis, limit))
# Latitude first: when both values are metres, the one read as latitude is
# the diagnostic the caller can act on.
for value_index, axis, limit in checks:
for point, value in enumerate(columns[value_index]):
if math.isfinite(value) and abs(value) > limit:
raise CoordinateOutOfRangeError(
f"point {point} has {axis.name.lower()} {value} "
f"{axis.unit_name}, outside the valid range "
f"[-{limit:g}, {limit:g}] for {_label(source)}; values are "
f"given in {source.value_axis_abbreviations} order and in "
f"{source.axis_units} -- projected coordinates in metres "
"need a projected CRS"
)
def _require_finite(
rows: tuple[tuple[float, ...], ...],
source: CoordinateReferenceSystem,
target: CoordinateReferenceSystem,
) -> None:
"""Refuse to return an infinite or undefined coordinate.
PROJ signals an unrepresentable result with infinity. Returned as-is it
would propagate silently into whatever consumes it.
"""
for index, row in enumerate(rows):
if not all(math.isfinite(value) for value in row):
raise TransformationFailedError(
f"point {index} has no finite representation transforming "
f"{_label(source)} to {_label(target)}; PROJ produced {row}"
)
def _label(crs: CoordinateReferenceSystem) -> str:
"""Short identification of a CRS for error messages."""
return crs.authority_code or crs.name
def _applied_label(definition: dict[str, Any]) -> str:
"""Short identification of the operation PROJ actually built."""
identifier = _identifier(definition)
name = definition.get("name") or "an unnamed operation"
if identifier is None:
return f"{name!r}"
return f"{identifier[0]}:{identifier[1]} ({name!r})"