CRS & Reprojection Helpers #

CRS-handling helpers in pyramids.base.crs — the single source of truth for osr.SpatialReference construction, WKT/Proj4 → EPSG resolution, and coordinate reprojection. The most-commonly-used ones are also re-exposed as FeatureCollection static methods for ergonomic continuity (e.g. FeatureCollection.reproject_coordinates delegates to pyramids.base.crs.reproject_coordinates).

Functions #

`pyramids.base.crs.sr_from_epsg(epsg)` #

Build an :class:osr.SpatialReference from an EPSG code.

Parameters:

Name	Type	Description	Default
`epsg`	`int`	EPSG code; cast to `int` before being handed to :meth:`osr.SpatialReference.ImportFromEPSG`.	required

Returns:

Type	Description
`SpatialReference`	osr.SpatialReference: The constructed SRS.

Raises:

Type	Description
`ValueError`	If GDAL cannot resolve the EPSG code (the non-zero return path from `ImportFromEPSG` — usually propagates as a GDAL exception when `gdal.UseExceptions()` is active, which pyramids installs at package import).

Source code in src/pyramids/base/crs.py

def sr_from_epsg(epsg: int) -> osr.SpatialReference:
    """Build an :class:`osr.SpatialReference` from an EPSG code.

    Args:
        epsg: EPSG code; cast to `int` before being handed to
            :meth:`osr.SpatialReference.ImportFromEPSG`.

    Returns:
        osr.SpatialReference: The constructed SRS.

    Raises:
        ValueError: If GDAL cannot resolve the EPSG code (the
            non-zero return path from `ImportFromEPSG` — usually
            propagates as a GDAL exception when
            `gdal.UseExceptions()` is active, which pyramids
            installs at package import).
    """
    sr = osr.SpatialReference()
    err = sr.ImportFromEPSG(int(epsg))
    if err != 0:
        raise ValueError(
            f"Failed to create SRS from EPSG:{epsg} (osr returned error {err})."
        )
    return sr

`pyramids.base.crs.sr_from_wkt(wkt)` #

Build an :class:osr.SpatialReference from a WKT string.

Thin wrapper around osr.SpatialReference(wkt=wkt) that gives the WKT path a consistent name alongside :func:sr_from_epsg and :func:create_sr_from_proj. Use this when you have a WKT (the most common case in the dataset stack — dataset.crs returns WKT) and want a typed SRS without re-typing the constructor's keyword argument every call site.

Parameters:

Name	Type	Description	Default
`wkt`	`str`	Well-Known Text representation of the spatial reference.	required

Returns:

Type	Description
`SpatialReference`	osr.SpatialReference: The constructed SRS.

Examples:

Round-trip an EPSG code through WKT:

>>> from osgeo import osr
>>> from pyramids.base.crs import sr_from_epsg, sr_from_wkt
>>> wkt = sr_from_epsg(4326).ExportToWkt()
>>> sr = sr_from_wkt(wkt)
>>> sr.IsGeographic()
1

Source code in src/pyramids/base/crs.py

def sr_from_wkt(wkt: str) -> osr.SpatialReference:
    """Build an :class:`osr.SpatialReference` from a WKT string.

    Thin wrapper around `osr.SpatialReference(wkt=wkt)` that gives
    the WKT path a consistent name alongside :func:`sr_from_epsg` and
    :func:`create_sr_from_proj`. Use this when you have a WKT (the
    most common case in the dataset stack — `dataset.crs` returns
    WKT) and want a typed SRS without re-typing the constructor's
    keyword argument every call site.

    Args:
        wkt: Well-Known Text representation of the spatial reference.

    Returns:
        osr.SpatialReference: The constructed SRS.

    Examples:
        - Round-trip an EPSG code through WKT:
            ```python
            >>> from osgeo import osr
            >>> from pyramids.base.crs import sr_from_epsg, sr_from_wkt
            >>> wkt = sr_from_epsg(4326).ExportToWkt()
            >>> sr = sr_from_wkt(wkt)
            >>> sr.IsGeographic()
            1

            ```
    """
    return osr.SpatialReference(wkt=wkt)

`pyramids.base.crs.create_sr_from_proj(prj, string_type=None)` #

Create an :class:osr.SpatialReference from a projection string.

Parameters:

Name	Type	Description	Default
`prj`	`str`	The projection string (WKT, ESRI WKT, or Proj4).	required
`string_type`	`str \| None`	One of `"WKT"`, `"ESRI wkt"`, `"PROj4"`, or `None` for auto-detect (default). Auto-detect uses WKT import and falls back to ESRI WKT or Proj4 based on the prefix.	`None`

Returns:

Type	Description
`SpatialReference`	osr.SpatialReference: The constructed spatial reference.

Examples:

Parse a standard EPSG:4326 WKT string and inspect the result:

>>> from osgeo import osr
>>> ref = osr.SpatialReference()
>>> _ = ref.ImportFromEPSG(4326)
>>> wkt = ref.ExportToWkt()
>>> srs = create_sr_from_proj(wkt)
>>> srs.IsGeographic()
1
>>> srs.GetName()
'WGS 84'

Parse a Proj4 string by passing string_type="PROJ4":

>>> srs = create_sr_from_proj(
...     "+proj=longlat +datum=WGS84 +no_defs", string_type="PROJ4"
... )
>>> srs.IsGeographic()
1
>>> srs.IsProjected()
0

Parse an EPSG:3857 WKT and confirm the axis order is projected:

>>> from osgeo import osr
>>> ref = osr.SpatialReference()
>>> _ = ref.ImportFromEPSG(3857)
>>> srs = create_sr_from_proj(ref.ExportToWkt())
>>> srs.IsProjected()
1
>>> srs.GetName()
'WGS 84 / Pseudo-Mercator'

Source code in src/pyramids/base/crs.py

def create_sr_from_proj(
    prj: str, string_type: str | None = None
) -> osr.SpatialReference:
    """Create an :class:`osr.SpatialReference` from a projection string.

    Args:
        prj (str):
            The projection string (WKT, ESRI WKT, or Proj4).
        string_type (str | None):
            One of `"WKT"`, `"ESRI wkt"`, `"PROj4"`, or `None`
            for auto-detect (default). Auto-detect uses WKT import and
            falls back to ESRI WKT or Proj4 based on the prefix.

    Returns:
        osr.SpatialReference: The constructed spatial reference.

    Examples:
        - Parse a standard EPSG:4326 WKT string and inspect the result:
            ```python
            >>> from osgeo import osr
            >>> ref = osr.SpatialReference()
            >>> _ = ref.ImportFromEPSG(4326)
            >>> wkt = ref.ExportToWkt()
            >>> srs = create_sr_from_proj(wkt)
            >>> srs.IsGeographic()
            1
            >>> srs.GetName()
            'WGS 84'

            ```
        - Parse a Proj4 string by passing `string_type="PROJ4"`:
            ```python
            >>> srs = create_sr_from_proj(
            ...     "+proj=longlat +datum=WGS84 +no_defs", string_type="PROJ4"
            ... )
            >>> srs.IsGeographic()
            1
            >>> srs.IsProjected()
            0

            ```
        - Parse an EPSG:3857 WKT and confirm the axis order is projected:
            ```python
            >>> from osgeo import osr
            >>> ref = osr.SpatialReference()
            >>> _ = ref.ImportFromEPSG(3857)
            >>> srs = create_sr_from_proj(ref.ExportToWkt())
            >>> srs.IsProjected()
            1
            >>> srs.GetName()
            'WGS 84 / Pseudo-Mercator'

            ```
    """
    srs = osr.SpatialReference()
    if string_type is None:
        srs.ImportFromWkt(prj)
    elif prj.startswith(("PROJCS", "GEOGCS")):
        srs.ImportFromESRI([prj])
    else:
        srs.ImportFromProj4(prj)
    return srs

`pyramids.base.crs.get_epsg_from_prj(prj)` #

Return the EPSG code identified by a projection string.

Resolves the EPSG of the root CRS object in three steps: :meth:osr.SpatialReference.AutoIdentifyEPSG (tags recognisable CRSes), then the root AUTHORITY code, then a confident, unambiguous :meth:osr.SpatialReference.FindMatches PROJ-database lookup for well-known CRSes whose WKT lacks a root authority (e.g. a UTM PROJCS from GDAL's AAIGrid driver). The code of a child unit/datum node is never returned as if it were a CRS — that bug (issue #403) made GRIB rasters resolve to the degree-unit EPSG:9122 and UTM ASCII grids to the WGS_1984 datum EPSG:6326.

An empty input string is no longer silently mapped to 4326; that legacy default masked real configuration errors. Callers that genuinely want a fallback should handle the CRSError themselves, or use :func:epsg_from_wkt which accepts an explicit default.

Parameters:

Name	Type	Description	Default
`prj`	`str`	Projection string.	required

Returns:

Name	Type	Description
`int`	`int`	The resolved EPSG code.

Raises:

Type	Description
`CRSError`	If `prj` is an empty string, or if its root CRS carries no EPSG authority and matches no PROJ-database entry (e.g. a custom spherical-earth GRIB GEOGCS). The unit/datum codes of child nodes are never returned as a CRS.

Examples:

Resolve EPSG:4326 from its standard WKT representation:

>>> from osgeo import osr
>>> ref = osr.SpatialReference()
>>> _ = ref.ImportFromEPSG(4326)
>>> get_epsg_from_prj(ref.ExportToWkt())
4326

Resolve EPSG:3857 (Web Mercator) from its WKT representation:

>>> from osgeo import osr
>>> ref = osr.SpatialReference()
>>> _ = ref.ImportFromEPSG(3857)
>>> get_epsg_from_prj(ref.ExportToWkt())
3857

A well-known CRS whose WKT carries no root authority still resolves, via a confident PROJ-database match (here a "WGS 84 / UTM zone 18N" PROJCS with its root authority stripped):

>>> from osgeo import osr
>>> ref = osr.SpatialReference()
>>> _ = ref.ImportFromEPSG(32618)
>>> wkt = ref.ExportToWkt()
>>> wkt = wkt[: wkt.rfind(",AUTHORITY")] + "]"
>>> osr.SpatialReference(wkt=wkt).GetAuthorityCode(None) is None
True
>>> get_epsg_from_prj(wkt)
32618

A genuinely custom CRS that matches no database entry raises CRSError rather than returning a child unit/datum code (here GDAL's spherical-earth GRIB GEOGCS, whose only authority node is the degree unit EPSG:9122):

>>> grib_wkt = (
...     'GEOGCS["Coordinate System imported from GRIB file",'
...     'DATUM["unnamed",SPHEROID["Sphere",6371229,0]],'
...     'PRIMEM["Greenwich",0],'
...     'UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],'
...     'AXIS["Latitude",NORTH],AXIS["Longitude",EAST]]'
... )
>>> get_epsg_from_prj(grib_wkt)
Traceback (most recent call last):
    ...
pyramids.base._errors.CRSError: get_epsg_from_prj could not resolve an EPSG code ...

An empty projection string raises CRSError (a ValueError subclass):

>>> get_epsg_from_prj("")
Traceback (most recent call last):
    ...
pyramids.base._errors.CRSError: get_epsg_from_prj received an empty projection string. ...

Source code in src/pyramids/base/crs.py

def get_epsg_from_prj(prj: str) -> int:
    """Return the EPSG code identified by a projection string.

    Resolves the EPSG of the *root* CRS object in three steps:
    :meth:`osr.SpatialReference.AutoIdentifyEPSG` (tags recognisable
    CRSes), then the root ``AUTHORITY`` code, then a confident,
    unambiguous :meth:`osr.SpatialReference.FindMatches` PROJ-database
    lookup for well-known CRSes whose WKT lacks a root authority (e.g. a
    UTM ``PROJCS`` from GDAL's AAIGrid driver). The code of a child
    unit/datum node is never returned as if it were a CRS — that bug
    (issue #403) made GRIB rasters resolve to the degree-unit EPSG:9122
    and UTM ASCII grids to the WGS_1984 datum EPSG:6326.

    An empty input string is no longer silently mapped to `4326`; that
    legacy default masked real configuration errors. Callers that
    genuinely want a fallback should handle the `CRSError` themselves,
    or use :func:`epsg_from_wkt` which accepts an explicit `default`.

    Args:
        prj (str): Projection string.

    Returns:
        int: The resolved EPSG code.

    Raises:
        CRSError: If `prj` is an empty string, or if its root CRS carries
            no EPSG authority *and* matches no PROJ-database entry (e.g. a
            custom spherical-earth GRIB GEOGCS). The unit/datum codes of
            child nodes are never returned as a CRS.

    Examples:
        - Resolve EPSG:4326 from its standard WKT representation:
            ```python
            >>> from osgeo import osr
            >>> ref = osr.SpatialReference()
            >>> _ = ref.ImportFromEPSG(4326)
            >>> get_epsg_from_prj(ref.ExportToWkt())
            4326

            ```
        - Resolve EPSG:3857 (Web Mercator) from its WKT representation:
            ```python
            >>> from osgeo import osr
            >>> ref = osr.SpatialReference()
            >>> _ = ref.ImportFromEPSG(3857)
            >>> get_epsg_from_prj(ref.ExportToWkt())
            3857

            ```
        - A well-known CRS whose WKT carries no root authority still
          resolves, via a confident PROJ-database match (here a
          "WGS 84 / UTM zone 18N" PROJCS with its root authority stripped):
            ```python
            >>> from osgeo import osr
            >>> ref = osr.SpatialReference()
            >>> _ = ref.ImportFromEPSG(32618)
            >>> wkt = ref.ExportToWkt()
            >>> wkt = wkt[: wkt.rfind(",AUTHORITY")] + "]"
            >>> osr.SpatialReference(wkt=wkt).GetAuthorityCode(None) is None
            True
            >>> get_epsg_from_prj(wkt)
            32618

            ```
        - A genuinely custom CRS that matches no database entry raises
          `CRSError` rather than returning a child unit/datum code (here
          GDAL's spherical-earth GRIB GEOGCS, whose only authority node is
          the degree unit EPSG:9122):
            ```python
            >>> grib_wkt = (
            ...     'GEOGCS["Coordinate System imported from GRIB file",'
            ...     'DATUM["unnamed",SPHEROID["Sphere",6371229,0]],'
            ...     'PRIMEM["Greenwich",0],'
            ...     'UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],'
            ...     'AXIS["Latitude",NORTH],AXIS["Longitude",EAST]]'
            ... )
            >>> get_epsg_from_prj(grib_wkt)
            Traceback (most recent call last):
                ...
            pyramids.base._errors.CRSError: get_epsg_from_prj could not resolve an EPSG code ...

            ```
        - An empty projection string raises `CRSError` (a `ValueError` subclass):
            ```python
            >>> get_epsg_from_prj("")
            Traceback (most recent call last):
                ...
            pyramids.base._errors.CRSError: get_epsg_from_prj received an empty projection string. ...

            ```
    """
    if prj == "":
        raise CRSError(
            "get_epsg_from_prj received an empty projection string. "
            "An empty projection is ambiguous and is no longer "
            "silently defaulted to EPSG:4326. If you want "
            "a fallback EPSG, catch CRSError (also a ValueError) "
            "and supply it at the call site, or call "
            "epsg_from_wkt(prj, default=...)."
        )
    srs = create_sr_from_proj(prj)
    try:
        # AutoIdentifyEPSG attaches a root EPSG authority when it can
        # recognise the CRS; we ignore its return code and read the root
        # authority below. It raises "Unsupported SRS" for custom CRSes
        # it cannot identify (e.g. GDAL's spherical-earth GRIB GEOGCS).
        srs.AutoIdentifyEPSG()
    except RuntimeError:
        pass

    # Resolve the EPSG of the *root* CRS object only. Do NOT fall back to
    # GetAttrValue("AUTHORITY", 1): that walks the WKT tree depth-first and
    # returns the first AUTHORITY node, which for a CRS whose root carries
    # no authority is a child unit/datum code (e.g. the degree-unit
    # EPSG:9122 inside a GRIB GEOGCS, or the WGS_1984 datum EPSG:6326 inside
    # a UTM PROJCS) — a non-CRS code that breaks every downstream
    # sr_from_epsg() call. See issue #403.
    # The authority must be EPSG, not merely present. `GetAuthorityCode` returns the
    # root node's code whatever the authority is, so an ESRI-authority CRS (Robinson
    # is ESRI:54030) otherwise reports 54030 from a property called `epsg` -- a number
    # that is not an EPSG code, that pyproj cannot resolve, and that GDAL happens to
    # accept, so it round-trips on one side and fails on the other. See issue #965.
    authority = srs.GetAuthorityName(None)
    code = srs.GetAuthorityCode(None) if authority == "EPSG" else None
    if not (code and str(code).isdigit()):
        # No usable root EPSG code: either absent (AutoIdentifyEPSG could not tag
        # the root — e.g. a UTM PROJCS whose WKT lacks an AUTHORITY node), or a
        # *non-numeric* authority code — notably OGC:CRS84 (the lon/lat WGS 84 that
        # WMS/WMTS layers report), where GetAuthorityCode returns "CRS84". Try an
        # exact PROJ-database match before giving up. A UTM PROJCS resolves to its
        # numeric code here; CRS84's best match is CRS84 itself (still non-numeric),
        # so it is dropped to None below and raises CRSError — the soft
        # epsg_from_wkt() then supplies its default. The point of this branch is to
        # never crash on int("CRS84"), not to equate CRS84 with EPSG:4326.
        match = _epsg_from_db_match(srs)
        code = match if (match and str(match).isdigit()) else None
    if code is None:
        raise CRSError(
            "get_epsg_from_prj could not resolve an EPSG code from the "
            "projection: its root CRS carries no EPSG authority and matches "
            "no PROJ-database entry. This is expected for genuinely custom "
            "CRSes such as GDAL's spherical-earth GRIB GEOGCS. Catch CRSError "
            "(also a ValueError) and supply a fallback, or call "
            "epsg_from_wkt(prj, default=...)."
        )
    return int(code)

`pyramids.base.crs.epsg_from_wkt(wkt, default=4326)` #

Resolve an EPSG code from a WKT / Proj string with a fallback.

Wraps :func:get_epsg_from_prj to absorb the get_epsg_from_prj(wkt) if wkt else default idiom that was previously open-coded in four places across the dataset stack. Returns default when wkt is empty (or None), and also when get_epsg_from_prj cannot resolve an EPSG from a non-empty wkt (it raises :class:CRSError for a custom CRS whose root carries no EPSG authority — e.g. a spherical-earth GRIB GEOGCS); otherwise delegates to :func:get_epsg_from_prj.

A CRS whose authority is not EPSG — Robinson is ESRI:54030 — resolves to no EPSG code at all (issue #965), so it takes the default here just as an empty projection does. That is this function's contract, but it means the default can stand in for a real, named projection rather than only for a missing one. When that distinction matters, use :func:epsg_of_crs, which reports None, and read the CRS itself from .crs.

Use this in places where an empty projection should be treated as a soft "unknown CRS, assume WGS84" rather than a hard error — for example the Dataset.epsg property on a freshly-built in-memory raster that has no projection metadata yet. Use :func:get_epsg_from_prj directly when you want the strict behaviour where an empty projection raises.

Parameters:

Name	Type	Description	Default
`wkt`	`str \| None`	Projection string (WKT, ESRI WKT, or Proj4). An empty string or `None` returns `default`.	required
`default`	`int`	EPSG code to return when `wkt` is empty / `None`, or when its CRS cannot be resolved to an EPSG. Defaults to `4326` (the historical pyramids default).	`4326`

Returns:

Name	Type	Description
`int`	`int`	EPSG code resolved from `wkt`, or `default` when `wkt` is
	`int`	empty or its CRS carries no resolvable EPSG.

Examples:

Empty input falls back to the supplied default:

>>> from pyramids.base.crs import epsg_from_wkt
>>> epsg_from_wkt("")
4326
>>> epsg_from_wkt("", default=3857)
3857

Non-empty WKT delegates to :func:get_epsg_from_prj:

>>> from osgeo import osr
>>> from pyramids.base.crs import epsg_from_wkt
>>> ref = osr.SpatialReference()
>>> _ = ref.ImportFromEPSG(3857)
>>> epsg_from_wkt(ref.ExportToWkt())
3857

An unresolvable custom CRS falls back to default instead of raising (here GDAL's spherical-earth GRIB GEOGCS):

>>> from pyramids.base.crs import epsg_from_wkt
>>> grib_wkt = (
...     'GEOGCS["Coordinate System imported from GRIB file",'
...     'DATUM["unnamed",SPHEROID["Sphere",6371229,0]],'
...     'PRIMEM["Greenwich",0],'
...     'UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],'
...     'AXIS["Latitude",NORTH],AXIS["Longitude",EAST]]'
... )
>>> epsg_from_wkt(grib_wkt)
4326
>>> epsg_from_wkt(grib_wkt, default=3857)
3857

Source code in src/pyramids/base/crs.py

def epsg_from_wkt(wkt: str | None, default: int = 4326) -> int:
    """Resolve an EPSG code from a WKT / Proj string with a fallback.

    Wraps :func:`get_epsg_from_prj` to absorb the
    `get_epsg_from_prj(wkt) if wkt else default` idiom that was
    previously open-coded in four places across the dataset stack.
    Returns `default` when `wkt` is empty (or `None`), and also when
    `get_epsg_from_prj` cannot resolve an EPSG from a non-empty `wkt`
    (it raises :class:`CRSError` for a custom CRS whose root carries no
    EPSG authority — e.g. a spherical-earth GRIB GEOGCS); otherwise
    delegates to :func:`get_epsg_from_prj`.

    A CRS whose authority is not EPSG — Robinson is `ESRI:54030` — resolves to no
    EPSG code at all (issue #965), so it takes the `default` here just as an empty
    projection does. That is this function's contract, but it means the default can
    stand in for a real, named projection rather than only for a missing one. When
    that distinction matters, use :func:`epsg_of_crs`, which reports `None`, and read
    the CRS itself from `.crs`.

    Use this in places where an empty projection should be treated as
    a soft "unknown CRS, assume WGS84" rather than a hard error — for
    example the `Dataset.epsg` property on a freshly-built
    in-memory raster that has no projection metadata yet. Use
    :func:`get_epsg_from_prj` directly when you want the strict
    behaviour where an empty projection raises.

    Args:
        wkt: Projection string (WKT, ESRI WKT, or Proj4). An empty
            string or `None` returns `default`.
        default: EPSG code to return when `wkt` is empty / `None`, or
            when its CRS cannot be resolved to an EPSG. Defaults to
            `4326` (the historical pyramids default).

    Returns:
        int: EPSG code resolved from `wkt`, or `default` when `wkt` is
        empty or its CRS carries no resolvable EPSG.

    Examples:
        - Empty input falls back to the supplied default:
            ```python
            >>> from pyramids.base.crs import epsg_from_wkt
            >>> epsg_from_wkt("")
            4326
            >>> epsg_from_wkt("", default=3857)
            3857

            ```
        - Non-empty WKT delegates to :func:`get_epsg_from_prj`:
            ```python
            >>> from osgeo import osr
            >>> from pyramids.base.crs import epsg_from_wkt
            >>> ref = osr.SpatialReference()
            >>> _ = ref.ImportFromEPSG(3857)
            >>> epsg_from_wkt(ref.ExportToWkt())
            3857

            ```
        - An unresolvable custom CRS falls back to `default` instead of
          raising (here GDAL's spherical-earth GRIB GEOGCS):
            ```python
            >>> from pyramids.base.crs import epsg_from_wkt
            >>> grib_wkt = (
            ...     'GEOGCS["Coordinate System imported from GRIB file",'
            ...     'DATUM["unnamed",SPHEROID["Sphere",6371229,0]],'
            ...     'PRIMEM["Greenwich",0],'
            ...     'UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],'
            ...     'AXIS["Latitude",NORTH],AXIS["Longitude",EAST]]'
            ... )
            >>> epsg_from_wkt(grib_wkt)
            4326
            >>> epsg_from_wkt(grib_wkt, default=3857)
            3857

            ```
    """
    if not wkt:
        result = default
    else:
        try:
            result = get_epsg_from_prj(wkt)
        except CRSError:
            # Non-empty but unresolvable CRS (e.g. a custom spherical-earth
            # GRIB GEOGCS that carries no root EPSG authority). Treat it as
            # the same soft "unknown CRS" case as empty input rather than
            # propagating the hard error to property reads like Dataset.epsg.
            result = default
    return result

`pyramids.base.crs.epsg_of_crs(wkt)` #

Resolve the EPSG code a CRS declares, or None when it declares none.

None means "this CRS has no EPSG code", which covers two situations that :func:epsg_from_wkt conflates by substituting 4326 for both:

No CRS at all — an empty (or None) projection. The dataset is not georeferenced, so there is nothing to report. Fabricating WGS 84 here would claim a georeference the data does not have.
A CRS with no EPSG authority — a real projection the EPSG register does not name: an orthographic or geostationary projection, a rotated pole, a spherical-earth GRIB GEOGCS. Reporting 4326 for these claimed WGS 84 for grids that are not WGS 84 — an orthographic frame is not lat/lon at all, and a spherical datum differs from the WGS 84 ellipsoid by up to ~20 km.

The CRS itself is not lost in the second case: .crs still returns the WKT, and :func:crs_spec falls back to it, so reprojection and every other CRS-consuming operation keeps working. Only the code is absent, because there genuinely is not one.

This mirrors standard CRS-handling behaviour — to_epsg() returns None rather than guessing.

Parameters:

Name	Type	Description	Default
`wkt`	`str \| None`	Projection string (WKT, ESRI WKT, or Proj4), possibly empty/`None`.	required

Returns:

Type	Description
`int \| None`	int \| None: The EPSG code, or `None` when `wkt` is empty, `None`, or
`int \| None`	names a CRS that carries no EPSG authority.

Examples:

An empty projection means "no CRS", not WGS 84:

>>> from pyramids.base.crs import epsg_of_crs
>>> epsg_of_crs("") is None
True
>>> epsg_of_crs(None) is None
True

A real projection resolves to its EPSG code:

>>> from osgeo import osr
>>> from pyramids.base.crs import epsg_of_crs
>>> ref = osr.SpatialReference()
>>> _ = ref.ImportFromEPSG(3857)
>>> epsg_of_crs(ref.ExportToWkt())
3857

A CRS the EPSG register does not name has no code to report:

>>> from pyramids.base.crs import create_sr_from_proj, epsg_of_crs
>>> ortho = create_sr_from_proj(
...     "+proj=ortho +lat_0=45 +lon_0=9 +datum=WGS84 +units=m +no_defs",
...     string_type="PROJ4",
... )
>>> epsg_of_crs(ortho.ExportToWkt()) is None
True

`pyramids.base.crs.crs_spec(epsg, wkt)` #

Best usable CRS specification for a dataset, or None when it has none.

Replaces the dataset.epsg or dataset.crs idiom. That expression looks total but is not: once epsg propagates None for an ungeoreferenced raster, it evaluates to the empty CRS string, which every downstream constructor rejects with an opaque "Invalid projection: ''". Returning None instead makes the absence explicit, so callers either pass it on (producing an ungeoreferenced result) or reject it deliberately via :func:require_crs_spec.

The word usable is load-bearing, and it is why the EPSG code is not blindly preferred. Most of this function's consumers hand the result to a library that resolves it with pyproj — geopandas.GeoDataFrame.set_crs is the common one — and pyproj's bundled PROJ database is routinely older than the one GDAL vendors. An EPSG code pyramids obtained from GDAL can therefore be one pyproj cannot look up, and returning it would hand every consumer a specification that raises "crs not found" (issue #943). When that is the case and a WKT is available, the WKT is returned instead: it describes the same CRS, and pyproj parses it happily — only the catalogue lookup is missing, never the projection itself. The code is still preferred whenever it works, which is the overwhelming majority of the time and is checked once per code and cached.

Parameters:

Name	Type	Description	Default
`epsg`	`int \| None`	EPSG code, or `None` for a CRS that carries no EPSG authority.	required
`wkt`	`str \| None`	Projection WKT, or an empty string / `None` when there is no CRS.	required

Returns:

Type	Description
`int \| str \| None`	int \| str \| None: The EPSG code when there is one and it is resolvable
`int \| str \| None`	downstream, else the WKT, else the code anyway when there is no WKT to fall
`int \| str \| None`	back to, else `None`.

Examples:

An EPSG code is preferred when present:

>>> from pyramids.base.crs import crs_spec
>>> crs_spec(4326, 'GEOGCS["WGS 84"]')
4326

A CRS with no EPSG authority falls back to its WKT:

>>> from pyramids.base.crs import crs_spec
>>> crs_spec(None, 'GEOGCS["custom"]')
'GEOGCS["custom"]'

A code the downstream CRS library cannot resolve yields the WKT, so the specification stays usable:

>>> from pyramids.base.crs import crs_spec
>>> wkt = 'GEOGCS["WGS 84"]'
>>> crs_spec(999999, wkt) is wkt  # no database carries 999999
True

No CRS at all is reported as None, not as an empty string:

>>> from pyramids.base.crs import crs_spec
>>> crs_spec(None, "") is None
True

`pyramids.base.crs.require_crs_spec(epsg, wkt, operation)` #

Like :func:crs_spec, but raise when the dataset has no CRS.

Use at the point of an operation that genuinely cannot proceed without a CRS — reprojection, a coordinate transform, a spatial join against a vector. Mirrors standard CRS-handling behaviour: a missing CRS propagates quietly until something actually needs it, and then fails with a message naming the fix.

Parameters:

Name	Type	Description	Default
`epsg`	`int \| None`	EPSG code, or `None`.	required
`wkt`	`str \| None`	Projection WKT, or an empty string / `None`.	required
`operation`	`str`	Short description of what needs the CRS, used in the error.	required

Returns:

Type	Description
`int \| str`	int \| str: The EPSG code when there is one, else the WKT.

Raises:

Type	Description
`CRSError`	Neither an EPSG code nor a WKT is available.

Examples:

Resolves exactly as :func:crs_spec when a CRS is present:

>>> from pyramids.base.crs import require_crs_spec
>>> require_crs_spec(3857, "", "reproject")
3857

Refuses, naming the operation, when there is none:

>>> from pyramids.base.crs import require_crs_spec
>>> try:
...     require_crs_spec(None, "", "reproject")
... except ValueError as exc:
...     print("reproject" in str(exc))
True

Source code in src/pyramids/base/crs.py

def require_crs_spec(epsg: int | None, wkt: str | None, operation: str) -> int | str:
    """Like :func:`crs_spec`, but raise when the dataset has no CRS.

    Use at the point of an operation that genuinely cannot proceed without a
    CRS — reprojection, a coordinate transform, a spatial join against a vector.
    Mirrors standard CRS-handling behaviour: a missing
    CRS propagates quietly until something actually needs it, and then fails
    with a message naming the fix.

    Args:
        epsg: EPSG code, or `None`.
        wkt: Projection WKT, or an empty string / `None`.
        operation: Short description of what needs the CRS, used in the error.

    Returns:
        int | str: The EPSG code when there is one, else the WKT.

    Raises:
        CRSError: Neither an EPSG code nor a WKT is available.

    Examples:
        - Resolves exactly as :func:`crs_spec` when a CRS is present:
            ```python
            >>> from pyramids.base.crs import require_crs_spec
            >>> require_crs_spec(3857, "", "reproject")
            3857

            ```
        - Refuses, naming the operation, when there is none:
            ```python
            >>> from pyramids.base.crs import require_crs_spec
            >>> try:
            ...     require_crs_spec(None, "", "reproject")
            ... except ValueError as exc:
            ...     print("reproject" in str(exc))
            True

            ```
    """
    spec = crs_spec(epsg, wkt)
    if spec is None:
        raise CRSError(
            f"cannot {operation}: the raster involved has no CRS. Set one first "
            f"(e.g. "
            "`dataset.epsg = <code>`, or `gdal_edit.py -a_srs EPSG:<code> "
            "<file>` on disk); pyramids does not assume WGS 84 for an "
            "ungeoreferenced raster."
        )
    return spec

`pyramids.base.crs.cf_geographic_wkt(units, axis_units=None)` #

WGS 84 WKT when CF axis units describe a lat/lon grid, else "".

CF-1.x lets a data variable carry no grid_mapping; when its coordinate axes are in degrees east/north the file is geographic and CF simply leaves the datum implicit. GDAL reports an empty projection for those, and the whole ecosystem reads them as WGS 84 — so inferring EPSG:4326 from this evidence is a convention-backed reading of the metadata, not the blanket "assume WGS 84 for anything unprojected" default that ARC-26 removed. A raster with no CRS and no such evidence still reports no CRS.

Parameters:

Name	Type	Description	Default
`units`	`set[str]`	Lower-cased unit strings from every coordinate array, including the 2-D auxiliary lat/lon a curvilinear grid uses.	required
`axis_units`	`set[str] \| None`	Lower-cased unit strings from the true horizontal dimension axes only. A unit here that belongs to a projected or rotated frame (`m`, `km`, unqualified `degrees` / `degree`, `rad`, ...) vetoes the inference: the grid is projected, rotated-pole or geostationary, and its degrees arrays are auxiliary coordinates. See :data:`PROJECTED_AXIS_UNITS` for the exact set. Defaults to `None` (no veto).	`None`

Returns:

Name	Type	Description
`str`	`str`	WGS 84 WKT when both a longitude and a latitude axis are in degrees,
	`str`	otherwise `""`.

Examples:

Degrees on both axes identify a geographic grid:

>>> from pyramids.base.crs import cf_geographic_wkt
>>> wkt = cf_geographic_wkt({"degrees_east", "degrees_north"})
>>> "WGS 84" in wkt
True

The CF singular spellings are accepted too:

>>> from pyramids.base.crs import cf_geographic_wkt
>>> bool(cf_geographic_wkt({"degree_east", "degree_north"}))
True

One axis alone, or non-degree units, is not evidence:

>>> from pyramids.base.crs import cf_geographic_wkt
>>> cf_geographic_wkt({"degrees_east", "meter"})
''

Source code in src/pyramids/base/crs.py

def cf_geographic_wkt(units: set[str], axis_units: set[str] | None = None) -> str:
    """WGS 84 WKT when CF axis units describe a lat/lon grid, else ``""``.

    CF-1.x lets a data variable carry no ``grid_mapping``; when its coordinate
    axes are in degrees east/north the file *is* geographic and CF simply leaves
    the datum implicit. GDAL reports an empty projection for those, and the whole
    ecosystem reads them as WGS 84 — so inferring EPSG:4326 from this evidence is
    a convention-backed reading of the metadata, not the blanket "assume WGS 84
    for anything unprojected" default that ARC-26 removed. A raster with no CRS
    and no such evidence still reports no CRS.

    Args:
        units: Lower-cased unit strings from every coordinate array, including
            the 2-D auxiliary lat/lon a curvilinear grid uses.
        axis_units: Lower-cased unit strings from the true *horizontal* dimension
            axes only. A unit here that belongs to a projected or rotated frame
            (`m`, `km`, unqualified `degrees` / `degree`, `rad`, ...) vetoes the
            inference: the grid is projected, rotated-pole or geostationary, and
            its degrees arrays are auxiliary coordinates. See
            :data:`PROJECTED_AXIS_UNITS` for the exact set. Defaults to `None`
            (no veto).

    Returns:
        str: WGS 84 WKT when both a longitude and a latitude axis are in degrees,
        otherwise ``""``.

    Examples:
        - Degrees on both axes identify a geographic grid:
            ```python
            >>> from pyramids.base.crs import cf_geographic_wkt
            >>> wkt = cf_geographic_wkt({"degrees_east", "degrees_north"})
            >>> "WGS 84" in wkt
            True

            ```
        - The CF singular spellings are accepted too:
            ```python
            >>> from pyramids.base.crs import cf_geographic_wkt
            >>> bool(cf_geographic_wkt({"degree_east", "degree_north"}))
            True

            ```
        - One axis alone, or non-degree units, is not evidence:
            ```python
            >>> from pyramids.base.crs import cf_geographic_wkt
            >>> cf_geographic_wkt({"degrees_east", "meter"})
            ''

            ```
    """
    has_lon = any(u.startswith(LON_UNIT_PREFIXES) for u in units)
    has_lat = any(u.startswith(LAT_UNIT_PREFIXES) for u in units)
    # A linear unit on a real *axis* means the grid is projected, and any degrees
    # arrays are auxiliary lat/lon coordinates rather than the CRS. Without this
    # a CF file with metre x/y plus 2-D aux lat/lon and no grid_mapping is
    # reported as WGS 84 on a metre geotransform. The check looks only at
    # `axis_units` because a data variable may legitimately be in metres (a ROMS
    # bathymetry or sea-surface height) on an otherwise geographic grid.
    projected = any(u in PROJECTED_AXIS_UNITS for u in (axis_units or set()))
    geographic = has_lon and has_lat and not projected
    return sr_from_epsg(4326).ExportToWkt() if geographic else ""

`pyramids.base.crs.reproject_coordinates(x, y, *, from_crs=4326, to_crs=3857, precision=6)` #

Reproject parallel x / y coordinate lists between CRSes.

Argument and return order is (x, y) throughout; accepts any CRS form :meth:pyproj.Transformer.from_crs understands (EPSG int, EPSG string, WKT, Proj4, :class:pyproj.CRS).

Parameters:

Name	Type	Description	Default
`x`	`list[float]`	X-coordinates in the source CRS (longitudes when `from_crs` is geographic).	required
`y`	`list[float]`	Y-coordinates in the source CRS (latitudes when `from_crs` is geographic).	required
`from_crs`	`Any`	Source CRS. Accepts anything :meth:`pyproj.Transformer.from_crs` accepts: EPSG integer (`4326`), authority string (`"EPSG:4326"`), WKT, Proj4, or a :class:`pyproj.CRS` instance. Default `4326`.	`4326`
`to_crs`	`Any`	Target CRS, same forms as `from_crs`. Default `3857`.	`3857`
`precision`	`int \| None`	Decimal places to round each returned coordinate to, using Python's built-in `round` — correctly rounded in decimal, which is not the same as `numpy.round` on values that are not exactly representable (`round(2.675, 2)` is `2.67`, `numpy.round(2.675, 2)` is `2.68`). Pass `None` to disable rounding and get the transformer's full output. Default `6`.	`6`

Returns:

Type	Description
`tuple[list[float], list[float]]`	tuple[list[float], list[float]]: `(x, y)` in the target CRS.

Raises:

Type Description

ValueError

If len(x)!= len(y).

CRSError

If :meth:pyproj.Transformer.from_crs raises one of pyproj.exceptions.CRSError (malformed WKT / proj string), TypeError (input is not CRS-like — e.g. a bare object()), or ValueError (out-of-range EPSG integer). The wrapper converts each into pyramids' :class:pyramids.base._errors.CRSError so callers do not need to import pyproj to catch bad-CRS failures, and the message names both CRSes plus the underlying explanation. Other exception types (AttributeError, ImportError, …) propagate unchanged — they signal a real bug, not a bad user input.

Examples:

Reproject a WGS84 point into Web Mercator:

>>> from pyramids.base.crs import reproject_coordinates
>>> x, y = reproject_coordinates(
...     [31.0], [30.0], from_crs=4326, to_crs=3857
... )
>>> round(x[0])
3450904
>>> round(y[0])
3503550

Source code in src/pyramids/base/crs.py

def reproject_coordinates(
    x: list[float],
    y: list[float],
    *,
    from_crs: Any = 4326,
    to_crs: Any = 3857,
    precision: int | None = 6,
) -> tuple[list[float], list[float]]:
    """Reproject parallel x / y coordinate lists between CRSes.

    Argument and return order is `(x, y)` throughout; accepts any
    CRS form :meth:`pyproj.Transformer.from_crs` understands (EPSG
    int, EPSG string, WKT, Proj4, :class:`pyproj.CRS`).

    Args:
        x (list[float]):
            X-coordinates in the source CRS (longitudes when
            `from_crs` is geographic).
        y (list[float]):
            Y-coordinates in the source CRS (latitudes when
            `from_crs` is geographic).
        from_crs:
            Source CRS. Accepts anything
            :meth:`pyproj.Transformer.from_crs` accepts: EPSG integer
            (`4326`), authority string (`"EPSG:4326"`), WKT, Proj4,
            or a :class:`pyproj.CRS` instance. Default `4326`.
        to_crs:
            Target CRS, same forms as `from_crs`. Default `3857`.
        precision (int | None):
            Decimal places to round each returned coordinate to, using
            Python's built-in `round` — correctly rounded in decimal,
            which is *not* the same as `numpy.round` on values that are
            not exactly representable (`round(2.675, 2)` is `2.67`,
            `numpy.round(2.675, 2)` is `2.68`). Pass `None` to disable
            rounding and get the transformer's full output. Default `6`.

    Returns:
        tuple[list[float], list[float]]: `(x, y)` in the target CRS.

    Raises:
        ValueError: If `len(x)!= len(y)`.
        CRSError: If :meth:`pyproj.Transformer.from_crs` raises one
            of `pyproj.exceptions.CRSError` (malformed WKT / proj
            string), `TypeError` (input is not CRS-like — e.g. a
            bare `object()`), or `ValueError` (out-of-range EPSG
            integer). The wrapper converts each into pyramids'
            :class:`pyramids.base._errors.CRSError` so callers do not
            need to import pyproj to catch bad-CRS failures, and the
            message names both CRSes plus the underlying explanation.
            Other exception types (`AttributeError`, `ImportError`,
            …) propagate unchanged — they signal a real bug, not a bad
            user input.

    Examples:
        - Reproject a WGS84 point into Web Mercator:
            ```python
            >>> from pyramids.base.crs import reproject_coordinates
            >>> x, y = reproject_coordinates(
            ...     [31.0], [30.0], from_crs=4326, to_crs=3857
            ... )
            >>> round(x[0])
            3450904
            >>> round(y[0])
            3503550

            ```
    """
    if len(x) != len(y):
        raise ValueError(
            f"x and y must have equal length; got len(x)={len(x)} vs. len(y)={len(y)}."
        )
    try:
        # Through `crs_from_user_input`, not `Transformer.from_crs` directly, so a code
        # only GDAL's PROJ database knows still builds a transformer. See issue #943.
        transformer = Transformer.from_crs(
            crs_from_user_input(from_crs), crs_from_user_input(to_crs), always_xy=True
        )
    except (pyproj.exceptions.CRSError, TypeError, ValueError) as exc:
        raise CRSError(
            f"reproject_coordinates failed to parse CRS "
            f"(from_crs={from_crs!r}, to_crs={to_crs!r}): {exc}"
        ) from exc
    # One vectorized call over the whole arrays rather than one call per point:
    # `Transformer.transform` accepts array input and does the loop in PROJ, so a
    # polygon ring with thousands of vertices costs one Python call, not thousands.
    xs, ys = transformer.transform(
        np.asarray(x, dtype=float), np.asarray(y, dtype=float)
    )
    out_x = np.asarray(xs, dtype=float).tolist()
    out_y = np.asarray(ys, dtype=float).tolist()
    if precision is not None:
        # Round with the built-in, NOT `np.round`. They are not interchangeable:
        # `round` is correctly rounded in decimal, while `np.round` scales by
        # `10**precision`, rounds, and divides back, so the two disagree on values
        # that are not exactly representable -- `round(2.675, 2)` is `2.67` but
        # `np.round(2.675, 2)` is `2.68`, and at the default `precision=6` they
        # differ on roughly 1 in 3000 Web-Mercator-magnitude coordinates. Keeping
        # the built-in preserves the per-point implementation's output exactly;
        # the expensive part was the PROJ round trip, which is already vectorized.
        out_x = [round(value, precision) for value in out_x]
        out_y = [round(value, precision) for value in out_y]
    return out_x, out_y