Merge pull request #10837 from OSGeo/backport-10761-to-release/3.9

[Backport release/3.9] Overview: fix nearest and mode resamplings to be exact with all data types

autotest/gcore/rasterio.py

@@ -884,29 +884,153 @@ def test_rasterio_12():
 # Test cubic resampling with masking
 
 
-def test_rasterio_13():
+@pytest.mark.parametrize(
+    "dt",
+    [
+        "Byte",
+        "Int8",
+        "Int16",
+        "UInt16",
+        "Int32",
+        "UInt32",
+        "Int64",
+        "UInt64",
+        "Float32",
+        "Float64",
+    ],
+)
+def test_rasterio_13(dt):
     numpy = pytest.importorskip("numpy")
 
-    for dt in [gdal.GDT_Byte, gdal.GDT_UInt16, gdal.GDT_UInt32]:
+    dt = gdal.GetDataTypeByName(dt)
+    mem_ds = gdal.GetDriverByName("MEM").Create("", 4, 3, 1, dt)
+    mem_ds.GetRasterBand(1).SetNoDataValue(0)
+    if dt == gdal.GDT_Int8:
+        x = (1 << 7) - 1
+    elif dt == gdal.GDT_Byte:
+        x = (1 << 8) - 1
+    elif dt == gdal.GDT_Int16:
+        x = (1 << 15) - 1
+    elif dt == gdal.GDT_UInt16:
+        x = (1 << 16) - 1
+    elif dt == gdal.GDT_Int32:
+        x = (1 << 31) - 1
+    elif dt == gdal.GDT_UInt32:
+        x = (1 << 32) - 1
+    elif dt == gdal.GDT_Int64:
+        x = (1 << 63) - 1
+    elif dt == gdal.GDT_UInt64:
+        x = (1 << 64) - 2048
+    elif dt == gdal.GDT_Float32:
+        x = 1.5
+    else:
+        x = 1.23456
+    mem_ds.GetRasterBand(1).WriteArray(
+        numpy.array([[0, 0, 0, 0], [0, x, 0, 0], [0, 0, 0, 0]])
+    )
 
-        mem_ds = gdal.GetDriverByName("MEM").Create("", 4, 3, 1, dt)
-        mem_ds.GetRasterBand(1).SetNoDataValue(0)
-        mem_ds.GetRasterBand(1).WriteArray(
-            numpy.array([[0, 0, 0, 0], [0, 255, 0, 0], [0, 0, 0, 0]])
-        )
+    ar_ds = mem_ds.ReadAsArray(
+        0, 0, 4, 3, buf_xsize=8, buf_ysize=3, resample_alg=gdal.GRIORA_Cubic
+    )
 
-        ar_ds = mem_ds.ReadAsArray(
-            0, 0, 4, 3, buf_xsize=8, buf_ysize=3, resample_alg=gdal.GRIORA_Cubic
-        )
+    expected_ar = numpy.array(
+        [
+            [0, 0, 0, 0, 0, 0, 0, 0],
+            [0, x, x, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0],
+        ]
+    )
+    assert numpy.array_equal(ar_ds, expected_ar)
 
-        expected_ar = numpy.array(
-            [
-                [0, 0, 0, 0, 0, 0, 0, 0],
-                [0, 255, 255, 0, 0, 0, 0, 0],
-                [0, 0, 0, 0, 0, 0, 0, 0],
-            ]
-        )
-        assert numpy.array_equal(ar_ds, expected_ar), (ar_ds, dt)
+
+###############################################################################
+# Test nearest and mode resampling
+
+
+@pytest.mark.parametrize(
+    "dt",
+    [
+        "Byte",
+        "Int8",
+        "Int16",
+        "UInt16",
+        "Int32",
+        "UInt32",
+        "Int64",
+        "UInt64",
+        "Float32",
+        "Float64",
+        "CInt16",
+        "CInt32",
+        "CFloat32",
+        "CFloat64",
+    ],
+)
+@pytest.mark.parametrize(
+    "resample_alg", [gdal.GRIORA_NearestNeighbour, gdal.GRIORA_Mode]
+)
+@pytest.mark.parametrize("use_nan", [True, False])
+def test_rasterio_nearest_or_mode(dt, resample_alg, use_nan):
+    numpy = pytest.importorskip("numpy")
+    gdal_array = pytest.importorskip("osgeo.gdal_array")
+
+    dt = gdal.GetDataTypeByName(dt)
+    mem_ds = gdal.GetDriverByName("MEM").Create("", 4, 4, 1, dt)
+    if dt == gdal.GDT_Int8:
+        x = (1 << 7) - 1
+    elif dt == gdal.GDT_Byte:
+        x = (1 << 8) - 1
+    elif dt == gdal.GDT_Int16 or dt == gdal.GDT_CInt16:
+        x = (1 << 15) - 1
+    elif dt == gdal.GDT_UInt16:
+        x = (1 << 16) - 1
+    elif dt == gdal.GDT_Int32 or dt == gdal.GDT_CInt32:
+        x = (1 << 31) - 1
+    elif dt == gdal.GDT_UInt32:
+        x = (1 << 32) - 1
+    elif dt == gdal.GDT_Int64:
+        x = (1 << 63) - 1
+    elif dt == gdal.GDT_UInt64:
+        x = (1 << 64) - 1
+    elif dt == gdal.GDT_Float32 or dt == gdal.GDT_CFloat32:
+        x = float("nan") if use_nan else 1.5
+    else:
+        x = float("nan") if use_nan else 1.234567890123
+
+    if gdal.DataTypeIsComplex(dt):
+        val = complex(x, x)
+    else:
+        val = x
+
+    dtype = gdal_array.flip_code(dt)
+    mem_ds.GetRasterBand(1).WriteArray(numpy.full((4, 4), val, dtype=dtype))
+
+    ar_ds = mem_ds.ReadAsArray(
+        0, 0, 4, 4, buf_xsize=1, buf_ysize=1, resample_alg=resample_alg
+    )
+
+    expected_ar = numpy.array([[val]]).astype(dtype)
+    if math.isnan(x):
+        if gdal.DataTypeIsComplex(dt):
+            assert math.isnan(ar_ds[0][0].real) and math.isnan(ar_ds[0][0].imag)
+        else:
+            assert math.isnan(ar_ds[0][0])
+    else:
+        assert numpy.array_equal(ar_ds, expected_ar)
+
+    resample_alg_mapping = {
+        gdal.GRIORA_NearestNeighbour: "NEAR",
+        gdal.GRIORA_Mode: "MODE",
+    }
+    mem_ds.BuildOverviews(resample_alg_mapping[resample_alg], [4])
+    ar_ds = mem_ds.GetRasterBand(1).GetOverview(0).ReadAsArray()
+    if math.isnan(x):
+        if gdal.DataTypeIsComplex(dt):
+            assert math.isnan(ar_ds[0][0].real) and math.isnan(ar_ds[0][0].imag)
+        else:
+            assert math.isnan(ar_ds[0][0])
+    else:
+        assert numpy.array_equal(ar_ds, expected_ar)
 
 
 ###############################################################################

gcore/gdalnodatamaskband.cpp

@@ -140,7 +140,8 @@ static GDALDataType GetWorkDataType(GDALDataType eDataType)
             eWrkDT = eDataType;
             break;
 
-        default:
+        case GDT_Unknown:
+        case GDT_TypeCount:
             CPLAssert(false);
             eWrkDT = GDT_Float64;
             break;