OPTIMIZE_SIZE Warper option from gdalwarp gives different re-projection output

[aafaque33]

What is the bug?

Recently we tried to upgrade the GDAL version in the docker container.

After upgrade of the GDAL from 3.5.3 to GDAL 3.9.1 we have found that one of the re-projection output was having differences (not in all pixels but kind of scan lines with slight pixel shift)
Verified that same shift happens in GDAL version 3.8

First we thought that it may be because of the rasterio as the issue was while using the rio warp but not with gdalwarp so logged the issue / discussion in raserio

• Look at detailed discussion here: `rio warp` different reprojection result for different `gdal` versions (Slight Pixel Shifts) rasterio/rasterio#3150

After discussion it was found that same issue haapens in gdalwarp if the -co TILED=YES -co BLOCKXSIZE=512 -co BLOCKYSIZE=512 which rio warp is using by default.

It turns out after so much trial and error that it was because of the one of the warper option OPTIMIZE_SIZE that seems like before GDAL 3.8 was always FALSE for everything, But since GDAL 3.8 is automatically enabled (may be for TILED / BLOCKSIZE / Compressed options) when it is safe to do so, and that produces a slight pixel shift / difference in reprojected output

Testing the gdalwarp post GDAL 3.8 with that warper option disabled gives the same reprojected output as it was in older versions of GDAL

Not sure if it is expected behavior but shouldn't the reprojected output be same with or without that warper option?

Steps to reproduce the issue

Note: Reprojection is from UTM Zone 22N to Geographic (EPSG:4326)

Data Files

Tif file with no reprojection done (Currently in UTM Zone 22N). Providing the google drive link as can't upload more than 25 MB file

• https://drive.google.com/file/d/1X3eP7O4nihHM7_mONIlekP6o5Dv7S1ut/view?usp=drive_link

Run a docker container

> docker run --rm -it -v ~/docker:/data ghcr.io/osgeo/gdal:ubuntu-full-3.9.2 /bin/bash

Run GDAL WARP with different scenarios**

> root@120bf665500e:/# gdalinfo --version
GDAL 3.9.2, released 2024/08/13

GDAL WARP - Non Tiled creation option with default Optimize Size

> root@120bf665500e:/# gdalwarp -t_srs 'EPSG:4326' -r near -tr 0.00026922090197700754 0.00026922090197700754 -tap /data/rasters/HLS_UTM_Zone_22N_Mosaiced.tif /data/rasters/HLS_UTM_Zone_22N_GDALWARP_GDAL3.9.2_Non_TILED_Reprojection.tif

Creating output file that is 16473P x 7055L.

GDAL WARP - Non Tiled creation option with Optimize Size disabled

> root@120bf665500e:/# gdalwarp -t_srs 'EPSG:4326' -r near -tr 0.00026922090197700754 0.00026922090197700754 -tap -wo OPTIMIZE_SIZE=FALSE /data/rasters/HLS_UTM_Zone_22N_Mosaiced.tif /data/rasters/HLS_UTM_Zone_22N_GDALWARP_GDAL3.9.2_Non_TILED_Reprojection_No_Optimize_Size.tif

Creating output file that is 16473P x 7055L.

GDAL WARP - Non Tiled creation option with Optimize Size enabled

> root@120bf665500e:/# gdalwarp -t_srs 'EPSG:4326' -r near -tr 0.00026922090197700754 0.00026922090197700754 -tap -wo OPTIMIZE_SIZE=TRUE /data/rasters/HLS_UTM_Zone_22N_Mosaiced.tif /data/rasters/HLS_UTM_Zone_22N_GDALWARP_GDAL3.9.2_Non_TILED_Reprojection_Optimize_Size.tif

Creating output file that is 16473P x 7055L.

GDAL WARP - Tiled creation option with default Optimize Size

> root@120bf665500e:/# gdalwarp -t_srs 'EPSG:4326' -r near -tr 0.00026922090197700754 0.00026922090197700754 -tap -co TILED=YES -co BLOCKXSIZE=512 -co BLOCKYSIZE=512  /data/raste
rs/HLS_UTM_Zone_22N_Mosaiced.tif /data/rasters/HLS_UTM_Zone_22N_GDALWARP_GDAL3.9.2_TILED_Reprojection.tif

Creating output file that is 16473P x 7055L.

GDAL WARP - Tiled creation option with Optimize Size disabled

> root@120bf665500e:/# gdalwarp -t_srs 'EPSG:4326' -r near -tr 0.00026922090197700754 0.00026922090197700754 -tap -wo OPTIMIZE_SIZE=FALSE -co TILED=YES -co BLOCKXSIZE=512 -co BLOCKYSIZE=512 /data/rasters/HLS_UTM_Zone_22N_Mosaiced.tif /data/rasters/HLS_UTM_Zone_22N_GDALWARP_GDAL3.9.2_TILED_Reprojection_No_Optimize_Size.tif

Creating output file that is 16473P x 7055L.

GDAL WARP - Tiled creation option with Optimize Size enabled

> root@120bf665500e:/# gdalwarp -t_srs 'EPSG:4326' -r near -tr 0.00026922090197700754 0.00026922090197700754 -tap -wo OPTIMIZE_SIZE=TRUE -co TILED=YES -co BLOCKXSIZE=512 -co BLOCKYSIZE=512 /data/rasters/HLS_UTM_Zone_22N_Mosaiced.tif /data/rasters/HLS_UTM_Zone_22N_GDALWARP_GDAL3.9.2_TILED_Reprojection_Optimize_Size.tif

Creating output file that is 16473P x 7055L.

Differences

The output that ran with -co TILED=YES -co BLOCKXSIZE=512 -co BLOCKYSIZE=512 with Default OPTIMIZE_SIZE that auto enables or with OPTIMIZE_SIZE as TRUE differs from the gdalwarp that didn't use the -co TILED=YES -co BLOCKXSIZE=512 -co BLOCKYSIZE=512 or used the OPTIMIZE_SIZE as FALSE

Opened the outputs from above re-projections in QGIS and ran the raster diff

• The outputs displayed on QGIS
  • 
• The diff of the Non Tiled vs the Tiled (that uses the auto OPTIMIZE_SIZE). From Zoom out it seems no differences
  • 
• Zoomed in Version of the diff you can see scan lines of difference
  • 
• Further Zoom, you can see how the pixel of Non Tiled / TILED with OPTIMIZE Size as FALSE differ from the TILED with Auto OPTIMIZE Size and Optimize size as True
  • 
  • Look at the video below for a demonstration of pixel change / shift (Old video used for rasterio issue but relatable for gdalwarp tiled vs non tiled)
  • https://github.com/user-attachments/assets/71225242-18f5-4fae-9e9a-1046377440a2

Versions and provenance

• GDAL versions
  • 3.9.2
  • 3.9.1
• Python version (3.11.9 (main, Jul 25 2024, 20:22:46) [GCC 8.5.0 20210514 (Red Hat 8.5.0-22)])
• Operation System Information (Rocky Linux 8.10 (Green Obsidian))

Additional context

Same details discussed on the rasterio repo which may have some extra details

• `rio warp` different reprojection result for different `gdal` versions (Slight Pixel Shifts) rasterio/rasterio#3150

[rouault]

Not sure if it is expected behavior

Yes

but shouldn't the reprojected output be same with or without that warper option?

In theory yes, but in practice, this is hard to achieve. If your reprojection is relative uniform over the whole dataset (i.e. the ratio of source pixels needed to compute a target pixel) is approximately constant over it, you could specify "-wo XSCALE= -wo YSCALE=" to the gdalwarp command line. If you don't change resolution, the value would be 1. If you downsample by a factor of a half, that would be 0.5, etc.
But in general that ratio is not constant, so we update it it for each "warping chunk". So the size of the warping chunk has an influence on the output.

[aafaque33]

I know there is a Change in 3.8 that makes OPTIMIZE_SIZE to be auto enabled when needed so, just curious if there a specific reason for that ?

Wondering if that can be default to FALSE again somehow? The reason being the other libs using gdal , i.e rasterio or any other python packages gives different results with different versions of GDAL only (same library / package version itself) and that was so confusing and time taking to find out why.

Or if there's way to define in the documentation that this option that is now defaults to True (when ever needed) can yield to different results so make sure to disable it if the olde results are needed?
https://gdal.org/en/latest/api/gdalwarp_cpp.html#_CPPv4N15GDALWarpOptions16papszWarpOptionsE:~:text=within%20the%20polygon.-,OPTIMIZE_SIZE%3A,-This%20defaults%20to

[rouault]

just curious if there a specific reason for that ?

yes, cf commit 9bd899a: "Warper: auto-enable OPTIMIZE_SIZE warping option when reasonable (fixes #7761)"

Wondering if that can be default to FALSE again somehow?

Not unless there would a major reason to do so, like bad results, very poor performance. Not just little pixel difference.

The reason being the other libs using gdal , i.e rasterio or any other python packages gives different results with different versions of GDAL only (same library / package version itself) and that was so confusing and time taking to find out why.

gdalwarp results shouldn't be expected to remain identical among GDAL (and PROJ) versions, due to various adjustments in heuristics, precision of computations, etc.

Or if there's way to define in the documentation that this option that is now defaults to True (when ever needed) can yield to different results so make sure to disable it if the olde results are needed?

Here you are: #10814

[aafaque33]

Thanks @rouault . I think I got my answer and we are able to get the correct or can say previous expected results via the OPTIMIZE_SIZE=FALSE after gdal upgrade which was great.

Thanks for all the explanations !!