For "local UTM" output CRS, write "backup" 326XX EPSG code in geotiff header

[dshean]

Issue encountered during Monday tag-up when downloading the PCD sample geotiff products. A sample *mos.tif geotiff was downloaded, but I forgot to also download the sidecar .aux.xml containing the WKT2 3D CRS definition, which cannot currently be stored in the geotiff header (see uw-cryo/3D_CRS_Transformation_Resources#7 and opengeospatial/geotiff#56).

gdalinfo *mos.tif
Driver: GTiff/GeoTIFF
Files: AZ_PimaCo_2_2021-DTM_fill_window_size_4_mos.tif
Size is 7335, 13497
Origin = (491902.000000000000000,3534741.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  491902.000, 3534741.000)
Lower Left  (  491902.000, 3521244.000)
Upper Right (  499237.000, 3534741.000)
Lower Right (  499237.000, 3521244.000)
Center      (  495569.500, 3527992.500)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
  NoData Value=-9999
  Overviews: 3668x6749, 1834x3375, 917x1688, 459x844, 230x422, 115x211

The geotransform is properly set, but there is no CRS information to interpret the origin coordinates. When the file was added into QGIS project with a different project CRS (a different projected CRS), the geotiff was rendered in the wrong location without any warning!

A potential hack to avoid this with lidar_tools is to do all processing, write out the aux.xml, but then set the CRS of the output geotiff using the corresponding UTM CRS using the WGS84 ensemble (EPSG:326XX). This could be accomplished with a final gdal_edit -a_srs EPSG:326XX *mos.tif call or a few lines with GDAL API.

If the sidecar xml with the full 3D CRS definition using specific WGS84 realization is present, that is read and used...

gdalinfo *mos.tif
Driver: GTiff/GeoTIFF
Files: AZ_PimaCo_2_2021-DTM_fill_window_size_4_mos.tif
       AZ_PimaCo_2_2021-DTM_fill_window_size_4_mos.tif.aux.xml
Size is 7335, 13497
Coordinate System is:
PROJCRS["WGS 84 (G2139) / UTM 12N",
    BASEGEOGCRS["WGS 84 (G2139)",
        DYNAMIC[
            FRAMEEPOCH[2016]],
        DATUM["World Geodetic System 1984 (G2139)",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]]],
    CONVERSION["UTM zone 12N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-111,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,3],
        AXIS["(E)",north,
            MERIDIAN[90,
                ANGLEUNIT["degree",0.0174532925199433]],
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            MERIDIAN[0,
                ANGLEUNIT["degree",0.0174532925199433]],
            ORDER[2],
            LENGTHUNIT["metre",1]],
        AXIS["ellipsoidal height (h)",up,
            ORDER[3],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
Data axis to CRS axis mapping: 1,2,3
Origin = (491902.000000000000000,3534741.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  491902.000, 3534741.000) (111d 5' 8.47"W, 31d56'54.94"N)
Lower Left  (  491902.000, 3521244.000) (111d 5' 8.06"W, 31d49'36.57"N)
Upper Right (  499237.000, 3534741.000) (111d 0'29.06"W, 31d56'55.04"N)
Lower Right (  499237.000, 3521244.000) (111d 0'29.03"W, 31d49'36.67"N)
Center      (  495569.500, 3527992.500) (111d 2'48.66"W, 31d53'15.83"N)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
  Min=671.114 Max=1120.684
  NoData Value=-9999
  Overviews: 3668x6749, 1834x3375, 917x1688, 459x844, 230x422, 115x211

But if not, at least we have something reasonable in the tif

gdalinfo test_32612.tif
Driver: GTiff/GeoTIFF
Files: test_32612.tif
Size is 7335, 13497
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 12N",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 12N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-111,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Navigation and medium accuracy spatial referencing."],
        AREA["Between 114°W and 108°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - Alberta; Northwest Territories (NWT); Nunavut; Saskatchewan. Mexico. United States (USA)."],
        BBOX[0,-114,84,-108]],
    ID["EPSG",32612]]
Data axis to CRS axis mapping: 1,2
Origin = (491902.000000000000000,3534741.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  491902.000, 3534741.000) (111d 5' 8.47"W, 31d56'54.94"N)
Lower Left  (  491902.000, 3521244.000) (111d 5' 8.06"W, 31d49'36.57"N)
Upper Right (  499237.000, 3534741.000) (111d 0'29.06"W, 31d56'55.04"N)
Lower Right (  499237.000, 3521244.000) (111d 0'29.03"W, 31d49'36.67"N)
Center      (  495569.500, 3527992.500) (111d 2'48.66"W, 31d53'15.83"N)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
  NoData Value=-9999
  Overviews: 3668x6749, 1834x3375, 917x1688, 459x844, 230x422, 115x211

This is unsatisfying and is a long-term solution, but it potentially avoids some inevitable user confusion/questions (I know people will forget to preserve the aux.xml files, as we've all done this several times, potentially leading to hours of unnecessary troublshooting down the line).

Because of the ~2.0 m uncertainty of the WGS84 ensemble, PROJ (and QGIS) will not actually show a difference between the geolocation of these two files when loaded in the same project.

[dshean]

The alternative is to drop the Ellipsoid height definition for the compound 3D CRS definition, as Even suggests: https://lists.osgeo.org/pipermail/gdal-dev/2025-January/060027.html

And to just be clear in the metadata that all heights are relative to the WGS84(G2139) ellipsoid.

[scottyhq]

The alternative is to drop the Ellipsoid height definition for the compound 3D CRS definition, as Even suggests: https://lists.osgeo.org/pipermail/gdal-dev/2025-January/060027.html And to just be clear in the metadata that all heights are relative to the WGS84(G2139) ellipsoid.

I think I favor this because dealing with sidecar files or storing conflicting information (ensemble + specific realization) is more confusing. Agreed that embedding a comment as a geotiff tag would be good and also documenting promotion to a 3D CRS for proper warping if needed.