CRS round-trip crashes for a valid EPSG code in GDAL's PROJ DB but absent from pyproj's (e.g. EPSG:10857 / BDC Albers)

[MAfarrag]

Summary

Dataset resolves a raster's CRS with GDAL/osr but reconstructs it with pyproj, and the two ship
different PROJ databases. When a raster carries a valid custom WKT whose closest match is an EPSG code that
exists in GDAL's PROJ DB but not in pyproj's bundled one, pyramids.base.crs first identifies the code
(via osr.FindMatches) and then fails to rebuild it (via pyproj.CRS.from_user_input), raising:

pyramids.base._errors.CRSError: could not interpret 10857 as a CRS:
    Invalid projection: EPSG:10857: (Internal Proj Error: proj_create: crs not found: EPSG:10857)

This crashes any read of such a raster. The concrete trigger is the Brazil Data Cube Albers-equal-area COGs
(EPSG:10857 = SIRGAS 2000 / Brazil Albers), but the defect is generic to any code present in GDAL's PROJ
database but absent from pyproj's — it will recur every time GDAL's vendored PROJ is newer than pyproj's.

Environment (as installed)

Component	Version	Knows EPSG:10857?
GDAL (vendored in pyramids-gis)	3.13.1	✅ yes (SIRGAS 2000 / Brazil Albers)
pyproj	3.7.2	❌ no
PROJ bundled by pyproj	9.5.1	❌ no

EPSG:10857 is a real, recently-added EPSG code, not a bogus one — it simply postdates pyproj's PROJ 9.5.1
database.

Root cause

The CRS resolve path and the CRS build path use two different libraries whose PROJ databases can disagree:

1. Resolve — a raster whose WKT carries no root authority is identified through
   pyramids/base/crs.py::_epsg_from_db_match → osr.FindMatches() (GDAL's PROJ DB). For the BDC Albers WKT
   this returns EPSG:10857 at confidence 70 (_MIN_EPSG_MATCH_CONFIDENCE = 70, so it is accepted), and
   get_epsg_from_prj adopts 10857 as the raster's definitive EPSG. The match is correct — 10857 really is
   this grid.

2. Rebuild — later, the reprojection inside Dataset.crop(...) reconstructs that CRS through
   pyramids/base/crs.py::sr_from_user_input (and epsg_from_user_input), whose body is:

   # pyramids/base/crs.py  (sr_from_user_input, ~L628-635)
   try:
       wkt = CRS.from_user_input(crs).to_wkt()   # <-- pyproj, NOT osr
   except (pyproj.exceptions.CRSError, TypeError, ValueError) as exc:
       raise CRSError(f"could not interpret {crs!r} as a CRS: {exc}") from exc
   sr = osr.SpatialReference()
   sr.ImportFromWkt(wkt)

   pyproj.CRS.from_user_input(10857) calls PROJ 9.5.1's proj_create("EPSG:10857"), which does not exist
   in pyproj's database → CRSError. Note there is no int fast-path here despite the docstring example
   sr_from_user_input(3857) — every EPSG int is rebuilt through pyproj.

So pyramids resolves a code with one PROJ database and reconstructs it with another. When the databases are out
of sync, a code that was legitimately found cannot be built, and the read crashes.

The asymmetry, proven on the installed stack:

from osgeo import osr
osr.SpatialReference().ImportFromEPSG(10857)          # OK  -> "SIRGAS 2000 / Brazil Albers"
osr.SpatialReference().SetFromUserInput("EPSG:10857") # OK

import pyproj
pyproj.CRS.from_epsg(10857)
# CRSError: Invalid projection: EPSG:10857:
#   (Internal Proj Error: proj_create: crs not found: EPSG:10857)

Reproduction

A. Minimal, offline — no network, no raster (just the DB skew)

import pyproj
from osgeo import osr, gdal

# The Brazil Data Cube CBERS-4 WFI grid: a valid custom Albers WKT, NO root authority code.
wkt = (
    'PROJCS["unknown",GEOGCS["unknown",'
    'DATUM["Unknown based on GRS80 ellipsoid",'
    'SPHEROID["GRS 1980",6378137,298.257222101004,AUTHORITY["EPSG","7019"]]],'
    'PRIMEM["Greenwich",0],'
    'UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],'
    'PROJECTION["Albers_Conic_Equal_Area"],'
    'PARAMETER["latitude_of_center",-12],PARAMETER["longitude_of_center",-54],'
    'PARAMETER["standard_parallel_1",-2],PARAMETER["standard_parallel_2",-22],'
    'PARAMETER["false_easting",5000000],PARAMETER["false_northing",10000000],'
    'UNIT["metre",1,AUTHORITY["EPSG","9001"]],'
    'AXIS["Easting",EAST],AXIS["Northing",NORTH]]'
)

srs = osr.SpatialReference(wkt=wkt)
match = srs.FindMatches()                       # GDAL's PROJ DB
print(match[0][0].GetAuthorityCode(None), match[0][1])   # -> 10857  70

pyproj.CRS.from_epsg(10857)                      # pyproj's PROJ DB -> CRSError

B. Through pyramids, against real BDC data

Real asset (public, anonymous — INPE Brazil Data Cube, CBERS-4 WFI 16-day NDVI composite):

/vsicurl/https://data.inpe.br/bdc/data/cbers4-wfi-16d/v2/007/007/2024/12/18/CB4-16D_V2_007007_20241218_NDVI.tif

import pyramids  # activates vendored GDAL
from pyramids.dataset import Dataset

url = ("/vsicurl/https://data.inpe.br/bdc/data/cbers4-wfi-16d/v2/007/007/"
       "2024/12/18/CB4-16D_V2_007007_20241218_NDVI.tif")

ds = Dataset.read_file(url)
print(ds.epsg)                    # -> 10857   (resolved via GDAL FindMatches)
ds.crop(bbox=[-46.8, -23.7, -46.3, -23.2], epsg=4326, touch=True)
# pyramids.base._errors.CRSError: could not interpret 10857 as a CRS:
#   Invalid projection: EPSG:10857:
#   (Internal Proj Error: proj_create: crs not found: EPSG:10857)

The raster's WKT is a perfectly valid Albers definition — GDAL can reproject from it directly. The crash comes
only from the round-trip through an EPSG int rebuilt by a different PROJ database.

Suggested fixes (in rough order of robustness)

1. Build with the same library that resolved. In sr_from_user_input / sr_from_epsg, try
   osr.ImportFromEPSG(code) first for a bare EPSG int (it succeeds for 10857) and only fall back to pyproj —
   or, conversely, resolve and rebuild both through osr — so a code found by GDAL is always rebuildable by the
   same GDAL PROJ DB. This is the smallest, most direct fix.
2. Never discard a valid WKT for a non-round-trippable code. In get_epsg_from_prj / _epsg_from_db_match,
   before adopting a FindMatches result, verify it round-trips through the same builder used downstream; if
   it doesn't, keep and reproject from the original WKT instead of reducing to an int. This also fixes the more
   general "custom projection, no EPSG" case rather than just this one code.
3. Align the PROJ databases. Ensure pyproj uses a PROJ database at least as new as GDAL's (pin/upgrade
   pyproj's PROJ, or point PROJ_DATA at GDAL's newer proj.db), so the two never disagree. Complementary to
   (1)/(2) but does not by itself remove the cross-library fragility.

Options (1) or (2) fix the crash without a dependency bump; (2) is the most future-proof.

Impact / cross-reference

This blocks reading the entire Brazil Data Cube endpoint from the downstream earthlens STAC backend, whose
e2e-stac lane fails on exactly this trace. Filed upstream from:

• bug(stac): the Brazil Data Cube's EPSG:10857 is not a resolvable CRS earthlens#934 — bug(stac): the Brazil Data Cube's EPSG:10857 is not a resolvable CRS

Fixing this here unblocks that endpoint (and any other provider shipping COGs whose CRS is newer than pyproj's
PROJ database).