feat: Add user-defined Oblique Stereographic CRS fixture

.github/workflows/ci.yml

@@ -18,5 +18,8 @@ jobs:
       - name: Check generated files are up to date
         run: pixi run check
 
+      - name: Run tests
+        run: pixi run test
+
       # Use ruff-action so we get annotations in the Github UI
       - uses: astral-sh/ruff-action@v4.1.0

pixi.toml

@@ -10,6 +10,7 @@ generate = "python generate_fixtures.py"
 info = "bash scripts/info.sh"
 check = "bash scripts/check.sh"
 generate-npy = "python scripts/generate_npy.py"
+test = "pytest tests"
 
 [dependencies]
 gdal = ">=3.12.1,<4"
@@ -19,6 +20,8 @@ python = ">=3.12"
 rasterio = ">=1.5.0,<2"
 tifffile = ">=2024.1.30"
 affine = ">=2.4.0,<3"
+pytest = ">=8"
 
 [pypi-dependencies]
 rio-cogeo = ">=7.0.1, <8"
+async-tiff = ">=0.7.2, <0.8"

real_data/source-coop-dataforcanada/O2308000_7586000_cog_info.md

@@ -0,0 +1,61 @@
+```
+Driver: GTiff
+File: real_data/source-coop-dataforcanada/O2308000_7586000_cog.tif
+COG: True
+Compression: LZW
+ColorSpace: None
+
+Profile
+    Width:            128
+    Height:           128
+    Bands:            3
+    Tiled:            False
+    Dtype:            uint8
+    NoData:           None
+    Alpha Band:       False
+    Internal Mask:    False
+    Interleave:       PIXEL
+    ColorMap:         False
+    ColorInterp:      ('red', 'green', 'blue')
+    Scales:           (1.0, 1.0, 1.0)
+    Offsets:          (0.0, 0.0, 0.0)
+
+Geo
+    Crs:              EPSG:2953
+    Origin:           (2307999.0, 7588002.0)
+    Resolution:       (0.1, -0.1)
+    BoundingBox:      (2307999.0, 7587989.2, 2308011.8, 7588002.0)
+    MinZoom:          20
+    MaxZoom:          20
+
+Image Metadata
+    DataType: Generic
+    AREA_OR_POINT: Area
+
+Image Structure
+    LAYOUT: COG
+    COMPRESSION: LZW
+    INTERLEAVE: PIXEL
+
+Band 1
+    ColorInterp: red
+    Metadata:
+        BandName: Band_1
+        RepresentationType: ATHEMATIC
+
+Band 2
+    ColorInterp: green
+    Metadata:
+        BandName: Band_2
+        RepresentationType: ATHEMATIC
+
+Band 3
+    ColorInterp: blue
+    Metadata:
+        BandName: Band_3
+        RepresentationType: ATHEMATIC
+
+IFD
+    Id      Size           BlockSize     Decimation
+    0       128x128        128x128       0
+```

real_data/source-coop-dataforcanada/README.md

@@ -0,0 +1,54 @@
+Includes data from
+
+https://source.coop/dataforcanada
+
+## User-defined Oblique Stereographic CRS
+
+New Brunswick Stereographic, equivalent to EPSG:2953, but stored as a
+user-defined projected CRS (`ProjectedCSTypeGeoKey = 32767`) with
+`ProjCoordTransGeoKey = 16` (`CT_ObliqueStereographic`). The projection origin
+is stored in `ProjNatOriginLatGeoKey`, `ProjNatOriginLongGeoKey` and
+`ProjScaleAtNatOriginGeoKey`, not the `ProjCenter*` keys.
+
+See https://github.com/source-cooperative/cog-viewer/issues/40.
+
+`listgeo -no_norm O2308000_7586000_cog.tif`:
+
+```
+   Keyed_Information:
+      GTModelTypeGeoKey (Short,1): ModelTypeProjected
+      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
+      GTCitationGeoKey (Ascii,42): "NAD83(CSRS) / New Brunswick Stereographic"
+      GeographicTypeGeoKey (Short,1): User-Defined
+      GeogCitationGeoKey (Ascii,43): "GCS Name = NAD83(CSRS)|Primem = Greenwich|"
+      GeogGeodeticDatumGeoKey (Short,1): Code-6140 (NAD83 Canadian Spatial Reference System)
+      GeogAngularUnitsGeoKey (Short,1): Angular_Degree
+      GeogSemiMajorAxisGeoKey (Double,1): 6378137
+      GeogInvFlatteningGeoKey (Double,1): 298.257222101
+      GeogPrimeMeridianLongGeoKey (Double,1): 0
+      ProjectedCSTypeGeoKey (Short,1): User-Defined
+      ProjectionGeoKey (Short,1): User-Defined
+      ProjCoordTransGeoKey (Short,1): CT_ObliqueStereographic
+      ProjLinearUnitsGeoKey (Short,1): Linear_Meter
+      ProjNatOriginLongGeoKey (Double,1): -66.5
+      ProjNatOriginLatGeoKey (Double,1): 46.5
+      ProjFalseEastingGeoKey (Double,1): 2500000
+      ProjFalseNorthingGeoKey (Double,1): 7500000
+      ProjScaleAtNatOriginGeoKey (Double,1): 0.999912
+      End_Of_Keys.
+```
+
+GDAL identifies the source CRS as EPSG:2953 by name when reading, so a plain
+`gdal_translate` writes `ProjectedCSTypeGeoKey = 2953` instead. Assigning the
+same WKT without the `PROJCS` and `GEOGCS` authority codes makes GDAL write the
+same user-defined geo keys as the source file:
+
+```bash
+pixi run gdal_translate \
+  -srcwin 0 0 128 128 \
+  -a_srs 'PROJCS["NAD83(CSRS) / New Brunswick Stereographic",GEOGCS["NAD83(CSRS)",DATUM["NAD83_Canadian_Spatial_Reference_System",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6140"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Oblique_Stereographic"],PARAMETER["latitude_of_origin",46.5],PARAMETER["central_meridian",-66.5],PARAMETER["scale_factor",0.999912],PARAMETER["false_easting",2500000],PARAMETER["false_northing",7500000],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]' \
+  -of COG \
+  -co BLOCKSIZE=128 \
+  /vsicurl/https://data.source.coop/dataforcanada/test/O2308000_7586000_cog.tif \
+  O2308000_7586000_cog.tif
+```

tests/test_real_data.py

@@ -0,0 +1,58 @@
+"""Check the raw GeoKeys of real_data fixtures.
+
+Some fixtures exist to exercise a specific GeoKey layout, which regenerating
+them with GDAL can silently change: GDAL may identify a user-defined CRS as an
+EPSG code and write only that code. Reading the raw keys with async-tiff
+catches that.
+"""
+
+import asyncio
+from pathlib import Path
+from typing import Any
+
+from async_tiff import TIFF
+from async_tiff.store import LocalStore
+
+REPO_ROOT = Path(__file__).parent.parent
+
+
+def read_geo_keys(path: str) -> dict[str, Any]:
+    """Read the GeoKeys of the first IFD of a TIFF, relative to the repo root."""
+
+    async def _read() -> dict[str, Any]:
+        tiff = await TIFF.open(path, store=LocalStore(REPO_ROOT))
+        gkd = tiff.ifds[0].geo_key_directory
+        return {key: gkd[key] for key in gkd.keys()}
+
+    return asyncio.run(_read())
+
+
+def test_dataforcanada_user_defined_oblique_stereographic():
+    """Assert that CRS is stored in custom geo keys, not simply as `EPSG:2953`.
+
+    See real_data/source-coop-dataforcanada/README.md.
+    """
+    geo_keys = read_geo_keys(
+        "real_data/source-coop-dataforcanada/O2308000_7586000_cog.tif"
+    )
+    assert geo_keys == {
+        "model_type": 1,
+        "raster_type": 1,
+        "citation": "NAD83(CSRS) / New Brunswick Stereographic",
+        "geographic_type": 32767,
+        "geog_citation": "GCS Name = NAD83(CSRS)|Primem = Greenwich|",
+        "geog_geodetic_datum": 6140,
+        "geog_angular_units": 9102,
+        "geog_semi_major_axis": 6378137.0,
+        "geog_inv_flattening": 298.257222101,
+        "geog_prime_meridian_long": 0.0,
+        "projected_type": 32767,
+        "projection": 32767,
+        "proj_coord_trans": 16,
+        "proj_linear_units": 9001,
+        "proj_nat_origin_long": -66.5,
+        "proj_nat_origin_lat": 46.5,
+        "proj_false_easting": 2500000.0,
+        "proj_false_northing": 7500000.0,
+        "proj_scale_at_nat_origin": 0.999912,
+    }