gdal warp changes 0's into 1's when called with a specific combination of options

[vincentschut]

What is the bug?

Not too long ago we upgraded our systems from gdal 3.9.1 to 3.10.1. After the upgrade, we found that all our binary masks were filled with 1's.

This appeared to be caused by a very specific combination of warp options. A minimal reproducer is given below. It starts with an uint8 array full of zeros, and after warping, all the data has become 1's. Note that all the given warp options seem necessary to reproduce the issue, where the resampling can be anything except nearest, the dstNodata must be 0, and the outputType can be any unsigned integer type, to reproduce the issue.

Steps to reproduce the issue

import numpy as np
from osgeo import gdal, gdal_array

gdal.UseExceptions()

data = np.zeros((3, 4), np.uint8)
print('before:', data)
ds = gdal_array.SaveArray(data, '/vsimem/data.tif')
ds.SetGeoTransform((0, 1, 0, 0, 0, 1))
warpoptions = {
    'format': 'GTiff',
    'resampleAlg': gdal.GRA_Bilinear,  # any non-nearest resampling
    "dstNodata": 0,
    "outputType": gdal.GDT_Byte,  # any unsigned int type
}
ds2 = gdal.Warp('/vsimem/data2.tif', ds, **warpoptions)
data2 = ds2.GetRasterBand(1).ReadAsArray()
print('after:', data2)

For me, this outputs:

before: [[0 0 0 0]
 [0 0 0 0]
 [0 0 0 0]]
after: [[1 1 1 1]
 [1 1 1 1]
 [1 1 1 1]]

Versions and provenance

OS: ubuntu (gdal's own ubuntu-small docker)
version where this worked: 3.9.1
version with this bug: 3.10.1

Additional context

No response

[rouault]

This is not a bug, but a feature (I would tend to think this behavior would be much older than a 3.10 change). Your source data includes 0 as a valid (non nodata value) value, so once warped, it would be incorrect to map it to the dest noData value, hence it is replaced from 0 to the closed non-zero value, 1.

Either remove dstNodata = 0 or add srcNodata = 0.

The fact that 0 is not replaced by 1 in the nearest neighbour case is actually a bug

[vincentschut]

Ah, right, seems legit. Though a warning would have been appreciated ;-) (just like we get a warning when specifying nan as nodata value for a integer dataset).
This being a change between 3.9 and 3.10 might also be due to that our code was a bit more complex than the mmre above, and the result of multiple refactorings over time.
Anyway, the current situation is at least not consistent. Not only regarding the resampling, but als regarding the output type. If I change that to GDT_Float32 in the mmre above, I get zero's as output, even with the dstNodata=0.