Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

GDAL Warp behavior when facing NoDataValue in the multigrid-GTG transformation grid #11650

Closed
mohammadashkezari-noaa opened this issue Jan 14, 2025 · 23 comments · Fixed by OSGeo/PROJ#4382
Closed

GDAL Warp behavior when facing NoDataValue in the multigrid-GTG transformation grid #11650

mohammadashkezari-noaa opened this issue Jan 14, 2025 · 23 comments · Fixed by OSGeo/PROJ#4382
Assignees
@rouault

Comments

@mohammadashkezari-noaa
Copy link

mohammadashkezari-noaa commented Jan 14, 2025 •
edited
Loading

What is the bug?

We are using GDAL Warp and custom GTG grids to apply vertical transformations on raster files. The elevation values are transformed as expected except over the regions where the underlying transformation grid has NoDataValues, in which case the transformed values remain identical to the input. Ideally, the transformed values should be NoDataValue when the underlying transformation grid is NoDataValue. Is it be possible to get this behavior?

Thanks.

01_Input
02_Output
03_Grid
04_Grid_and_Output

Steps to reproduce the issue

The input and output (transformed) raster files, our custom GTG grid, and the updated proj.db can be found here:

https://www.dropbox.com/scl/fo/3hg6szwq3gjywvj05daun/AHDL2ooWDQqCFOF6QU_7SiI?rlkey=vd3qexqtiw4gpm3e1eixfbpeq&st=vwmx3mfp&dl=0

We use the following warp options (in python) to apply the vertical shifts:

  warp_kwargs = {
                 "outputType": gdal.gdalconst.GDT_Float32,
                 "srcBands": [1],
                 "dstBands": [1],
                 "warpOptions": ["APPLY_VERTICAL_SHIFT=YES",
                                 "SAMPLE_GRID=YES",
                                 "SAMPLE_STEPS=ALL"
                                 ],
                 "errorThreshold": 0,
                 }

And this is the essence of the function that applies the transformation:

def warp(input_file: str,
         output_file: str,
         apply_vertical: bool,
         crs_from: Union[pp.CRS, str],
         crs_to: Union[pp.CRS, str],
         input_metadata: dict,
         warp_kwargs: Optional[dict] = None
         ):
    """
    A gdal-warp wrapper to transform an NBS raster (GTiff) file with 3 bands:
    Elevation, Uncertainty, and Contributors.

    Parameters
    ----------
    input_file: str
        Path to the input raster file (gdal supported).
    output_file: str
        Path to the transformed raster file.
    apply_vertical: bool
        Apply GDAL vertical shift.
    crs_from: pyproj.crs.CRS or input used to create one
        Projection of input data.
    crs_to: pyproj.crs.CRS or input used to create one
        Projection of output data.
    input_metadata: dict
        Dictionary containing metadata generated by vyperdatum.raster_utils.raster_metadata()
    warp_kwargs: dict
        gdal kwargs.

    Returns
    --------
    str:
        Absolute path to the transformed file.
    """
    if isinstance(crs_from, pp.CRS):
        crs_from = crs_to_code_auth(crs_from)
    if isinstance(crs_to, pp.CRS):
        crs_to = crs_to_code_auth(crs_to)

    gdal.Warp(destNameOrDestDS=output_file,
              srcDSOrSrcDSTab=input_file,
              dstSRS=crs_to,
              srcSRS=crs_from,
              **(warp_kwargs or {})
              )

    if apply_vertical:
        if isinstance(warp_kwargs.get("srcBands"), list):
            ds_in = gdal.Open(input_file, gdal.GA_ReadOnly)
            ds_out = gdal.Open(output_file, gdal.GA_ReadOnly)
            # combine the vertically transformed bands and
            # the non-transformed ones into a new raster
            driver = gdal.GetDriverByName(input_metadata["driver"])
            mem_path = f"/vsimem/{os.path.splitext(os.path.basename(output_file))[0]}.tiff"
            ds_temp = driver.Create(mem_path,
                                    ds_in.RasterXSize,
                                    ds_in.RasterYSize,
                                    ds_in.RasterCount,
                                    gdal.GDT_Float32
                                    )
            ds_temp.SetGeoTransform(ds_out.GetGeoTransform())
            ds_temp.SetProjection(ds_out.GetProjection())
            for b in range(1, ds_in.RasterCount+1):
                if b in warp_kwargs.get("srcBands"):
                    out_shape = ds_out.GetRasterBand(b).ReadAsArray().shape
                    in_shape = ds_in.GetRasterBand(b).ReadAsArray().shape
                    if out_shape != in_shape:
                        logger.error(f"Band {b} dimensions has changed from"
                                     f"{in_shape} to {out_shape}")
                    ds_temp.GetRasterBand(b).WriteArray(ds_out.GetRasterBand(b).ReadAsArray())
                    ds_temp.GetRasterBand(b).SetDescription(ds_out.GetRasterBand(b).GetDescription())
                    ds_temp.GetRasterBand(b).SetNoDataValue(ds_out.GetRasterBand(b).GetNoDataValue())
                else:
                    ds_temp.GetRasterBand(b).WriteArray(ds_in.GetRasterBand(b).ReadAsArray())
                    ds_temp.GetRasterBand(b).SetDescription(ds_in.GetRasterBand(b).GetDescription())
                    ds_temp.GetRasterBand(b).SetNoDataValue(ds_in.GetRasterBand(b).GetNoDataValue())
            ds_in, ds_out = None, None
            ds_temp.FlushCache()
            driver.CreateCopy(output_file, ds_temp)
            ds_temp = None
            gdal.Unlink(mem_path)

    return output_file

Versions and provenance

GDLA 3.8.4

pyproj info:
    pyproj: 3.6.1
      PROJ: 9.4.0
PROJ DATA (recommended version): 1.19
PROJ Database: 1.4
EPSG Database: v11.016 [2024-08-31]
ESRI Database: ArcGIS Pro 3.3 [2024-08-16]
IGNF Database: 3.1.0 [2019-05-24]

System:
    python: 3.11.11 | packaged by conda-forge | (main, Dec  5 2024, 14:06:23) [MSC v.1942 64 bit (AMD64)]
executable: C:\Users\mohammad.ashkezari\.conda\envs\proj94\python.exe
   machine: Windows-10-10.0.19045-SP0

Python deps:
   certifi: 2024.12.14
    Cython: None
setuptools: None
       pip: 24.3.1

Additional context

No response
@jratike80
Copy link
Collaborator

I would not call the current behavior (to keep the original heights) as a bug.

I would post-process the result for example with gdal_calc https://gdal.org/en/stable/programs/gdal_calc.html#gdal-calc. Copy the pixels from the target image, except write nodata where the transformation grid has nodata. Alternatively update the target image with gdal_rasterize https://gdal.org/en/stable/programs/gdal_rasterize.html. Convert the nodata areas first into polygons and then burn nodata value into the pixels which are covered by the polygons.

Questions should go to the gdal-dev mailing list at https://lists.osgeo.org/mailman/listinfo/gdal-dev or other support forums such as https://gis.stackexchange.com/. GitHub issues are for bug reports and suggestions for new features.

@GlenRice-NOAA
Copy link
Contributor

The GTG behavior is expected to be different from single resolution transformation grids? A transformation to a new file with single resolution transformation grid results in no data in the areas where there is no transformation. This is why we thought this might be a bug. The user may not be aware of the type of transformation grid used by a pipeline and thus would not expect differing results requiring further postprocessing.

@jratike80
Copy link
Collaborator

The GTG behavior is expected to be different from single resolution transformation grids?

Sorry, I notice that I wrote my opinions about a thing that I do not understand. Especially I did not notice the different behavior with multiresolution vs. single resolution grids. Reopening.

If I read the documentation https://proj.org/en/stable/specifications/geodetictiffgrids.html correctly, Proj should always use the highest resolution grids

If a low-resolution grid is available, it should be put before subgrids of higher-resolution in the chain of IFD linking. On reading, PROJ will use the value from the highest-resolution grid that contains the point of interest.

@jratike80 jratike80 reopened this Jan 14, 2025
@GlenRice-NOAA
Copy link
Contributor

No worries! We appreciate your engagement.

The note from PROJ is what made us think that if a nodata value is returned from the transformation GTG the resulting transformation should also be nodata for a new file, similar to a single resolution grid.

@jratike80
Copy link
Collaborator

I may not understand the tiffinfo output right, but maybe the higher-resolution grids are first, not last. What do you think?
tiffinfo us_noaa_nos_MLLW-ITRF2020_2020.0_(nwldatum_4.7.0_20240621).tif

@barry-gallagher
Copy link

barry-gallagher commented Jan 14, 2025 •
edited
Loading

Hello, Agreed, the highest res grids are used and we wrote a pre-processor that takes data from the lower res and inserts it to the higher res as needed. That part works fine.

This issue should be independent of the resolution order since there is a nodata value in all three underlying grids and yet we still get an output value (like -39 from the original raster) at those positions unlike if the grid_transform file is a single res tiff then we get a nodata value at that position.

@rouault
Copy link
Member

rouault commented Jan 14, 2025 •
edited
Loading

More a PROJ issue (or feature) than a GDAL one. There's a strong ambiguity in the community of use of PROJ:

  • some users don't know that PROJ exist at all and hardly know anything about CRS, and hardly care about "correct" results. They might perfer "something possibly incorrectly located" to "nothing if unsure". This is especially true if they use it for example in QGIS at scales where datum shifts don't matter that much.
  • picky geodesists that prefer exactness.

PROJ defaults tend to serve more the former than the later. So what you get by default is:

$ echo 34.528078 -77.355236 0| PROJ_DATA=/tmp PROJ_DEBUG=2 cs2cs -d 8 "ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height" ITRF2020
pj_open_lib(proj.ini): call fopen(/tmp/proj.ini) - failed
pj_open_lib(proj.db): call fopen(/tmp/proj.db) - succeeded
pj_open_lib(us_noaa_nos_MLLW-ITRF2020_2020.0_(nwldatum_4.7.0_20240621).tif): call fopen(/tmp/us_noaa_nos_MLLW-ITRF2020_2020.0_(nwldatum_4.7.0_20240621).tif) - succeeded
pj_open_lib(us_noaa_nos_MLLW-ITRF2020_2020.0_(nwldatum_4.7.0_20240621).tif): call fopen(/tmp/us_noaa_nos_MLLW-ITRF2020_2020.0_(nwldatum_4.7.0_20240621).tif) - succeeded
Using coordinate operation Inverse of 'ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020'
pj_open_lib(us_noaa_nos_MLLW-ITRF2020_2020.0_(nwldatum_4.7.0_20240621).tif): call fopen(/tmp/us_noaa_nos_MLLW-ITRF2020_2020.0_(nwldatum_4.7.0_20240621).tif) - succeeded
Partially intersecting grids found!
[... snip ...]
Partially intersecting grids found!
Attempt to use coordinate operation Inverse of 'ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020' failed. This might become an error in a future PROJ major release. Set the ONLY_BEST option to YES or NO. This warning will no longer be emitted (for the current transformation instance).
Coordinate to transform falls into a grid cell that evaluates to nodata
Did not result in valid result. Attempting a retry with another operation.
Using coordinate operation Transformation from MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 (ballpark vertical transformation, without ellipsoid height to vertical height correction)
34.52807800	-77.35523600 0.00000000

PROJ does notice that the grid evaluates to nodata, and thus fallbacks to using the ballpark operation that does not vertical shift.

you can likely get the result you want by setting the PROJ_ONLY_BEST_DEFAULT=YES environment variable (environment variable, not GDAL configuration option), or by setting "only_best_default=on" in PROJ's proj.ini file (cf https://proj.org/en/stable/resource_files.html#proj-ini)

$ echo 34.528078 -77.355236 0| PROJ_ONLY_BEST_DEFAULT=YES PROJ_DATA=/tmp PROJ_DEBUG=2 bin/cs2cs "ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height" ITRF2020
pj_open_lib(proj.ini): call fopen(/tmp/proj.ini) - failed
pj_open_lib(proj.db): call fopen(/tmp/proj.db) - succeeded
pj_open_lib(us_noaa_nos_MLLW-ITRF2020_2020.0_(nwldatum_4.7.0_20240621).tif): call fopen(/tmp/us_noaa_nos_MLLW-ITRF2020_2020.0_(nwldatum_4.7.0_20240621).tif) - succeeded
pj_open_lib(us_noaa_nos_MLLW-ITRF2020_2020.0_(nwldatum_4.7.0_20240621).tif): call fopen(/tmp/us_noaa_nos_MLLW-ITRF2020_2020.0_(nwldatum_4.7.0_20240621).tif) - succeeded
Using coordinate operation Inverse of 'ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020'
pj_open_lib(us_noaa_nos_MLLW-ITRF2020_2020.0_(nwldatum_4.7.0_20240621).tif): call fopen(/tmp/us_noaa_nos_MLLW-ITRF2020_2020.0_(nwldatum_4.7.0_20240621).tif) - succeeded
Partially intersecting grids found!
[ ... snip ... ]
Partially intersecting grids found!
Attempt to use coordinate operation Inverse of 'ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020' failed.
*	* inf

I didn't check how gdalwarp behaved given your reproducer lacks a bit the exact argument used to invoke it, but hopefully that should help it. You might need to specify the output bounds and resolution ,because in that strict mode, it can struggle getting valid sample points to guess it !

@rouault rouault self-assigned this Jan 14, 2025
@mohammadashkezari-noaa
Copy link
Author

Thank you.
Setting os.environ["PROJ_ONLY_BEST_DEFAULT"] = "YES" and Warp's xRes, yRes, outputBounds params didn't change the warp's final output.
Although, with PROJ_ONLY_BEST_DEFAULT=YES in place, GDAL now explicitly raises -many- PROJ errors that read like:

ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.

@rouault
Copy link
Member

rouault commented Jan 14, 2025

Setting os.environ["PROJ_ONLY_BEST_DEFAULT"] = "YES"

as far as I know, that is too late. You need to set it before starting python

@mohammadashkezari-noaa
Copy link
Author

Put it at the system level

> echo $Env:PROJ_ONLY_BEST_DEFAULT
YES

Still getting the same warp output and the same PROJ error messages as above.

@rouault
Copy link
Member

rouault commented Jan 14, 2025

Still getting the same warp output and the same PROJ error messages as above.

please provide a gdalwarp command line to reproduce the issue

@jratike80
Copy link
Collaborator

jratike80 commented Jan 14, 2025 •
edited
Loading

I made an experiment by:

  • write the proj.db and the gridshift file from the test data into my proj_data directory
  • do set PROJ_ONLY_BEST_DEFAULT=YES
  • run gdalwarp -wo apply_vertical_shif=yes -wo sample_grid=yes -wo sample_steps=all input.tif test.tif --debug on

I thing the test was not successful because of
PROJ: proj_as_wkt: C:\OSGeo4W\share\proj\proj.db contains DATABASE.LAYOUT.VERSION.MINOR = 3 whereas a number >= 4 is expected. It comes from another PROJ installation.

I won't have a try with older proj because I do not know if the test makes sense. And with the knowledge that I have it would be faster for me to burn the nodata polygons into the result.

@jratike80
Copy link
Collaborator

Two questions:

  • Could it be better to define the transformation with a Proj pipeline instead of using a tailored proj.db? It might be easier to understand and control what happens.
  • If the single resolution grid file gives a desired result, could you share that as well?

@mohammadashkezari-noaa
Copy link
Author

@rouault & @jratike80 Thanks for your help.
For simplicity I down-sampled and cropped the input file (see small_input.tif).
I ran the following command both usig proj9.4 and 9.5:
gdalwarp -s_srs EPSG:9989 -t_srs EPSG:9990+NOAA:98 -wo apply_vertical_shift=yes -wo sample_grid=yes -wo sample_steps=all small_input.tif output_proj94.tif --debug on

I also wrote an equivalent python script to the warp operation. They all appear to behave the same. Please see the attachment.

gdalwarp

nodata.zip

@mohammadashkezari-noaa
Copy link
Author

BTW, here is the gdalwarp stdout:

> gdalwarp -s_srs EPSG:9989 -t_srs EPSG:9990+NOAA:98 -wo apply_vertical_shift=yes -wo sample_grid=yes -wo sample_steps=all small_input.tif output_proj94.tif --debug on

GDAL: GDALOpen(small_input.tif, this=0000018CB27CAF40) succeeds as GTiff.
GDAL: GDALOpen(output_proj94.tif, this=0000018CB2803130) succeeds as GTiff.
GDALWARP: Defining SKIP_NOSOURCE=YES
Processing small_input.tif [1/1] : 0WARP: Copying metadata from first source to destination dataset
GDAL: Computing area of interest: -77.3764, 34.5167, -77.326, 34.55
OGRCT: Wrap source at -77.3512.
GTiff: ScanDirectories()
GDAL: GDALDefaultOverviews::OverviewScan()
Using internal nodata values (e.g. -9999) for image small_input.tif.
WARP: band=0 dstNoData=-9999.000000
WARP: Applying MULT_FACTOR_VERTICAL_SHIFT=1
WARP: Applying MULT_FACTOR_VERTICAL_SHIFT_PIPELINE=1
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
OGRCT: No operation matching criteria found for coordinate
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
OGRCT: No operation matching criteria found for coordinate
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
OGRCT: No operation matching criteria found for coordinate
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
OGRCT: No operation matching criteria found for coordinate
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
OGRCT: No operation matching criteria found for coordinate
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
OGRCT: No operation matching criteria found for coordinate
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
OGRCT: No operation matching criteria found for coordinate
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
OGRCT: No operation matching criteria found for coordinate
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
OGRCT: No operation matching criteria found for coordinate
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
OGRCT: No operation matching criteria found for coordinate
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
OGRCT: No operation matching criteria found for coordinate
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
OGRCT: No operation matching criteria found for coordinate
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
OGRCT: No operation matching criteria found for coordinate
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
OGRCT: No operation matching criteria found for coordinate
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
OGRCT: No operation matching criteria found for coordinate
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
OGRCT: No operation matching criteria found for coordinate
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
OGRCT: No operation matching criteria found for coordinate
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
OGRCT: No operation matching criteria found for coordinate
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
OGRCT: No operation matching criteria found for coordinate
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
OGRCT: Reprojection failed, err = 2051, further errors will be suppressed on the transform object.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
ERROR 1: PROJ: Attempt to use coordinate operation ITRF2020 + MLLW (National_Water_Level_Datum/nwldatum_4.7.0_20240621) height to ITRF2020 failed.
GDAL: GDAL_CACHEMAX = 801 MB
GDAL: GDALWarpKernel()::GWKRealCase() Src=0,0,50x33 Dst=0,0,50x33
...10...20...30...40...50...60...70...80...90...100 - done.
GDAL: Flushing dirty blocks: 0...10...20...30...40...50...60...70...80...90...100 - done.
GTiff: Adjusted bytes to write from 8000 to 5200.
GDAL: GDALClose(output_proj94.tif, this=0000018CB2803130)
GDAL: GDALClose(small_input.tif, this=0000018CB27CAF40)

and

> echo $Env:PROJ_ONLY_BEST_DEFAULT
YES

@jratike80
Copy link
Collaborator

I am looking at the black areas in the transformed results and for me they look uniformly black. Not like being the original pixels copied as-is. I suppose that those pixels have value -9999 (nodata). Isn't that the result that you wanted to achieve?

@mohammadashkezari-noaa
Copy link
Author

They look uniformly black because they are much larger than the transformed area (30m, compared to 0.5m). But they are not nodata.

@GlenRice-NOAA
Copy link
Contributor

We tried conducting a transformation with a single resolution grid using the database. This resulted in the same behavior we were seeing with the GTG where a region with a nodata value in the transformation grid passed through the original data value. This was surprising to us since we had experienced a different behavior in this ticket which was using a pipeline transformation.

gdalwarp -t_srs EPSG:9990+NOAA:98 -wo apply_vertical_shift=yes -wo sample_grid=yes -wo sample_steps=all "C:\Data\vypertemp\orig\low_res_itrf_input.tif" "C:\Data\vypertemp\orig\low_res_itrf_output_NC02_931.tif"

We then tried doing a pipeline transformation (as in the previous ticket) with GDAL Warp with both a single resolution and GTG and we were successful in honoring the nodata transformation we expected.

gdalwarp -vshift -ct "+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +inv +proj=vgridshift +grids=us_noaa_nos_MLLW-ITRF2020_2020.0_(nwldatum_4.7.0_20240621).tif +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1" "C:\Data\vypertemp\orig\low_res_itrf_input.tif" "C:\Data\vypertemp\orig\low_res_itrf_out_pipeline_nwld_931.tif"

It seems the bug we thought we had (a difference between single res and GTG behavior) was actually a difference between the behavior with a pipeline vs the database. Should we close this ticket and bring this question up in the mailing list?

@rouault rouault mentioned this issue Jan 15, 2025
Merged
@rouault
Copy link
Member

rouault commented Jan 15, 2025

I've a fix in OSGeo/PROJ#4382 that makes the following to correctly work on areas where there is nodata in the grid with

PROJ_ONLY_BEST_DEFAULT=YES PROJ_DATA=noaa_issue_11650  gdalwarp -s_srs EPSG:9989 -t_srs EPSG:9990+NOAA:98 -wo apply_vertical_shift=yes -wo sample_grid=yes -wo sample_steps=all small_input.tif  out3.tif -overwrite

why does it work correctly without PROJ_ONLY_BEST_DEFAULT=YES when you specify an explicit pipeline? The reason is that in that situation PROJ has no choice of pipelines, so it doesn't try to fallback to the dummy / ballpark transformation. The fix OSGeo/PROJ#4382 makes PROJ_ONLY_BEST_DEFAULT=YES to work properly in that context since GDAL uses proj_clone() and that function didn't propagate properly the flag set initially by PROJ_ONLY_BEST_DEFAULT=YES on the source object to the copy, thus causing the fallback ballpark transformation to still be used

@barry-gallagher
Copy link

Yay! Thank you. I was just preparing some images to go with Glen's last comment and you made a fix before I could post the images -- awesome.

@jratike80
Copy link
Collaborator

jratike80 commented Jan 15, 2025 •
edited
Loading

So, if you want to have a control on what happens, use pipeline and don't trust on defaults, especially if it is hard to find out what the defaults are. (I don't mean that I like pipelines, usually if I write one it does not work).

@GlenRice-NOAA
Copy link
Contributor

Thank you Even and Jukka! We really appreciate your help despite that we didn't post this to the right place or even ask the right question. It looks like the update is tagged for PROJ 9.6 (March). Until then we will fall back to pipelines.

@rouault
Copy link
Member

rouault commented Jan 16, 2025

We really appreciate your help despite that we didn't post this to the right place or even ask the right question

that's fine. It was far from obvious which part of the software stack was at fault. That said things are midly satisfactory with PROJ_ONLY_BEST_DEFAULT=YES causing a CPLError() to be emitted for each point failing to transform, which is both wanted, and unwanted...

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees

@rouault rouault

Labels
None yet
Projects
None yet
Milestone
No milestone
5 participants
@rouault @jratike80 @barry-gallagher @GlenRice-NOAA @mohammadashkezari-noaa