CRS & Reprojection¶

Assign or read a coordinate reference system (via rioxarray / pyproj), reproject a dataset, convert between geographic (lon/lat) and projected (x/y) coordinates, and express coordinates in a metric local frame anchored at a site (UTM zone picked from the origin).

Operators¶

AssignLocalXY ¶

AssignLocalXY(frame: LocalFrame | dict[str, Any], *, lon: str = 'lon', lat: str = 'lat', x: str = 'x', y: str = 'y')

Bases: Operator

Add metric x/y coordinates relative to a :class:LocalFrame.

Wraps :func:xrops.geom.assign_local_xy — the point / swath counterpart of :class:ReprojectMatch. The frame is stored as its JSON payload, so get_config() round-trips through the constructor.

Parameters:

Name	Type	Description	Default
`frame`	`LocalFrame \| dict[str, Any]`	A :class:`xrops.geom.LocalFrame` or its `to_dict()` payload (`{"crs", "origin_lon", "origin_lat"}`).	required
`lon`	`str`	Longitude coordinate name.	`'lon'`
`lat`	`str`	Latitude coordinate name.	`'lat'`
`x`	`str`	Name of the output east coordinate.	`'x'`
`y`	`str`	Name of the output north coordinate.	`'y'`

Returns:

Type	Description
	The dataset with `x`/`y` coordinates in metres (2-D for grids
	and swaths, 1-D for along-track points) and `attrs["local_frame"]`
	set.

Source code in packages/xrtoolz-ops/src/xrops/geom/operators.py

def __init__(
    self,
    frame: _crs.LocalFrame | dict[str, Any],
    *,
    lon: str = "lon",
    lat: str = "lat",
    x: str = "x",
    y: str = "y",
):
    self.frame = (
        frame
        if isinstance(frame, _crs.LocalFrame)
        else _crs.LocalFrame.from_dict(frame)
    )
    self.lon = lon
    self.lat = lat
    self.x = x
    self.y = y

Local metric frames¶

local_frame(origin_lon, origin_lat) picks the UTM zone containing the origin and returns a frozen, JSON-serialisable LocalFrame whose to_xy gives metres east / north of the origin and whose to_lonlat inverts it. assign_local_xy stamps those coordinates onto a dataset (1-D rectilinear, 2-D swath, or paired 1-D along-track lon/lat on one dim); calc_latlon is the inverse direction for gridded rasters that already carry projected x/y axes.

utm_crs_for ¶

utm_crs_for(lon: float, lat: float, *, datum: str = 'WGS 84') -> str

CRS authority code of the UTM zone containing (lon, lat).

The zone number follows the standard 6° bands plus the Norway (zone 32 widened over 56–64°N) and Svalbard (zones 31/33/35/37 over 72–84°N) exceptions. A zone-boundary meridian belongs to the zone east of it (lon=-102.0 → zone 14), longitudes are wrapped into [-180, 180) first (so 0–360 grids work), and lat=0 counts as the northern hemisphere. The code itself comes from :func:pyproj.database.query_utm_crs_info for datum (a point :class:~pyproj.aoi.AreaOfInterest on the zone's central meridian); when the database has no entry the WGS 84 convention "EPSG:326NN" (north) / "EPSG:327NN" (south) is used.

Parameters:

Name	Type	Description	Default
`lon`	`float`	Longitude in degrees east.	required
`lat`	`float`	Latitude in degrees north.	required
`datum`	`str`	Datum name as spelled in the PROJ database (e.g. `"WGS 84"`, `"NAD83"`, `"ETRS89"`).	`'WGS 84'`

Returns:

Type	Description
`str`	`"<AUTH>:<code>"` string, e.g. `"EPSG:32613"`.

Raises:

Type	Description
`ValueError`	If `lat` is outside `[-80, 84]` (the UTM domain — use a polar stereographic CRS instead), or if `datum` is not `"WGS 84"` and the database has no UTM CRS for that zone.

Example

from xrops.geom import utm_crs_for utm_crs_for(-102.5, 31.5) 'EPSG:32613' utm_crs_for(151.2, -33.9) 'EPSG:32756'

Source code in packages/xrtoolz-ops/src/xrops/geom/_src/crs.py

def utm_crs_for(lon: float, lat: float, *, datum: str = "WGS 84") -> str:
    """CRS authority code of the UTM zone containing ``(lon, lat)``.

    The zone number follows the standard 6° bands plus the Norway
    (zone 32 widened over 56–64°N) and Svalbard (zones 31/33/35/37 over
    72–84°N) exceptions. A zone-boundary meridian belongs to the zone
    east of it (``lon=-102.0`` → zone 14), longitudes are wrapped into
    ``[-180, 180)`` first (so ``0–360`` grids work), and ``lat=0`` counts
    as the northern hemisphere. The code itself comes from
    :func:`pyproj.database.query_utm_crs_info` for ``datum`` (a point
    :class:`~pyproj.aoi.AreaOfInterest` on the zone's central meridian);
    when the database has no entry the WGS 84 convention
    ``"EPSG:326NN"`` (north) / ``"EPSG:327NN"`` (south) is used.

    Args:
        lon: Longitude in degrees east.
        lat: Latitude in degrees north.
        datum: Datum name as spelled in the PROJ database
            (e.g. ``"WGS 84"``, ``"NAD83"``, ``"ETRS89"``).

    Returns:
        ``"<AUTH>:<code>"`` string, e.g. ``"EPSG:32613"``.

    Raises:
        ValueError: If ``lat`` is outside ``[-80, 84]`` (the UTM domain —
            use a polar stereographic CRS instead), or if ``datum`` is not
            ``"WGS 84"`` and the database has no UTM CRS for that zone.

    Example:
        >>> from xrops.geom import utm_crs_for
        >>> utm_crs_for(-102.5, 31.5)
        'EPSG:32613'
        >>> utm_crs_for(151.2, -33.9)
        'EPSG:32756'

    """
    lat = float(lat)
    # Chained comparison is False for NaN, so NaN is rejected too.
    if not _UTM_LAT_MIN <= lat <= _UTM_LAT_MAX:
        raise ValueError(
            "utm_crs_for: UTM is only defined for "
            f"{_UTM_LAT_MIN}° <= lat <= {_UTM_LAT_MAX}°, got lat={lat}. "
            "Use a polar stereographic CRS for the poles."
        )
    lon = (float(lon) + 180.0) % 360.0 - 180.0
    zone = _utm_zone(lon, lat)
    hemisphere = "N" if lat >= 0.0 else "S"
    central_meridian = (zone - 1) * 6.0 - 180.0 + 3.0
    infos = query_utm_crs_info(
        datum_name=datum,
        area_of_interest=AreaOfInterest(central_meridian, lat, central_meridian, lat),
    )
    suffix = f"zone {zone}{hemisphere}"
    for info in infos:
        if info.name.endswith(suffix):
            return f"{info.auth_name}:{info.code}"
    if datum != "WGS 84":
        raise ValueError(
            f"utm_crs_for: no UTM {suffix} CRS found for datum {datum!r} in the "
            "PROJ database."
        )
    return f"EPSG:{326 if hemisphere == 'N' else 327}{zone:02d}"

LocalFrame `dataclass` ¶

LocalFrame(crs: str, origin_lon: float, origin_lat: float)

A metric frame: projected CRS plus an origin expressed in that CRS.

to_xy returns metres east / north of the origin in crs; to_lonlat inverts it. The frame holds only plain strings and floats, so it is hashable and JSON-serialisable (to_dict / from_dict); the :class:pyproj.Transformer pair is cached at module level per CRS.

Attributes:

Name	Type	Description
`crs`	`str`	A projected CRS with metre axes, as any specifier accepted by :class:`pyproj.CRS`, kept in string form (e.g. `"EPSG:32613"`).
`origin_lon`	`float`	WGS-84 longitude of the origin, degrees east.
`origin_lat`	`float`	WGS-84 latitude of the origin, degrees north.
`origin_xy`	`tuple[float, float]`	`(x0, y0)` of the origin projected into `crs`, computed at construction.

Raises:

Type	Description
`ValueError`	If `crs` is geographic or its axes are not in metres (e.g. `"EPSG:2263"`, US survey feet).

Example

from xrops.geom import local_frame frame = local_frame(-102.5, 31.5) frame.crs 'EPSG:32613' x, y = frame.to_xy(-102.5, 31.5) float(x), float(y) (0.0, 0.0)

to_xy ¶

to_xy(lon: float | Sequence[float] | ndarray, lat: float | Sequence[float] | ndarray) -> tuple[Float[ndarray, '...'], Float[ndarray, '...']]

Project WGS-84 lon/lat to metres relative to the origin.

Parameters:

Name	Type	Description	Default
`lon`	`float \| Sequence[float] \| ndarray`	Longitude(s) in degrees east (any shape).	required
`lat`	`float \| Sequence[float] \| ndarray`	Latitude(s) in degrees north (same shape as `lon`).	required

Returns:

Type	Description
`tuple[Float[ndarray, '...'], Float[ndarray, '...']]`	`(x, y)` arrays of metres east / north of the origin.

Source code in packages/xrtoolz-ops/src/xrops/geom/_src/crs.py

def to_xy(
    self,
    lon: float | Sequence[float] | np.ndarray,
    lat: float | Sequence[float] | np.ndarray,
) -> tuple[Float[np.ndarray, "..."], Float[np.ndarray, "..."]]:
    """Project WGS-84 ``lon``/``lat`` to metres relative to the origin.

    Args:
        lon: Longitude(s) in degrees east (any shape).
        lat: Latitude(s) in degrees north (same shape as ``lon``).

    Returns:
        ``(x, y)`` arrays of metres east / north of the origin.
    """
    forward, _ = _frame_transformers(self.crs)
    x, y = forward.transform(lon, lat)
    x0, y0 = self.origin_xy
    return np.asarray(x, dtype=float) - x0, np.asarray(y, dtype=float) - y0

to_lonlat ¶

to_lonlat(x: float | Sequence[float] | ndarray, y: float | Sequence[float] | ndarray) -> tuple[Float[ndarray, '...'], Float[ndarray, '...']]

Inverse of :meth:to_xy: metres from the origin to WGS-84 lon/lat.

Parameters:

Name	Type	Description	Default
`x`	`float \| Sequence[float] \| ndarray`	Metres east of the origin (any shape).	required
`y`	`float \| Sequence[float] \| ndarray`	Metres north of the origin (same shape as `x`).	required

Returns:

Type	Description
`tuple[Float[ndarray, '...'], Float[ndarray, '...']]`	`(lon, lat)` arrays in degrees.

Source code in packages/xrtoolz-ops/src/xrops/geom/_src/crs.py

def to_lonlat(
    self,
    x: float | Sequence[float] | np.ndarray,
    y: float | Sequence[float] | np.ndarray,
) -> tuple[Float[np.ndarray, "..."], Float[np.ndarray, "..."]]:
    """Inverse of :meth:`to_xy`: metres from the origin to WGS-84 lon/lat.

    Args:
        x: Metres east of the origin (any shape).
        y: Metres north of the origin (same shape as ``x``).

    Returns:
        ``(lon, lat)`` arrays in degrees.
    """
    _, inverse = _frame_transformers(self.crs)
    x0, y0 = self.origin_xy
    lon, lat = inverse.transform(
        np.asarray(x, dtype=float) + x0, np.asarray(y, dtype=float) + y0
    )
    return np.asarray(lon, dtype=float), np.asarray(lat, dtype=float)

to_dict ¶

to_dict() -> dict[str, Any]

JSON-serialisable {"crs", "origin_lon", "origin_lat"} payload.

Source code in packages/xrtoolz-ops/src/xrops/geom/_src/crs.py

def to_dict(self) -> dict[str, Any]:
    """JSON-serialisable ``{"crs", "origin_lon", "origin_lat"}`` payload."""
    return {
        "crs": self.crs,
        "origin_lon": float(self.origin_lon),
        "origin_lat": float(self.origin_lat),
    }

from_dict `classmethod` ¶

from_dict(d: dict[str, Any]) -> LocalFrame

Rebuild a frame from a :meth:to_dict payload.

Parameters:

Name	Type	Description	Default
`d`	`dict[str, Any]`	Mapping with `"crs"`, `"origin_lon"` and `"origin_lat"` keys; `origin_xy` is recomputed, never read.	required

Returns:

Type	Description
`LocalFrame`	The reconstructed :class:`LocalFrame`.

Source code in packages/xrtoolz-ops/src/xrops/geom/_src/crs.py

@classmethod
def from_dict(cls, d: dict[str, Any]) -> LocalFrame:
    """Rebuild a frame from a :meth:`to_dict` payload.

    Args:
        d: Mapping with ``"crs"``, ``"origin_lon"`` and ``"origin_lat"``
            keys; ``origin_xy`` is recomputed, never read.

    Returns:
        The reconstructed :class:`LocalFrame`.
    """
    return cls(
        crs=str(d["crs"]),
        origin_lon=float(d["origin_lon"]),
        origin_lat=float(d["origin_lat"]),
    )

local_frame ¶

local_frame(origin_lon: float, origin_lat: float, *, crs: str | None = None) -> LocalFrame

Build a :class:LocalFrame anchored at (origin_lon, origin_lat).

Parameters:

Name	Type	Description	Default
`origin_lon`	`float`	WGS-84 longitude of the origin, degrees east.	required
`origin_lat`	`float`	WGS-84 latitude of the origin, degrees north.	required
`crs`	`str \| None`	Projected CRS to use. `None` selects the UTM zone containing the origin via :func:`utm_crs_for`.	`None`

Returns:

Type	Description
`LocalFrame`	The frame; `frame.to_xy(origin_lon, origin_lat)` is `(0, 0)`.

Raises:

Type	Description
`ValueError`	If `crs` is `None` and `origin_lat` is outside `[-80, 84]`, or if `crs` is not a metric projected CRS.

Source code in packages/xrtoolz-ops/src/xrops/geom/_src/crs.py

def local_frame(
    origin_lon: float, origin_lat: float, *, crs: str | None = None
) -> LocalFrame:
    """Build a :class:`LocalFrame` anchored at ``(origin_lon, origin_lat)``.

    Args:
        origin_lon: WGS-84 longitude of the origin, degrees east.
        origin_lat: WGS-84 latitude of the origin, degrees north.
        crs: Projected CRS to use. ``None`` selects the UTM zone containing
            the origin via :func:`utm_crs_for`.

    Returns:
        The frame; ``frame.to_xy(origin_lon, origin_lat)`` is ``(0, 0)``.

    Raises:
        ValueError: If ``crs`` is ``None`` and ``origin_lat`` is outside
            ``[-80, 84]``, or if ``crs`` is not a metric projected CRS.
    """
    if crs is None:
        crs = utm_crs_for(origin_lon, origin_lat)
    return LocalFrame(
        crs=crs, origin_lon=float(origin_lon), origin_lat=float(origin_lat)
    )

assign_local_xy ¶

assign_local_xy(ds: Dataset, frame: LocalFrame, *, lon: str = 'lon', lat: str = 'lat', x: str = 'x', y: str = 'y') -> Dataset

Add metric x/y coordinates (metres from frame's origin).

Rectilinear grids, swaths and point tracks are supported. Because a UTM grid is not axis-aligned with lon/lat, gridded output is 2-D: for 1-D lon/lat on different dims they are broadcast onto (lat_dim, lon_dim); for 2-D lon/lat (e.g. (scanline, ground_pixel)) they share the input dims. Paired 1-D lon/lat on the same dim (along-track points, e.g. lon(obs)/lat(obs)) are transformed elementwise and yield 1-D x/y on that dim.

Parameters:

Name	Type	Description	Default
`ds`	`Dataset`	Dataset carrying geographic coordinates.	required
`frame`	`LocalFrame`	Local frame that defines the CRS and origin.	required
`lon`	`str`	Longitude coordinate name.	`'lon'`
`lat`	`str`	Latitude coordinate name.	`'lat'`
`x`	`str`	Name of the output east coordinate.	`'x'`
`y`	`str`	Name of the output north coordinate.	`'y'`

Returns:

Type	Description
`Dataset`	`ds` with `x`/`y` coordinates (`units="m"`; 2-D for grids
`Dataset`	and swaths, 1-D for point tracks) and
`Dataset`	`ds.attrs["local_frame"] = frame.to_dict()`.

Raises:

Type	Description
`ValueError`	If `lon`/`lat` are not both 1-D or both 2-D on the same dimensions.

Source code in packages/xrtoolz-ops/src/xrops/geom/_src/crs.py

def assign_local_xy(
    ds: xr.Dataset,
    frame: LocalFrame,
    *,
    lon: str = "lon",
    lat: str = "lat",
    x: str = "x",
    y: str = "y",
) -> xr.Dataset:
    """Add metric ``x``/``y`` coordinates (metres from ``frame``'s origin).

    Rectilinear grids, swaths and point tracks are supported. Because a
    UTM grid is not axis-aligned with lon/lat, gridded output is 2-D: for
    1-D ``lon``/``lat`` on *different* dims they are broadcast onto
    ``(lat_dim, lon_dim)``; for 2-D ``lon``/``lat`` (e.g.
    ``(scanline, ground_pixel)``) they share the input dims. Paired 1-D
    ``lon``/``lat`` on the *same* dim (along-track points, e.g.
    ``lon(obs)``/``lat(obs)``) are transformed elementwise and yield 1-D
    ``x``/``y`` on that dim.

    Args:
        ds: Dataset carrying geographic coordinates.
        frame: Local frame that defines the CRS and origin.
        lon: Longitude coordinate name.
        lat: Latitude coordinate name.
        x: Name of the output east coordinate.
        y: Name of the output north coordinate.

    Returns:
        ``ds`` with ``x``/``y`` coordinates (``units="m"``; 2-D for grids
        and swaths, 1-D for point tracks) and
        ``ds.attrs["local_frame"] = frame.to_dict()``.

    Raises:
        ValueError: If ``lon``/``lat`` are not both 1-D or both 2-D on the
            same dimensions.
    """
    lon_da, lat_da = ds[lon], ds[lat]
    if lon_da.ndim == 1 and lat_da.ndim == 1 and lon_da.dims == lat_da.dims:
        # Along-track points: one (lon, lat) pair per element, no meshgrid.
        lon2d, lat2d = lon_da.values, lat_da.values
        dims: tuple[str, ...] = (str(lon_da.dims[0]),)
    elif lon_da.ndim == 1 and lat_da.ndim == 1:
        lon2d, lat2d = np.meshgrid(lon_da.values, lat_da.values)
        dims = (str(lat_da.dims[0]), str(lon_da.dims[0]))
    elif lon_da.ndim == 2 and lon_da.dims == lat_da.dims:
        lon2d, lat2d = lon_da.values, lat_da.values
        dims = tuple(str(d) for d in lon_da.dims)
    else:
        raise ValueError(
            f"assign_local_xy expects {lon!r}/{lat!r} to be both 1-D or both 2-D "
            f"on the same dims, got {lon!r}: {lon_da.dims} and "
            f"{lat!r}: {lat_da.dims}."
        )
    xx, yy = frame.to_xy(lon2d, lat2d)
    out = ds.assign_coords(
        {
            x: (dims, xx, {"units": "m", "long_name": "distance east of origin"}),
            y: (dims, yy, {"units": "m", "long_name": "distance north of origin"}),
        }
    )
    return out.assign_attrs(local_frame=frame.to_dict())

Functions¶

These pure functions back the operators above; each takes xr.DataArray/xr.Dataset and a dim: argument.

assign_crs ¶

assign_crs(ds: Dataset, crs: str = 'EPSG:4326') -> Dataset

Attach a CRS to ds via :mod:rioxarray.

Parameters:

Name	Type	Description	Default
`ds`	`Dataset`	Input dataset.	required
`crs`	`str`	Any CRS specifier accepted by :class:`pyproj.CRS`.	`'EPSG:4326'`

Returns:

Type	Description
`Dataset`	Dataset with the CRS attached (`ds.rio.crs`).

Source code in packages/xrtoolz-ops/src/xrops/geom/_src/crs.py

def assign_crs(ds: xr.Dataset, crs: str = "EPSG:4326") -> xr.Dataset:
    """Attach a CRS to ``ds`` via :mod:`rioxarray`.

    Args:
        ds: Input dataset.
        crs: Any CRS specifier accepted by :class:`pyproj.CRS`.

    Returns:
        Dataset with the CRS attached (``ds.rio.crs``).
    """
    return ds.rio.write_crs(crs)

get_crs ¶

get_crs(ds: Dataset) -> CRS | None

Return ds.rio.crs if set, otherwise None.

Source code in packages/xrtoolz-ops/src/xrops/geom/_src/crs.py

def get_crs(ds: xr.Dataset) -> CRS | None:
    """Return ``ds.rio.crs`` if set, otherwise ``None``."""
    crs = ds.rio.crs
    return CRS(crs) if crs is not None else None

reproject ¶

reproject(ds: Dataset, target_crs: str, resolution: float | None = None, resampling: str = 'bilinear') -> Dataset

Reproject a raster dataset to target_crs via :mod:rioxarray.

Parameters:

Name	Type	Description	Default
`ds`	`Dataset`	Input dataset with a CRS attached.	required
`target_crs`	`str`	Any CRS specifier accepted by :class:`pyproj.CRS`.	required
`resolution`	`float \| None`	If provided, resample to this output cell size.	`None`
`resampling`	`str`	Name of the `rasterio.enums.Resampling` member.	`'bilinear'`

Returns:

Type	Description
`Dataset`	Reprojected dataset.

Source code in packages/xrtoolz-ops/src/xrops/geom/_src/crs.py

def reproject(
    ds: xr.Dataset,
    target_crs: str,
    resolution: float | None = None,
    resampling: str = "bilinear",
) -> xr.Dataset:
    """Reproject a raster dataset to ``target_crs`` via :mod:`rioxarray`.

    Args:
        ds: Input dataset with a CRS attached.
        target_crs: Any CRS specifier accepted by :class:`pyproj.CRS`.
        resolution: If provided, resample to this output cell size.
        resampling: Name of the ``rasterio.enums.Resampling`` member.

    Returns:
        Reprojected dataset.
    """
    from rasterio.enums import Resampling

    if not hasattr(Resampling, resampling):
        valid = [r.name for r in Resampling]
        raise ValueError(f"Unknown resampling {resampling!r}; expected one of {valid}.")
    return ds.rio.reproject(
        target_crs,
        resolution=resolution,
        resampling=getattr(Resampling, resampling),
    )

lonlat_to_xy ¶

lonlat_to_xy(crs: str, lon: Sequence[float] | ndarray, lat: Sequence[float] | ndarray) -> tuple[Float[ndarray, '...'], Float[ndarray, '...']]

Convert WGS-84 lon/lat to crs x/y coordinates (same shape as input).

Source code in packages/xrtoolz-ops/src/xrops/geom/_src/crs.py

def lonlat_to_xy(
    crs: str,
    lon: Sequence[float] | np.ndarray,
    lat: Sequence[float] | np.ndarray,
) -> tuple[Float[np.ndarray, "..."], Float[np.ndarray, "..."]]:
    """Convert WGS-84 lon/lat to ``crs`` x/y coordinates (same shape as input)."""
    transformer = Transformer.from_crs("EPSG:4326", crs, always_xy=True)
    x, y = transformer.transform(lon, lat)
    return np.asarray(x), np.asarray(y)

xy_to_lonlat ¶

xy_to_lonlat(crs: str, x: Sequence[float] | ndarray, y: Sequence[float] | ndarray) -> tuple[Float[ndarray, '...'], Float[ndarray, '...']]

Convert crs x/y coordinates back to WGS-84 lon/lat (same shape).

Source code in packages/xrtoolz-ops/src/xrops/geom/_src/crs.py

def xy_to_lonlat(
    crs: str,
    x: Sequence[float] | np.ndarray,
    y: Sequence[float] | np.ndarray,
) -> tuple[Float[np.ndarray, "..."], Float[np.ndarray, "..."]]:
    """Convert ``crs`` x/y coordinates back to WGS-84 lon/lat (same shape)."""
    transformer = Transformer.from_crs(crs, "EPSG:4326", always_xy=True)
    lon, lat = transformer.transform(x, y)
    return np.asarray(lon), np.asarray(lat)

CRS & Reprojection¶

Operators¶

AssignLocalXY ¶

Local metric frames¶

utm_crs_for ¶

LocalFrame dataclass ¶

to_xy ¶

to_lonlat ¶

to_dict ¶

from_dict classmethod ¶

local_frame ¶

assign_local_xy ¶

Functions¶

assign_crs ¶

get_crs ¶

reproject ¶

lonlat_to_xy ¶

xy_to_lonlat ¶

LocalFrame `dataclass` ¶

from_dict `classmethod` ¶