gdalwarp: add arbitrary percentile resampling (Pxx)

alg/gdalwarper.cpp

@@ -2486,18 +2486,34 @@ bool GDALGetWarpResampleAlg(const char *pszResampling,
         eResampleAlg = GRA_Q3;
     else if (EQUAL(pszResampling, "sum"))
         eResampleAlg = GRA_Sum;
-    else
+    else if ((pszResampling[0] == 'P' || pszResampling[0] == 'p') &&
+            CPLStrnlen(pszResampling, 16) > 1)
     {
-        if (bThrow)
-        {
-            throw std::invalid_argument("Unknown resampling method");
-        }
-        else
+        const char* pszValue = pszResampling + 1;
+        char* pszEnd = nullptr;
+        const double dfPercentile = CPLStrtod(pszValue, &pszEnd);
+
+        if (pszEnd == pszValue || *pszEnd != '\0' ||
+            dfPercentile <= 0.0 || dfPercentile >= 100.0)
         {
-            CPLError(CE_Failure, CPLE_IllegalArg,
-                     "Unknown resampling method: %s.", pszResampling);
-            return false;
+            if (bThrow)
+            {
+                throw std::invalid_argument(
+                    "Invalid percentile resampling method");
+            }
+            else
+            {
+                CPLError(CE_Failure, CPLE_IllegalArg,
+                        "Invalid percentile resampling method: %s.",
+                        pszResampling);
+                return false;
+            }
         }
+
+        eResampleAlg = GRA_Percentile;
+
+        // Store percentile for warp kernel
+        CPLSetConfigOption("GDAL_WARP_PERCENTILE",
+                        CPLSPrintf("%.17g", dfPercentile));
     }
-    return true;
 }

alg/gdalwarper.h

@@ -63,8 +63,10 @@ typedef enum
     /*! RMS (weighted root mean square (quadratic mean) of all non-NODATA
        contributing pixels) */
     GRA_RMS = 14,
+    /*! Percentile (select arbitrary percentile of all non-NODATA contributing pixels) */
+    GRA_Percentile = 15,
     /*! @cond Doxygen_Suppress */
-    GRA_LAST_VALUE = GRA_RMS
+    GRA_LAST_VALUE = GRA_Percentile
     /*! @endcond */
 } GDALResampleAlg;
 

alg/gdalwarpkernel.cpp

@@ -7337,15 +7337,40 @@ static void GWKAverageOrModeThread(void *pData)
     {
         quant = 0.75f;
     }
+    else if (poWK->eResample == GRA_Percentile)
+    {
+        const char* pszPercentile =
+            CSLFetchNameValue(poWK->papszWarpOptions, "PERCENTILE");
+
+        if (!pszPercentile)
+        {
+            CPLError(CE_Fatal, CPLE_AppDefined,
+                    "Percentile resampling selected but "
+                    "PERCENTILE warping option is not set");
+            return;
+        }
+
+        const double dfPercentile = CPLAtof(pszPercentile);
+        if (dfPercentile <= 0.0 || dfPercentile >= 100.0)
+        {
+            CPLError(CE_Fatal, CPLE_AppDefined,
+                    "Invalid PERCENTILE value: %s",
+                    pszPercentile);
+            return;
+        }
+
+        quant = static_cast<float>(dfPercentile / 100.0);
+    }
     else if (poWK->eResample != GRA_Average && poWK->eResample != GRA_RMS &&
-             poWK->eResample != GRA_Min && poWK->eResample != GRA_Max)
+            poWK->eResample != GRA_Min && poWK->eResample != GRA_Max)
     {
         // Other resample algorithms not permitted here.
         CPLError(CE_Fatal, CPLE_AppDefined,
-                 "GDALWarpKernel():GWKAverageOrModeThread() ERROR, "
-                 "illegal resample");
+                "GDALWarpKernel():GWKAverageOrModeThread() ERROR, "
+                "illegal resample");
     }
 
+
     CPLDebug("GDAL", "GDALWarpKernel():GWKAverageOrModeThread()");
 
     /* -------------------------------------------------------------------- */

autotest/utilities/test_gdalwarp.py

@@ -1606,3 +1606,83 @@ def test_gdalwarp_invalid_wm(gdalwarp_path, tmp_vsimem):
     )
     assert "Memory percentage" in err
     assert "Failed to parse value of -wm" in err
+
+def test_gdalwarp_percentile_equivalence(gdalwarp_path, tmp_vsimem):
+    """Test that P50, P25, P75 match med, q1, q3"""
+
+    src = f"{tmp_vsimem}/src.tif"
+    med = f"{tmp_vsimem}/med.tif"
+    p50 = f"{tmp_vsimem}/p50.tif"
+    q1 = f"{tmp_vsimem}/q1.tif"
+    p25 = f"{tmp_vsimem}/p25.tif"
+    q3 = f"{tmp_vsimem}/q3.tif"
+    p75 = f"{tmp_vsimem}/p75.tif"
+
+    # Create a simple raster
+    gdal.Translate(
+        src,
+        "../gcore/data/byte.tif",
+        width=5,
+        height=1,
+        resampleAlg="nearest",
+    )
+
+    ds = gdal.Open(src, gdal.GA_Update)
+    ds.GetRasterBand(1).WriteArray([[1, 2, 3, 4, 100]])
+    ds = None
+
+    gdaltest.runexternal(f"{gdalwarp_path} -r med {src} {med}")
+    gdaltest.runexternal(f"{gdalwarp_path} -r P50 {src} {p50}")
+
+    gdaltest.runexternal(f"{gdalwarp_path} -r q1 {src} {q1}")
+    gdaltest.runexternal(f"{gdalwarp_path} -r P25 {src} {p25}")
+
+    gdaltest.runexternal(f"{gdalwarp_path} -r q3 {src} {q3}")
+    gdaltest.runexternal(f"{gdalwarp_path} -r P75 {src} {p75}")
+
+    def read_val(path):
+        return gdal.Open(path).ReadAsArray()[0][0]
+
+    assert read_val(med) == read_val(p50)
+    assert read_val(q1) == read_val(p25)
+    assert read_val(q3) == read_val(p75)
+
+
+def test_gdalwarp_percentile_outlier(gdalwarp_path, tmp_vsimem):
+    """Test that P95 ignores a strong outlier"""
+
+    src = f"{tmp_vsimem}/src.tif"
+    out = f"{tmp_vsimem}/out.tif"
+
+    gdal.Translate(
+        src,
+        "../gcore/data/byte.tif",
+        width=5,
+        height=1,
+        resampleAlg="nearest",
+    )
+
+    ds = gdal.Open(src, gdal.GA_Update)
+    ds.GetRasterBand(1).WriteArray([[10, 11, 12, 13, 1000]])
+    ds = None
+
+    gdaltest.runexternal(f"{gdalwarp_path} -r P95 {src} {out}")
+
+    val = gdal.Open(out).ReadAsArray()[0][0]
+    assert val < 1000
+
+
+def test_gdalwarp_invalid_percentile(gdalwarp_path, tmp_vsimem):
+    """Test invalid percentile values"""
+
+    ret, err = gdaltest.runexternal_out_and_err(
+        f"{gdalwarp_path} -r P0 ../gcore/data/byte.tif {tmp_vsimem}/out.tif"
+    )
+    assert ret != 0
+    assert "Invalid percentile" in err
+
+    ret, err = gdaltest.runexternal_out_and_err(
+        f"{gdalwarp_path} -r P100 ../gcore/data/byte.tif {tmp_vsimem}/out.tif"
+    )
+    assert ret != 0
+    assert "Invalid percentile" in err

doc/source/programs/gdal_raster_calc.rst

@@ -84,6 +84,27 @@ Program-Specific Options
     Arguments to functions are passed within parentheses, like ``sum(k=1)``,
     ``min(propagateNoData=true)`` or ``interpolate_linear(t0=1,dt=1,t=10)``.
 
+    Built-in variables
+    ^^^^^^^^^^^^^^^^^^
+
+    The following built-in variables are available in expressions:
+
+    - ``_CENTER_X_``: X coordinate of the pixel center, expressed in the dataset CRS
+    - ``_CENTER_Y_``: Y coordinate of the pixel center, expressed in the dataset CRS
+
+    These variables are useful for calculations that depend on spatial position,
+    such as latitude-based corrections, solar angle approximations, or zonal masks.
+
+    Example:
+
+    ::
+
+        gdal raster calc -i A=input.tif \
+            --calc="sin(_CENTER_Y_ * 0.0174533)" \
+            -o output.tif
+
+    Note: To work with longitude/latitude values, the input dataset must be in a
+    geographic coordinate reference system (for example, EPSG:4326).
 
 .. option:: --dialect muparser|builtin