Four ways your GIS stack loses a planetary CRS

IAU_2015 is a versioned authority, and most readers drop the version

python
planetary
geospatial
gdal
proj
rasterio
crs
Published

2026-09-02

If you work with Earth data, a CRS round-trip is boring: you write EPSG:32618, you read it back, to_authority() returns ('EPSG', '32618'), and you never think about it again.

Planetary data is not boring. PROJ ships the IAU coordinate reference systems under the authority name IAU_2015 — a name that carries an edition year, because the IAU revises body parameters. Mars alone has

IAU_2015:49900   Mars (2015) - Sphere / Ocentric     a = b = 3396190, f = 0
IAU_2015:49901   Mars (2015) / Ographic              a = 3396190, b = 3376200
IAU_2015:49910   ... / Equirectangular, clon = 0

That versioned authority is the problem. Below are four round-trips through common formats. One of them works.

One thing to say up front, because it frames everything below. The IAU’s own work here is long-settled; what is new is the plumbing — PROJ only gained the IAU authority in version 8.2, and GeoParquet’s CRS metadata is younger still. Every library in this post is implementing a young standard, and most of them are getting most of it right. So read this as a status check, not a complaint: a snapshot of where planetary CRS identification can be relied on today, and where it cannot. I would expect several of these rows to change, and I would be glad if they did.

All output is from a single script (given at the end) on:

rasterio 1.5.1 | GDAL 3.13.2 | PROJ 9.8.1 | pyproj 3.7.2 | geopandas 1.1.4

The starting point: nothing is wrong yet

from pyproj import CRS
crs = CRS.from_user_input("IAU_2015:49910")
PROJJSON id    : {"authority": "IAU", "code": 49910, "version": 2015}
to_authority() : ('IAU_2015', '49910')

PROJJSON has three slots — authority, code, version — and all three are filled. Note the authority is stored as IAU with a separate version: 2015; the string IAU_2015 is a rendering of that pair, not what is stored.

Keep that structure in mind. Every failure below is a format or a reader with nowhere to put the third slot.

The test files

The four round-trips below all start from that same crs object and write it into three files — a raster, and one shared GeoDataFrame through two vector drivers, GeoPackage (.gpkg) and GeoParquet (.parquet):

import geopandas as gpd
import numpy as np
import rasterio
from rasterio.transform import from_origin
from shapely.geometry import box

with rasterio.open("x.tif", "w", driver="GTiff", height=8, width=8, count=1,
                   dtype="uint8", crs=crs,
                   transform=from_origin(1e5, 1e5, 10, 10)) as dst:
    dst.write(np.zeros((8, 8), "uint8"), 1)

gdf = gpd.GeoDataFrame({"a": [1]}, geometry=[box(0, 0, 1e4, 1e4)], crs=crs)
gdf.to_file("x.gpkg", driver="GPKG", layer="t")
gdf.to_parquet("x.parquet")

Eight-by-eight zeros and a single square. Nothing about the pixels or the geometry matters — the files exist only to carry the CRS from a writer to a reader. x.gpkg and x.parquet get theirs from the same gdf, which is what makes sections 3 and 4 a fair comparison: same CRS, same object, different driver.

1. GeoTIFF through GDAL: depends which accessor you use

dataset = gdal.Open("x.tif")
wkt1 = dataset.GetProjection()                                     # WKT1
wkt2 = dataset.GetSpatialRef().ExportToWkt(["FORMAT=WKT2_2019"])   # WKT2
GetProjection() tail        : NORTH],AUTHORITY["IAU","49910"]]
-> to_authority()           : None
GetSpatialRef() WKT2 tail   : tre",1]],ID["IAU",49910,2015]]
-> to_authority()           : ('IAU_2015', '49910')
GetAuthorityName/Code       : IAU_2015 49910

The file is fine. WKT2 carries ID["IAU",49910,2015] — authority, code, and version — and GetAuthorityName() returns the full IAU_2015.

WKT1 cannot. Its AUTHORITY["name","code"] node has exactly two slots and no third one to put the edition in, so AUTHORITY["IAU","49910"] is the most it can say. Nothing downstream can recover 2015 from that string, because it was never in it.

This is not a GDAL bug. It is a WKT1 format limitation, and GDAL documents GetProjectionRef() as “a compatibility layer around GetSpatialRef()” since GDAL 3.0. RFC 73 exists precisely to move callers off representation-dependent strings:

using “a SRS that is not dependent of the representation (WKT 1 or WKT 2), hence using a OGRSpatialReference object instead of a const char* string”

2. The same GeoTIFF through rasterio: lost

rio_crs = rasterio.open("x.tif").crs
crs.to_authority()          : None
crs.to_wkt() tail           : NORTH],AUTHORITY["IAU","49910"]]
to_wkt(version='WKT2_2019') : T["metre",1]],ID["IAU",49910]]
-> to_authority()           : None

Two things worth separating here.

rasterio reads the CRS through the WKT1 accessor. In rasterio/_base.pyx:

def read_crs(self):
    cdef const char *wkt = GDALGetProjectionRef(self.handle())   # WKT1
    return self._handle_crswkt(wkt)

So the version is destroyed before rasterio’s CRS object exists. GDALGetSpatialRef — the RFC 73 replacement — is not bound in rasterio/gdal.pxi at all, which suggests this is a GDAL-2-era carry-over rather than a decision.

Asking for WKT2 does not help, and by the paragraph above it cannot. to_wkt(version= "WKT2_2019") emits ID["IAU",49910] — WKT2 without the version field — because the CRS object it is exporting never had one. The loss already happened at read time; no choice of output format undoes it. It is worth showing only because requesting WKT2 is the obvious thing to reach for, and it hands back a well-formed WKT2 string that looks like it worked.

(Careful with the call: to_wkt’s first positional parameter is morph_to_esri_dialect, not the version. to_wkt("WKT2_2019") silently warns and returns WKT1.)

rasterio’s CRS class is perfectly capable — CRS.from_wkt(wkt2) on GDAL’s WKT2 string gives ('IAU_2015', '49910'), and so does CRS.from_user_input("IAU_2015:49910"). Only the file reading path loses it.

3. GeoPackage: the identifier survives, the version does not

This is the interesting case, because the format does better than the reader.

A GeoPackage is a SQLite database with a spec-mandated table gpkg_spatial_ref_sys, and each layer, each geometry column, and each geometry BLOB header points into it by an integer srs_id. The row has dedicated authority columns:

import sqlite3

con = sqlite3.connect("x.gpkg")
org, code, wkt = next(con.execute(
    "select organization, organization_coordsys_id, definition "
    "from gpkg_spatial_ref_sys where organization not in ('NONE','EPSG')"))

The filter skips the three rows the spec mandates in every GeoPackage — two NONE placeholders and EPSG:4326 — leaving the layer’s actual CRS:

org                         : 'IAU'
code                        : 49910
wkt tail                    : NORTH],AUTHORITY["IAU","49910"]]
CRS.from_authority(org, code)
-> to_authority()           : ('IAU_2015', '49910')   <-- resolves!
gpd.read_file("x.gpkg").crs
-> to_authority()           : None                    <-- but this is what you get

organization + organization_coordsys_id resolves — but it is worth being precise about why, because the reason is not that the file kept enough.

There is no authority called IAU in PROJ. proj.db’s authority_list holds EPSG, ESRI, IAU_2015, IGNF, NKG, NRCAN, OGC, PROJ, and every IAU system sits under IAU_2015. The bare name resolves through a dedicated mapping table, which today holds exactly one row:

-- versioned_auth_name_mapping
versioned_auth_name   auth_name   version   priority
'IAU_2015'            'IAU'       '2015'    1

So CRS.from_authority('IAU', 49910) does not recover 2015 from the GeoPackage. It supplies it locally, from a lookup table with one candidate in it. What the file actually says is IAU:49910 — ambiguous in principle, unambiguous today only because there is nothing yet to be ambiguous with. PROJ plainly expects that to change: priority, and the UNIQUE (auth_name, priority) constraint alongside it, exist to rank editions of the same authority. Ship an IAU_2025 and a bare IAU reference is settled by PROJ’s ranking rather than by what the file meant.

The reader does not consult those columns regardless — it reconstructs from the WKT definition, which is WKT1, which cannot carry the version. So the GeoPackage row is better than what the reader does with it, but it is not the lossless record it looks like.

GeoPackage 1.2 added a gpkg_crs_wkt extension with a definition_12_063 column for WKT2. Enabling it does not help either:

gdf.to_file("x2.gpkg", driver="GPKG", CRS_WKT_EXTENSION="YES")
with CRS_WKT_EXTENSION=YES  : "metre",1]],ID["IAU",49910]]
-> to_authority()           : None

The WKT2 written into the extension column is missing the version, even though GDAL’s GeoTIFF path writes ID["IAU",49910,2015] from the same CRS. So this one is a writer defect, not a format limitation.

4. GeoParquet: works

No special reader needed — just open it and ask the CRS what it is:

crs_back = gpd.read_parquet("x.parquet").crs
projjson = crs_back.to_json_dict()
crs_back.to_authority()      : ('IAU_2015', '49910')
projjson["id"]               : {"authority": "IAU", "code": 49910, "version": 2015}
projjson["base_crs"]["id"]   : {"authority": "IAU", "code": 49900, "version": 2015}

GeoParquet stores the CRS as PROJJSON, and PROJJSON has a first-class version field. Nothing has to be squeezed into a two-slot node or inferred from a string. It round-trips losslessly, and it is the only one of the four that does. The lineage survives too — the base geographic CRS carries its own versioned identifier, so the sphere the projection is built on is named as explicitly as the projection itself.

to_json_dict() is worth keeping in your pocket, because it is also the sharpest way to see what the other formats lost. The same two lines against the GeoPackage:

gpd.read_file("x.gpkg").crs.to_json_dict()["id"]
{"authority": "IAU", "code": 49910}

Same authority, same code, no version key at all. That is the whole difference between sections 3 and 4 in one line: not a different identifier, just a missing third slot — and with it gone, to_authority() has nothing to return but None.

It is not the authority name

A reasonable guess is that IAU versus IAU_2015 is the mismatch — the files say IAU, PROJ’s authority is called IAU_2015, so lookups fail. That guess is wrong:

CRS.from_authority('IAU', 49910)      -> ('IAU_2015', '49910')
CRS.from_authority('IAU_2015', 49910) -> ('IAU_2015', '49910')

PROJ accepts either — the mapping table in section 3 is what makes the short name work. So the name is not the problem. The problem is exclusively that readers reconstruct the CRS from a WKT string when a better reference is at hand: losslessly through a non-string API (GDAL’s GetSpatialRef), or, in GeoPackage, through authority columns that at least outlive the WKT1 the reader chooses to parse.

Summary

path carries the version? to_authority()
in memory (pyproj) yes, PROJJSON ('IAU_2015', '49910')
GeoTIFF → GDAL GetSpatialRef() WKT2 yes ('IAU_2015', '49910')
GeoTIFF → GDAL GetProjection() WKT1 no — format limit None
GeoTIFF → rasterio no — reader uses WKT1 None
GeoPackage → geopandas/GDAL no — authority columns only, unused None
GeoParquet → geopandas yes, PROJJSON ('IAU_2015', '49910')

Three distinct causes, worth keeping apart:

  1. A format limitation. WKT1 has no version slot. Nothing to fix; just do not round-trip through it.
  2. A reader choosing the lossy path. rasterio reads WKT1 when GetSpatialRef exists; the GeoPackage reader parses WKT when authority columns are right there.
  3. A writer dropping a field it has. GDAL’s GeoPackage WKT2 omits the version its GeoTIFF WKT2 includes.

Only (1) is inherent.

What to do meanwhile

  • Verify with GDAL’s authority accessors, not to_authority(): gdal.Open(f).GetSpatialRef().GetAuthorityName(None). On correctly-labelled files, rasterio ... .crs.to_authority() returning None means nothing about the data.
  • Prefer GeoParquet for tabular vector data that must carry a planetary CRS faithfully; keep GeoPackage as the QGIS-facing export.
  • For GeoPackage, read the authority columns directly:
org, code = con.execute("""
    select srs.organization, srs.organization_coordsys_id
    from gpkg_geometry_columns geom
    join gpkg_spatial_ref_sys srs using (srs_id)
    where geom.table_name = ?""", (layer,)).fetchone()
crs = CRS.from_authority(org, code)

This leans on PROJ resolving a bare IAU to its one installed edition. It is the best the format gives you, but record the edition yourself if you need it to stay pinned.

  • Do not conclude a file has a broken CRS because a reader says None. Every file in this post was written correctly.

This matters beyond tidiness. Mars has both a sphere (49900) and an ellipsoid (49901) in the same authority, differing by 20 km in polar radius, and planetary products are routinely unlabelled or ESRI-flavoured. When identification silently degrades to None, the usual next step is someone guessing — and a guess between 49900 and 49901 is a real geometric error, not a metadata nit.

The reproducer

import sqlite3, tempfile, warnings
warnings.filterwarnings("ignore")

import geopandas as gpd
import numpy as np
import rasterio
from osgeo import gdal
from pyproj import CRS
from rasterio.transform import from_origin
from shapely.geometry import box

gdal.UseExceptions()
crs = CRS.from_user_input("IAU_2015:49910")
tmpdir = tempfile.mkdtemp()
tif, gpkg, pqt = f"{tmpdir}/x.tif", f"{tmpdir}/x.gpkg", f"{tmpdir}/x.parquet"

with rasterio.open(tif, "w", driver="GTiff", height=8, width=8, count=1,
                   dtype="uint8", crs=crs,
                   transform=from_origin(1e5, 1e5, 10, 10)) as dst:
    dst.write(np.zeros((8, 8), "uint8"), 1)

gdf = gpd.GeoDataFrame({"a": [1]}, geometry=[box(0, 0, 1e4, 1e4)], crs=crs)
gdf.to_file(gpkg, driver="GPKG", layer="t")
gdf.to_parquet(pqt)

dataset = gdal.Open(tif)
wkt1 = dataset.GetProjection()
wkt2 = dataset.GetSpatialRef().ExportToWkt(["FORMAT=WKT2_2019"])
print("GDAL WKT1 :", CRS.from_wkt(wkt1).to_authority())
print("GDAL WKT2 :", CRS.from_wkt(wkt2).to_authority())
print("rasterio  :", rasterio.open(tif).crs.to_authority())
print("rio WKT2  :", CRS.from_wkt(
    rasterio.open(tif).crs.to_wkt(version="WKT2_2019")).to_authority())

con = sqlite3.connect(gpkg)
org, code = next(con.execute(
    "select organization, organization_coordsys_id from gpkg_spatial_ref_sys "
    "where organization not in ('NONE','EPSG')"))
print("GPKG cols :", CRS.from_authority(org, code).to_authority())
print("GPKG read :", gpd.read_file(gpkg).crs.to_authority())

gpq_crs = gpd.read_parquet(pqt).crs
print("GeoParquet:", gpq_crs.to_json_dict()["id"])
print("GPQ read  :", gpq_crs.to_authority())