Gdalwarp assigning zero values outside input raster extents

Question

I'm attempting to combine several GeoTIFF files into a single mosaic raster using gdalwarp:

import os
os.environ['AWS_NO_SIGN_REQUEST'] = 'YES'

file1 = "https://dea-public-data-dev.s3-ap-southeast-2.amazonaws.com/derivative/ga_s2ls_intertidal_cyear_3/1-0-0/x085/y140/2020--P1Y/ga_s2ls_intertidal_cyear_3_x085y140_2020--P1Y_final_elevation.tif"
file2 = "https://dea-public-data-dev.s3-ap-southeast-2.amazonaws.com/derivative/ga_s2ls_intertidal_cyear_3/1-0-0/x092/y144/2020--P1Y/ga_s2ls_intertidal_cyear_3_x092y144_2020--P1Y_final_elevation.tif"
!gdalwarp $file1 $file2 test.tif -b 1 -overwrite -multi -wm 80% -co NUM_THREADS=ALL_CPUS -of COG -co COMPRESS=ZSTD -co PREDICTOR=YES -co OVERVIEWS=FORCE_USE_EXISTING

While my output mosaic raster correctly includes my input data (red circles), areas outside my rasters get assigned zero values (grey pixels):

In the past, I've had similar issues when nodata values are set incorrectly. However, as far as I can tell from gdalinfo, my input data is correctly configured with a float32 datatype and NaN nodata values. I've tried many combinations of -srcnodata and -dstnodata to try and get these artificial zero values set to nodata (e.g. -dstnodata nan, -srcnodata nan), but nothing I have tried so far has had any effect. Another challenge is that my actual data spans from -3 to +3, and includes "correct" zero values I don't want to accidently lose.

This issue seems specific to rasters with a floating point data type - interestingly, I don't have the same issue when attempting to merge integer datatype rasters - in that case, areas outside my rasters are correctly set to nodata:

file1 = "https://dea-public-data-dev.s3-ap-southeast-2.amazonaws.com/derivative/ga_s2ls_intertidal_cyear_3/1-0-0/x085/y140/2020--P1Y/ga_s2ls_intertidal_cyear_3_x085y140_2020--P1Y_final_exposure.tif"
file2 = "https://dea-public-data-dev.s3-ap-southeast-2.amazonaws.com/derivative/ga_s2ls_intertidal_cyear_3/1-0-0/x092/y144/2020--P1Y/ga_s2ls_intertidal_cyear_3_x092y144_2020--P1Y_final_exposure.tif"
!gdalwarp $file1 $file2 test.tif -b 1 -overwrite -multi -wm 80% -co NUM_THREADS=ALL_CPUS -of COG -co COMPRESS=ZSTD -co PREDICTOR=YES -co OVERVIEWS=FORCE_USE_EXISTING

Answer 1

I do not know what happens, but warping directly into COG does not honour the nodata setting. BTW. are you sure about the nodata in the source images? My gdalinfo does not show nodata. Therefore I added -dstnodata into the gdalwarp command.

gdalwarp ga_s2ls_intertidal_cyear_3_x085y140_2020--P1Y_final_elevation.tif  ga_s2ls_intertidal_cyear_3_x092y144_2020--P1Y_final_elevation.tif test3.tif -b 1 -overwrite -multi -wm 80% -co NUM_THREADS=ALL_CPUS -of COG -co Compress=ZSTD  -dstnodata nan

Despite having -dstnodata in the command, the resulting COG does not have nodata. But if I write GeoTIFF then it does have nodata. Also when I converted this GeoTIFF into COG the nodata remains. Alternatively, defining nodata works if the mosaic is make first with gdalbuildvrt.

gdalbuildvrt test.vrt ga_s2ls_intertidal_cyear_3_x085y140_2020--P1Y_final_elevation.tif ga_s2ls_intertidal_cyear_3_x092y144_2020--P1Y_final_elevation.tif

gdalwarp test.vrt testvrt.tif -b 1 -overwrite -multi -wm 80% -co NUM_THREADS=ALL_CPUS -of COG -co COMPRESS=ZSTD -co PREDICTOR=YES -co OVERVIEWS=FORCE_USE_EXISTING -dstnodata nan

You may have found a bug so feel free to write mail into gdal-dev mailing list. However, I would not use your workflow. Because there is no I would create the mosaic with gdalbuildvrt and because the CRS remains the same I would use gdal_translate instead of gdalwarp.

gdalbuildvrt test.vrt ga_s2ls_intertidal_cyear_3_x085y140_2020--P1Y_final_elevation.tif ga_s2ls_intertidal_cyear_3_x092y144_2020--P1Y_final_elevation.tif -vrtnodata nan

gdal_translate -of COG test.vrt test5.tif