Merge pull request OSGeo#13103 from rouault/fix_13100

gdaldem: fix wrong results on non north-up src ds wih aspect/tpi/tri (and…

Author: rouault

apps/gdaldem_lib.cpp

@@ -3634,6 +3634,39 @@ GDALDatasetH GDALDEMProcessing(const char *pszDest, GDALDatasetH hSrcDataset,
 
     GDALDEMProcessingOptions *psOptions = psOptionsToFree.get();
 
+    double adfGeoTransform[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
+
+    // Keep that variable in that scope.
+    // A reference on that dataset will be taken by GDALGeneric3x3Dataset
+    // (and released at its destruction) if we go to that code path, and
+    // GDALWarp() (actually the VRTWarpedDataset) will itself take a reference
+    // on hSrcDataset.
+    std::unique_ptr<GDALDataset, GDALDatasetUniquePtrReleaser> poTmpSrcDS;
+
+    if (GDALGetGeoTransform(hSrcDataset, adfGeoTransform) == CE_None &&
+        // For following modes, detect non north-up datasets and autowarp
+        // them so they are north-up oriented.
+        (((eUtilityMode == ASPECT || eUtilityMode == TRI ||
+           eUtilityMode == TPI) &&
+          (adfGeoTransform[2] != 0 || adfGeoTransform[4] != 0 ||
+           adfGeoTransform[5] > 0)) ||
+         // For following modes, detect rotated datasets and autowarp
+         // them so they are north-up oriented. (south-up are fine for those)
+         ((eUtilityMode == SLOPE || eUtilityMode == HILL_SHADE) &&
+          (adfGeoTransform[2] != 0 || adfGeoTransform[4] != 0))))
+    {
+        const char *const apszWarpOptions[] = {"-of", "VRT", nullptr};
+        GDALWarpAppOptions *psWarpOptions = GDALWarpAppOptionsNew(
+            const_cast<char **>(apszWarpOptions), nullptr);
+        poTmpSrcDS.reset(GDALDataset::FromHandle(
+            GDALWarp("", nullptr, 1, &hSrcDataset, psWarpOptions, nullptr)));
+        GDALWarpAppOptionsFree(psWarpOptions);
+        if (!poTmpSrcDS)
+            return nullptr;
+        hSrcDataset = GDALDataset::ToHandle(poTmpSrcDS.get());
+        GDALGetGeoTransform(hSrcDataset, adfGeoTransform);
+    }
+
     const int nXSize = GDALGetRasterXSize(hSrcDataset);
     const int nYSize = GDALGetRasterYSize(hSrcDataset);
 
@@ -3730,12 +3763,8 @@ GDALDatasetH GDALDEMProcessing(const char *pszDest, GDALDatasetH hSrcDataset,
         return nullptr;
     }
 
-    double adfGeoTransform[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
-
     GDALRasterBandH hSrcBand = GDALGetRasterBand(hSrcDataset, psOptions->nBand);
 
-    GDALGetGeoTransform(hSrcDataset, adfGeoTransform);
-
     CPLString osFormat;
     if (psOptions->osFormat.empty())
     {

autotest/utilities/test_gdaldem_lib.py

@@ -910,3 +910,38 @@ def test_gdaldem_lib_dict_arguments():
     ind = opt.index("-co")
 
     assert opt[ind : ind + 4] == ["-co", "COMPRESS=DEFLATE", "-co", "LEVEL=4"]
+
+
+###############################################################################
+
+
+def flipped_ds():
+
+    src_ds = gdal.Open("../gdrivers/data/n43.tif")
+    gt = src_ds.GetGeoTransform()
+    ulx = "%.17g" % gt[0]
+    uly = "%.17g" % gt[3]
+    lrx = "%.17g" % (gt[0] + src_ds.RasterXSize * gt[1])
+    lry = "%.17g" % (gt[3] + src_ds.RasterYSize * gt[5])
+    return gdal.Warp("", src_ds, options=f"-of VRT -te {ulx} {uly} {lrx} {lry}")
+
+
+@pytest.mark.parametrize("alg", ["aspect", "TPI", "TRI"])
+def test_gdaldem_lib_flipped_aspect_tpi_tri(alg):
+
+    src_ds = gdal.Open("../gdrivers/data/n43.tif")
+    ref_ds = gdal.DEMProcessing("", src_ds, alg, format="MEM", computeEdges=True)
+    out_ds = gdal.DEMProcessing("", flipped_ds(), alg, format="MEM", computeEdges=True)
+    assert ref_ds.GetGeoTransform() == pytest.approx(out_ds.GetGeoTransform())
+    assert ref_ds.ReadRaster() == out_ds.ReadRaster()
+
+
+@pytest.mark.parametrize("alg", ["hillshade", "slope", "roughness"])
+def test_gdaldem_lib_flipped_hillshade_slope_roughness(alg):
+
+    src_ds = gdal.Open("../gdrivers/data/n43.tif")
+    ref_ds = gdal.DEMProcessing("", src_ds, alg, format="MEM", computeEdges=True)
+    out_ds = gdal.DEMProcessing("", flipped_ds(), alg, format="MEM", computeEdges=True)
+    out_ds = gdal.Warp("", out_ds, format="MEM")
+    assert ref_ds.GetGeoTransform() == pytest.approx(out_ds.GetGeoTransform())
+    assert ref_ds.ReadRaster() == out_ds.ReadRaster()