From 3a911bfa41fb4b9f463df56aedaf097cea20eb0a Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Fri, 11 Oct 2024 01:18:12 +0200
Subject: [PATCH] gdalwarp: make sure to check that angular unit is degree for
 heuristics related to geographic CRS

Fixes #10975
---
 alg/gdaltransformer.cpp                     |  25 +++++++++++++-------
 autotest/utilities/data/geog_arc_second.tif | Bin 0 -> 587 bytes
 autotest/utilities/test_gdalwarp_lib.py     |  16 +++++++++++++
 3 files changed, 32 insertions(+), 9 deletions(-)
 create mode 100644 autotest/utilities/data/geog_arc_second.tif

diff --git a/alg/gdaltransformer.cpp b/alg/gdaltransformer.cpp
index 3cd1c2edf6d9..ba2adc3dc73e 100644
--- a/alg/gdaltransformer.cpp
+++ b/alg/gdaltransformer.cpp
@@ -869,7 +869,7 @@ CPLErr CPL_STDCALL GDALSuggestedWarpOutput2(GDALDatasetH hSrcDS,
                  "transform.",
                  nFailedCount, nSamplePoints);
 
-    bool bIsGeographicCoords = false;
+    bool bIsGeographicCoordsDeg = false;
     if (bIsGDALGenImgProjTransform)
     {
         const GDALGenImgProjTransformInfo *pGIPTI =
@@ -899,15 +899,17 @@ CPLErr CPL_STDCALL GDALSuggestedWarpOutput2(GDALDatasetH hSrcDS,
                     CSLFetchNameValue(pGLTI->papszGeolocationInfo, "SRS");
                 if (pszGLSRS == nullptr)
                 {
-                    bIsGeographicCoords = true;
+                    bIsGeographicCoordsDeg = true;
                 }
                 else
                 {
                     OGRSpatialReference oSRS;
                     if (oSRS.SetFromUserInput(pszGLSRS) == OGRERR_NONE &&
-                        oSRS.IsGeographic())
+                        oSRS.IsGeographic() &&
+                        std::fabs(oSRS.GetAngularUnits() -
+                                  CPLAtof(SRS_UA_DEGREE_CONV)) < 1e-9)
                     {
-                        bIsGeographicCoords = true;
+                        bIsGeographicCoordsDeg = true;
                     }
                 }
             }
@@ -967,7 +969,7 @@ CPLErr CPL_STDCALL GDALSuggestedWarpOutput2(GDALDatasetH hSrcDS,
                      CPLAtof(SRS_UA_DEGREE_CONV)) < 1e-9 &&
                 (!poSourceCRS || !poSourceCRS->IsGeographic()))
             {
-                bIsGeographicCoords = true;
+                bIsGeographicCoordsDeg = true;
 
                 std::unique_ptr<CPLConfigOptionSetter> poSetter;
                 if (pGIPTI->bCheckWithInvertPROJ)
@@ -1194,7 +1196,7 @@ CPLErr CPL_STDCALL GDALSuggestedWarpOutput2(GDALDatasetH hSrcDS,
     /*      Recompute some bounds so that all return values are consistent  */
     /* -------------------------------------------------------------------- */
     double dfMaxXOutNew = dfMinXOut + (*pnPixels) * dfPixelSizeX;
-    if (bIsGeographicCoords &&
+    if (bIsGeographicCoordsDeg &&
         ((dfMaxXOut <= 180 && dfMaxXOutNew > 180) || dfMaxXOut == 180))
     {
         dfMaxXOut = 180;
@@ -1206,7 +1208,7 @@ CPLErr CPL_STDCALL GDALSuggestedWarpOutput2(GDALDatasetH hSrcDS,
     }
 
     double dfMinYOutNew = dfMaxYOut - (*pnLines) * dfPixelSizeY;
-    if (bIsGeographicCoords && dfMinYOut >= -90 && dfMinYOutNew < -90)
+    if (bIsGeographicCoordsDeg && dfMinYOut >= -90 && dfMinYOutNew < -90)
     {
         dfMinYOut = -90;
         dfPixelSizeY = (dfMaxYOut - dfMinYOut) / *pnLines;
@@ -1426,8 +1428,11 @@ static void InsertCenterLong(GDALDatasetH hDS, OGRSpatialReference *poSRS,
                              CPLStringList &aosOptions)
 
 {
-    if (!poSRS->IsGeographic())
+    if (!poSRS->IsGeographic() || std::fabs(poSRS->GetAngularUnits() -
+                                            CPLAtof(SRS_UA_DEGREE_CONV)) > 1e-9)
+    {
         return;
+    }
 
     if (poSRS->GetExtension(nullptr, "CENTER_LONG"))
         return;
@@ -2967,7 +2972,9 @@ int GDALTransformLonLatToDestGenImgProjTransformer(void *hTransformArg,
     double z = 0;
     int success = true;
     auto poSourceCRS = psReprojInfo->poForwardTransform->GetSourceCS();
-    if (poSourceCRS->IsGeographic())
+    if (poSourceCRS->IsGeographic() &&
+        std::fabs(poSourceCRS->GetAngularUnits() -
+                  CPLAtof(SRS_UA_DEGREE_CONV)) < 1e-9)
     {
         // Optimization to avoid creating a OGRCoordinateTransformation
         OGRAxisOrientation eSourceFirstAxisOrient = OAO_Other;
diff --git a/autotest/utilities/data/geog_arc_second.tif b/autotest/utilities/data/geog_arc_second.tif
new file mode 100644
index 0000000000000000000000000000000000000000..19c84cabfdbcde6d3de9611398cdaac9fc238804
GIT binary patch
literal 587
zcmebD)MDUZU|<krU|?isU<9(B0SQJZTMx))g0exXn4xTtIyNM>AR`NdFi;F+fGCo<
z7?dplloAKBL27ORu?&)$U?etQGY<nZP~Hb<PE!jH0~?T?0%WgmXJ*g<vZnyq8`>di
zo&eb!+nE^rf#QFF?8YTb3@kvlk*O(=90-uO;%tSQ)K&XSXYFK;3S6}(<c7D4>!UTV
z2_AfQb7j7f;ZFycdPZ<Oa58|r!^ptI01R9(4T>m6hX3`992?vHfQp&G{B06!V0wZw
zC&z|%kX~*u&A`A5rkNNza!aPkE-3<;@31T+!VzYVyK}IDUt(^mf~`V$eo;<}f_rLy
zN@_`FvO;iWaY<^ff}y2_NlideW-gFztKeRgnwnRhnVeDM7@C(^QVbMIEK1fbPEF3w
JOR1@0002LOND%-4

literal 0
HcmV?d00001

diff --git a/autotest/utilities/test_gdalwarp_lib.py b/autotest/utilities/test_gdalwarp_lib.py
index 8341b83c93bb..3cb97c04aee6 100755
--- a/autotest/utilities/test_gdalwarp_lib.py
+++ b/autotest/utilities/test_gdalwarp_lib.py
@@ -4314,3 +4314,19 @@ def test_gdalwarp_lib_average_ten_ten_to_one_one():
         "", src_ds, width=1, height=1, resampleAlg=gdal.GRIORA_Average, format="MEM"
     )
     assert out_ds.GetRasterBand(1).ComputeRasterMinMax() == (1, 1)
+
+
+###############################################################################
+# Test bugfix for https://github.com/OSGeo/gdal/issues/10975
+
+
+def test_gdalwarp_lib_src_is_geog_arc_second():
+
+    out_ds = gdal.Warp(
+        "",
+        "data/geog_arc_second.tif",
+        options="-f MEM -srcnodata 0 -te -81.5 28.5 -81.0 29.0 -t_srs EPSG:4326",
+    )
+    assert out_ds.RasterXSize == 5464
+    assert out_ds.RasterYSize == 5464
+    assert out_ds.GetRasterBand(1).Checksum() == 31856