OSGeo / gdal

GDAL is an open source MIT licensed translator library for raster and vector geospatial data formats.
https://gdal.org
Other
4.91k stars 2.55k forks source link

gdalwarp overview selection issue due to numerical instability #10873

Closed calogeromauceri closed 1 month ago

calogeromauceri commented 1 month ago

What is the bug?

Summary When reprojecting a file using gdalwarp, the overview selected can unexpectedly change by slightly shifting the output window. This results in visible differences in the output image, as lower resolution overviews might be chosen.

Example Consider the following example, where the same input file produces different output with a small shift in the target extent

gdalwarp -of PNG -ts 256 256 -t_srs "+proj=eqc +a=2439700 +b=2439700 +lon_0=0.000000 +lat_0=0.000000 +no_defs +x_0=0 +y_0=0" -te -3481372.584110 1173100.739750 -2713372.607150 1941100.716710 -overwrite test.tif out1.png

This selects the second overview. out1

shifting a little bit the extent gdalwarp -of PNG -ts 256 256 -t_srs "+proj=eqc +a=2439700 +b=2439700 +lon_0=0.000000 +lat_0=0.000000 +no_defs +x_0=0 +y_0=0" -te -3484372.584010 1173100.739750 -2716372.607050 1941100.716710 -overwrite test.tif out2.png

This selects the first, higher resolution overview out2

Source problem and suggested fix Upon debugging, I found the problem is due to numeric instability when checking ratio between the desired resolution and the available overviews. Specifically in the GDALWarpDirect there is a comparison for selecting the right overview.

Relevant lines of code: GDALWarpDirect

if (dfOvrRatio < dfTargetRatio &&
    dfNextOvrRatio > dfTargetRatio)
    break;

Using an epsilon to stabilize the ratio checks seems to fix the issue. For example

constexpr double EPS = 1e-1;
if (dfOvrRatio - dfTargetRatio < EPS &&
   dfNextOvrRatio - dfTargetRatio > EPS )
   break;

You can reproduce the issue using the following test file: https://files.actgate.com/temp/test.tif

Steps to reproduce the issue

gdalwarp -of PNG -ts 256 256 -t_srs "+proj=eqc +a=2439700 +b=2439700 +lon_0=0.000000 +lat_0=0.000000 +no_defs +x_0=0 +y_0=0" -te -3481372.584110 1173100.739750 -2713372.607150 1941100.716710 -overwrite test.tif out1.png

gdalwarp -of PNG -ts 256 256 -t_srs "+proj=eqc +a=2439700 +b=2439700 +lon_0=0.000000 +lat_0=0.000000 +no_defs +x_0=0 +y_0=0" -te -3484372.584010 1173100.739750 -2716372.607050 1941100.716710 -overwrite test.tif out2.png

Versions and provenance

3.9.2

Additional context

No response

jratike80 commented 1 month ago

I cannot judge if your suggestion is good, but if you want to be sure about what overview level gdalwarp takes it is best to define it with -ovrhttps://gdal.org/en/latest/programs/gdalwarp.html#cmdoption-gdalwarp-ovr.

calogeromauceri commented 1 month ago

I confirm I get a consistent result if I pass the overview to use with -ovr, but I would like to use the -ovr AUTO option and let gdalwarp select the right overview to use, otherwise I would need to implement/duplicate what is already done by gdalwarp for the automatic overview selection