multi-grid GTG with overlapping grids and nodata

barry-gallagher commented 1 month ago

What is the bug?

When trying to use our proj.db and transformation grids we get all zeros or a no-op for our GTG TIFFs. It could be related to OSGeo/gdal#10395 but we don't get the vertical shift where the non-gtg tiffs give a vertical shift. Perhaps the first sub-tiff has a nodata value in the location we are querying and a zero is reported rather than inspecting the next sub-tiff. Proj.transform does return the correct transform value but gdalwarp does not.

Steps to reproduce the issue

We can supply our proj.db and tiffs if needed

Versions and provenance

Windows 10/11, gdal 3.8.4, python 3.11

Additional context

No response

rouault commented 1 month ago

Please provide elements on how to exactly reproduce the issue

barry-gallagher commented 1 month ago

Sorry for the delay. In trying to reproduce the issue I found that our initial tests were likely using the first grid in a multi-grid GTG Tiff. When I tried the data in the other ticket proj returned INF in anything other than the first grids. It seems that this is a Proj question of if our files are correct and should be supported.

@rouault do you agree it shoud be a ticket at proj? Maybe you know this shouldn't work before I make more spam over there (like should we use NTv2 instead etc). And, of course, many thanks for all your time and efforts on both projects.

echo -76.37 38.08 0 0 | cct +proj=pipeline +step +inv +proj=vgridshift "+grids=..\datum_files\non-GTG\us_noaa_nos_ITRF2014_MLLW_LMSL_(DEdelches02_vdatum_4.4_20220315_1983-2001).tif"
-76.3700000000   38.0800000000       -0.2495        0.0000

single grid transform "non-gtg"

Same location using the multi-grid file (this happens to be in the fifth grid)

echo -76.37 38.08 0 0 | cct +proj=pipeline +step +inv +proj=vgridshift "+grids=..\datum_files\GTG\us_noaa_nos_ITRF2014_MLLW_LMSL_(DEdelches02_vdatum_4.4_20220315_1983-2001).tif"
# Record 0 TRANSFORMATION ERROR: -76.37 38.08 0 0
 (Coordinate to transform falls into a grid cell that evaluates to nodata)

Pick a location for the first grid in the file instead

echo -74.5 38.0 0 0 | cct +proj=pipeline +step +inv +proj=vgridshift "+grids=..\datum_files\GTG\us_noaa_nos_ITRF2014_MLLW_LMSL_(DEdelches02_vdatum_4.4_20220315_1983-2001).tif"
-74.5000000000   38.0000000000       -0.5160        0.0000

GTG multi grid transform file

rouault commented 1 month ago

@barry-gallagher The issue is that the point -76.37 38.08 falls into several partially-overlapping grids. The current logic in PROJ is that, in that situation, it uses the first grid whose extent contains the point, independently if the grid value is nodata or not. And here the first grid that contains the point is at nodata

$ gdallocationinfo -geoloc "GTIFF_DIR:1:us_noaa_nos_ITRF2014_MLLW_LMSL_(DEdelches02_vdatum_4.4_20220315_1983-2001).tif" -76.37 38.08
Report:
  Location: (-303P,592L)

Location is off this file! No further details to report.

$ gdallocationinfo -geoloc "GTIFF_DIR:2:us_noaa_nos_ITRF2014_MLLW_LMSL_(DEdelches02_vdatum_4.4_20220315_1983-2001).tif" -76.37 38.08
Report:
  Location: (-370P,1085L)

Location is off this file! No further details to report.

$ gdallocationinfo -geoloc "GTIFF_DIR:3:us_noaa_nos_ITRF2014_MLLW_LMSL_(DEdelches02_vdatum_4.4_20220315_1983-2001).tif" -76.37 38.08
Report:
  Location: (802P,113L)
  Band 1:
    Value: -32768

$ gdallocationinfo -geoloc "GTIFF_DIR:4:us_noaa_nos_ITRF2014_MLLW_LMSL_(DEdelches02_vdatum_4.4_20220315_1983-2001).tif" -76.37 38.08
Report:
  Location: (-586P,601L)

Location is off this file! No further details to report.

$ gdallocationinfo -geoloc "GTIFF_DIR:5:us_noaa_nos_ITRF2014_MLLW_LMSL_(DEdelches02_vdatum_4.4_20220315_1983-2001).tif" -76.37 38.08
Report:
  Location: (69P,1515L)
  Band 1:
    Value: -32768

$ gdallocationinfo -geoloc "GTIFF_DIR:6:us_noaa_nos_ITRF2014_MLLW_LMSL_(DEdelches02_vdatum_4.4_20220315_1983-2001).tif" -76.37 38.08
Report:
  Location: (1172P,1625L)
  Band 1:
    Value: -0.249500006437302

$ gdallocationinfo -geoloc "GTIFF_DIR:7:us_noaa_nos_ITRF2014_MLLW_LMSL_(DEdelches02_vdatum_4.4_20220315_1983-2001).tif" -76.37 38.08
Report:
  Location: (-2524P,2411L)

Location is off this file! No further details to report.

rouault commented 1 month ago

Looking at the GGXF spec at https://raw.githubusercontent.com/opengeospatial/CRS-Gridded-Geodetic-data-eXchange-Format/master/specification/GGXF%20v1-0%20OGC-22-051r7_2024-01-08.pdf , there's a a pargraph "5.7 Multiple grids" that deals with that. We don't necessarily implement their concept of gridPriority exactly, but in any case, I don't see any mention of trying several grids when a point falls into several one and/or when its value is nodata in one of the grid candidates. The priority rules just determine the only one that is used for a given point.

barry-gallagher commented 1 month ago

Thanks for the quick response. Is there a different multi-grid format that has the behavior of checking multiple grids on nodata (perhaps NTv2)? We had read that same 5.7 section and my recollection was that nodata would check further grids. Admittedly this was a while ago I read it on a flight with a coworker. I will check with my colleague who created the files and see if it was a different spec that applied to or if we are just mistaken.

rouault commented 1 month ago

Thanks for the quick response. Is there a different multi-grid format that has the behavior of checking multiple grids on nodata (perhaps NTv2)?

no, the NTv2 reader should behave the same, and NTv2 is only for horizontal adjustment grids, not vertical ones.

We had read that same 5.7 section and my recollection was that nodata would check further grids

could be worth filing that as a question to https://github.com/opengeospatial/CRS-Gridded-Geodetic-data-eXchange-Format to see how they envision that issue, so we have possibly a common solution for that use case (although checking for pixel values when deciding the candidate sub-grid might have performance implications or implementation complication). Generally speaking nodata values in geodetic grids are a bit of a pain (especially for horizontal adjument grids, in the reverse direction). Most data producers extend the fields to avoid data holes.

OSGeo / PROJ