[gdal-dev] Problem with file source projection doing a vertical datum reprojection

[Michael Smith]

I guess I thought GTIFF_REPORT_COMPD_CS would only change the display
of the output, not change the interpretation of the SRS? I can confirm
that setting GTIFF_REPORT_COMPD_CS=FALSE does indeed allow the
vertical datum reprojection to happen but I guess I have a poor
understanding of just what GTIFF_REPORT_COMPD_CS does.

Mike


On Wed, Apr 12, 2023 at 5:24 AM Even Rouault <even.rouault at spatialys.com> wrote:
>
> Michael,
>
> it is well possible this might have worked in the past due to slightly different behavior of GDAL and/or PROJ, but the current behavior seems consistent to me. You're trying to transform between 2 compound CRS, the source one with a unknown vertical datum and the target one with EGM2008 height. There's nothing reasonable that PROJ can do regarding the vertical transformation in that situation, in particular it would be quite a bold assumption to assume that the source CRS is actually meant to be refering to ellipsoidal heights.  If you override with -s_srs EPSG:32617 (or drop REPORT_COMPD_CS=TRUE), then GDAL will auto-promote the source CRS to a 3D CRS with ellipsoidal heights, which could also be questionable from a purist point of view, but is the behaviour we have had for long now.
>
> Even
>
> Le 12/04/2023 à 10:36, Michael Smith a écrit :
>
> I have a file with a projection shown, using REPORT_COMPD_CS=TRUE of:
>
>
>
> COMPOUNDCRS["WGS 84 / UTM zone 17N + unknown",
>
>     PROJCRS["WGS 84 / UTM zone 17N",
>
>         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 17N",
>
>             METHOD["Transverse Mercator",
>
>                 ID["EPSG",9807]],
>
>             PARAMETER["Latitude of natural origin",0,
>
>                 ANGLEUNIT["degree",0.0174532925199433],
>
>                 ID["EPSG",8801]],
>
>             PARAMETER["Longitude of natural origin",-81,
>
>                 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["easting",east,
>
>                 ORDER[1],
>
>                 LENGTHUNIT["metre",1]],
>
>             AXIS["northing",north,
>
>                 ORDER[2],
>
>                 LENGTHUNIT["metre",1]],
>
>         ID["EPSG",32617]],
>
>     VERTCRS["unknown",
>
>         VDATUM["unknown"],
>
>         CS[vertical,1],
>
>             AXIS["up",up,
>
>                 LENGTHUNIT["metre",1,
>
>                     ID["EPSG",9001]]]]]
>
>
>
> When trying to gdalwarp to EGM2008 using gdalwarp -t_srs epsg:32617+3855, vertical datum reprojection does not occur unless I overwride the source srs with -s_srs epsg:32617. Using GDAL 3.6.2 and proj 9.1.1. With proj_debug, I can see that the PROJ_TRACE: vgridshift: is not occurring unless the srs is overridden. I’m fairly sure that this has worked in the past but I can’t say exactly which version. This is an older geotiff that was written about 7 years ago.
>
>
>
>
>
> --
>
> Michael Smith
>
> US Army Corps / Remote Sensing GIS Center
>
>
>
>
>
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> --
> http://www.spatialys.com
> My software is free, but my time generally not.