Bad gdalwarp interpolation when using resampling other than nearest

[calogeromauceri]

What is the bug?

gdalwarp is producing a bad output when using interpolation with the following command

gdalwarp -of PNG -r bilinear -t_srs '+proj=eqc +a=2439700 +b=2439700 +lon_0=0.000000 +lat_0=0.000000 +x_0=0 +y_0=0 +units=m +no_defs' -ts 256 256 -te -958067.949620 -958067.950079 0.000460 0.000000 -overwrite msgr_dem_low_res.tif tile.png

The output image is kind of blurred/low resolution

On the other side, it looks good if nearest interpolation is used

gdalwarp -of PNG -t_srs '+proj=eqc +a=2439700 +b=2439700 +lon_0=0.000000 +lat_0=0.000000 +x_0=0 +y_0=0 +units=m +no_defs' -ts 256 256 -te -958067.949620 -958067.950079 0.000460 0.000000 -overwrite msgr_dem_low_res.tif tile.png tile.png

The problem seems related to the fact that the right edge of the tile is around longitude 0, while the source image is a full globe/360 degrees image centered at 180 longitude

Avoiding the 0 longitude at the edge, for example reducing the tile size to 255x255 (and adjusting the window coordinates accordingly), the interpolated output file looks much better

gdalwarp -of PNG -r bilinear -t_srs '+proj=eqc +a=2439700 +b=2439700 +lon_0=0.000000 +lat_0=0.000000 +x_0=0 +y_0=0 +units=m +no_defs' -ts 255 255 -te -958067.949620 -950583.044220 -7484.905400 0.000000 -overwrite msgr_dem_low_res.tif tile.png

I'm using the mercury dem from here
https://pdsimage2.wr.usgs.gov/Messenger/MESSDEM_1001/DEM/GLOBAL/IMG/
Here is a lower resolution of the same dem for testing purpose
https://files.actgate.com/terrain/mercury/msgr_dem_low_res.tif

Steps to reproduce the issue

gdalwarp -of PNG -r bilinear -t_srs '+proj=eqc +a=2439700 +b=2439700 +lon_0=0.000000 +lat_0=0.000000 +x_0=0 +y_0=0 +units=m +no_defs' -ts 256 256 -te -958067.949620 -958067.950079 0.000460 0.000000 -overwrite msgr_dem_low_res.tif tile.png

Versions and provenance

GDAL 3.9.1

Additional context

No response

[calogeromauceri]

Using undocumented XSCALE and YSCALE options seems to fix the issue, but I'm not sure what those options mean and if they can cause other side effects.
We are exploring using gdalwarp C++ API in a GIS application that renders any GDAL supported file. It is important it works at any resolution on any valid window.

gdalwarp -of PNG -r bilinear -t_srs '+proj=eqc +a=2439700 +b=2439700 +lon_0=0.000000 +lat_0=0.000000 +x_0=0 +y_0=0 +units=m +no_defs' -ts 256 256 -te -958067.949620 -958067.950079 0.000460 0.000000 -overwrite -wo XSCALE=1 -wo YSCALE=1 msgr_dem_low_res.tif tile.png

[jjimenezshaw]

Using undocumented XSCALE and YSCALE

Looking at the comment in the code it is probably better to let the algorithm to compute the sampling with just
-wo XSCALE=FROM_GRID_SAMPLING

    // If the xscale is significantly lower than the yscale, this is highly
    // suspicious of a situation of wrapping a very large virtual file in
    // geographic coordinates with left and right parts being close to the
    // antimeridian. In that situation, the xscale computed by the above method
    // is completely wrong. Prefer doing an average of a few sample points
    // instead

[calogeromauceri]

Very interesting. Thanks for pointing out that solution. It seems to work. I guess that's another hidden setting as I was not able to find any documentation to it.

Digging into the code after that comment though, I found out that when the rendering of this image works the
ratio
dfYScale / dfXScale is near ~1.0. When the tile is misrendered, then the ratio is much higher, ~16. The threshold for entering that block is set to 100, that's why that logic is not executed by default.

[jjimenezshaw]

It is the never ending problem with the antimeridian.
In the Earth most of the cases do not cross it.

You do not have an empty Pacific Ocean in Mercury ;)

As you mention, maybe the threshold to enter into that code should be much lower... Or even do it always? I do not know the implications.

[calogeromauceri]

Now I wonder I should always use that option even when the nearest resampling is used. In that case, the rendering looks better but the ratio dfYScale / dfXScale is still wrong in the case I mentioned above.