unknown CRS in geotiff causes errors while transform_bounds

[JesJehle]

To get the bounds during upload we use:

# Open with rasterio and get bounds
with rasterio.open(str(target_path), 'r') as src:
    bounds = src.bounds
    transformed_bounds = rasterio.warp.transform_bounds(src.crs, 'EPSG:4326', *bounds)

I cases where CRS is unknows like in spain_14_09_2023_sierra_de_baza_2_ortho.tif

Driver: GTiff/GeoTIFF
Warning 1: The definition of projected CRS EPSG:3857 got from GeoTIFF keys is not the same as the one from the EPSG registry, which may cause issues during reprojection operations. Set GTIFF_SRS_SOURCE configuration option to EPSG to use official parameters (overriding the ones from GeoTIFF keys), or to GEOKEYS to use custom values from GeoTIFF keys and drop the EPSG code.
Files: 00ef0583-7bf2-49f4-88eb-3607a8bceabb_spain_14_09_2023_sierra_de_baza_2_ortho.tif
Size is 35064, 12549
Coordinate System is:
ENGCRS["WGS 84 / Pseudo-Mercator",
    EDATUM["Unknown engineering datum"],
    CS[Cartesian,2],
        AXIS["easting",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    ID["EPSG",3857]]
Data axis to CRS axis mapping: 1,2
Origin = (-321403.103884378680959,4499542.093287520110607)
Pixel Size = (0.037308065466493,-0.037466616489492)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( -321403.104, 4499542.093) 
Lower Left  ( -321403.104, 4499071.925) 
Upper Right ( -320094.934, 4499542.093) 
Lower Right ( -320094.934, 4499071.925) 
Center      ( -320749.019, 4499307.009) 
Band 1 Block=35064x32 Type=Byte, ColorInterp=Red
  Overviews: 17532x6275, 8766x3138, 4383x1569, 2192x785, 1096x393, 548x197, 274x99, 137x50
  Unit Type: metre
Band 2 Block=35064x32 Type=Byte, ColorInterp=Green
  Overviews: 17532x6275, 8766x3138, 4383x1569, 2192x785, 1096x393, 548x197, 274x99, 137x50
  Unit Type: metre
Band 3 Block=35064x32 Type=Byte, ColorInterp=Blue
  Overviews: 17532x6275, 8766x3138, 4383x1569, 2192x785, 1096x393, 548x197, 274x99, 137x50
  Unit Type: metre

the following error occurs:

api-1      |   File "rasterio/_warp.pyx", line 1560, in rasterio._warp._transform_bounds
api-1      |   File "rasterio/_err.pyx", line 359, in rasterio._err.exc_wrap_pointer
api-1      | rasterio._err.CPLE_NotSupportedError: Cannot find coordinate operations from 'ENGCRS["WGS 84 / Pseudo-Mercator",EDATUM[""],CS[Cartesian,2],AXIS["easting",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["northing",north,ORDER[2],LENGTHUNIT["metre",1]],ID["EPSG",3857]]' to 'EPSG:4326'

@cmosig how to handle such cases?

[cmosig]

This is weird. QGIS is able to handle this CRS just fine.

[cmosig]

https://github.com/rasterio/rasterio/issues/1563 does the transform work maybe with higher GDAL versions?

[JesJehle]

this seems to do the trick:

with Env(GTIFF_SRS_SOURCE='EPSG'):
    with rasterio.open(str(target_path), 'r') as src:
        try:
                       transformed_bounds = rasterio.warp.transform_bounds(src.crs, 'EPSG:4326', *src.bounds)
        except Exception as e:
               logger.error(f'No CRS found for {target_path}: {e}')
               return HTTPException(status_code=400, detail=f'No CRS found for {target_path}: {e}')

Setting GTIFF_SRS_SOURCE='EPSG': This configuration option forces GDAL to use the official EPSG parameters for the CRS, which can resolve inconsistencies between the CRS defined in the GeoTIFF keys and the EPSG registry.

[cmosig]

Okay, that is very interesting.

[JesJehle]

I also run into the error in processing the cog: Warning 1: The definition of projected CRS EPSG:3857 got from GeoTIFF keys is not the same as the one from the EPSG registry, which may cause issues during reprojection operations. Set GTIFF_SRS_SOURCE configuration option to EPSG to use official parameters (overriding the ones from GeoTIFF keys), or to GEOKEYS to use custom values from GeoTIFF keys and drop the EPSG code.

This is interesting: https://gis.stackexchange.com/questions/439555/crs-registry-definition-issues-epsg3035

[JesJehle]

ENGCRS["WGS 84 / Pseudo-Mercator",EDATUM[""],CS[Cartesian,2],AXIS["easting",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["northing",north,ORDER[2],LENGTHUNIT["metre",1]],ID["EPSG",3857]]' to `EPSG:3857

Error: ERROR 6: Cannot find coordinate operations from ENGCRS["WGS 84 / Pseudo-Mercator",EDATUM[""],CS[Cartesian,2],AXIS["easting",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["northing",north,ORDER[2],LENGTHUNIT["metre",1]],ID["EPSG",3857]]' toEPSG:3857'

[cmosig]

This is interesting: https://gis.stackexchange.com/questions/439555/crs-registry-definition-issues-epsg3035

Hmm but this is a different EPSG. Maybe still related though

[cmosig]

I wonder how we ended up with this weird ENGCRS as this was produced by us an an export from Agisoft Metashape, which is what everyone is using.

[cmosig]

If we want to make our lives easy, we could simply reject non PROJ / EPSG CRS definitions and have the user to the proper reprojection. But this would need to be made clear to the user in the error message.

It's weird however that QGIS is able to load this data just fine, but rasterio not. Both use GDAL in the background I'd expect

[JesJehle]

the issue is that the registry entry of the crs is not the same as the official registry. gdalinfo of the file actually explain it all:

Warning 1: The definition of projected CRS EPSG:3857 got from GeoTIFF keys is not the same as the one from the EPSG registry, which may cause issues during reprojection operations. Set GTIFF_SRS_SOURCE configuration option to EPSG to use official parameters (overriding the ones from GeoTIFF keys), or to GEOKEYS to use custom values from GeoTIFF keys and drop the EPSG code.

if we believe the local entry is wrong, set GTIFF_SRS_SOURCE=EPSG like I did with:

with Env(GTIFF_SRS_SOURCE='EPSG'):

Otherwise, we use a custom CRS based on the GEOKEYS I would strongly prefer the first version.

[cmosig]

Hmm, let's do that, but can we somehow add a log message or similar if this happens, so that the data audit person can visually revalidate the georeferencing?