VRT Pixel Functions: Add function to evaluate arbitrary expression

.github/workflows/alpine/Dockerfile.ci

@@ -41,6 +41,7 @@ RUN apk add \
     lz4-dev \
     make \
     mariadb-connector-c-dev \
+    muparser-dev \
     netcdf-dev \
     odbc-cpp-wrapper-dev \
     ogdi-dev \

.github/workflows/cmake_builds.yml

@@ -328,7 +328,7 @@ jobs:
             base-devel git mingw-w64-x86_64-toolchain mingw-w64-x86_64-cmake mingw-w64-x86_64-ccache
             mingw-w64-x86_64-pcre mingw-w64-x86_64-xerces-c mingw-w64-x86_64-zstd mingw-w64-x86_64-libarchive
             mingw-w64-x86_64-geos mingw-w64-x86_64-libspatialite mingw-w64-x86_64-proj
-            mingw-w64-x86_64-cgal mingw-w64-x86_64-libfreexl mingw-w64-x86_64-hdf5 mingw-w64-x86_64-netcdf mingw-w64-x86_64-poppler  mingw-w64-x86_64-podofo mingw-w64-x86_64-postgresql
+            mingw-w64-x86_64-cgal mingw-w64-x86_64-libfreexl mingw-w64-x86_64-hdf5 mingw-w64-x86_64-muparser mingw-w64-x86_64-netcdf mingw-w64-x86_64-poppler  mingw-w64-x86_64-podofo mingw-w64-x86_64-postgresql
             mingw-w64-x86_64-libgeotiff mingw-w64-x86_64-libpng mingw-w64-x86_64-libtiff mingw-w64-x86_64-openjpeg2
             mingw-w64-x86_64-python-pip mingw-w64-x86_64-python-numpy mingw-w64-x86_64-python-pytest mingw-w64-x86_64-python-setuptools mingw-w64-x86_64-python-lxml mingw-w64-x86_64-swig mingw-w64-x86_64-python-psutil mingw-w64-x86_64-blosc mingw-w64-x86_64-libavif
       - name: Setup cache
@@ -434,7 +434,7 @@ jobs:
             cfitsio freexl geotiff libjpeg-turbo libpq libspatialite libwebp-base pcre pcre2 postgresql \
             sqlite tiledb zstd cryptopp cgal doxygen librttopo openssl liblzma-devel \
             openjdk ant qhull armadillo blas blas-devel libblas libcblas liblapack liblapacke blosc libarchive \
-            libarrow pyarrow libaec libheif libavif cmake fsspec
+            libarrow pyarrow libaec libheif libavif muparser cmake fsspec
     - name: Check CMake version
       shell: bash -l {0}
       run: |

.github/workflows/fedora_rawhide/Dockerfile.ci

@@ -19,6 +19,7 @@ RUN dnf install -y clang make diffutils ccache cmake \
       mdbtools-devel mdbtools-odbc unixODBC-devel \
       armadillo-devel qhull-devel \
       hdf-devel hdf5-devel netcdf-devel \
+      muParser-devel \
       libpq-devel \
       libavif-devel \
       python3-setuptools python3-pip python3-devel python3-lxml swig \

.github/workflows/ubuntu_20.04/Dockerfile.ci

@@ -49,6 +49,7 @@ RUN apt-get update -y \
     liblz4-dev \
     liblzma-dev \
     libmono-system-drawing4.0-cil \
+    libmuparser-dev \
     libmysqlclient-dev \
     libnetcdf-dev \
     libogdi-dev \
@@ -274,6 +275,11 @@ RUN if test "${OPENDRIVE_VERSION}" != ""; then ( \
     && rm -rf libOpenDRIVE-${OPENDRIVE_VERSION} \
     ); fi
 
+# Install exprtk
+RUN wget -q https://www.partow.net/downloads/exprtk.zip && \
+    unzip -j -d /usr/local/include exprtk.zip exprtk/exprtk.hpp && \
+    rm exprtk.zip
+
 RUN ldconfig
 
 COPY requirements.txt /tmp/

.github/workflows/ubuntu_20.04/build.sh

@@ -15,6 +15,7 @@ cmake "${GDAL_SOURCE_DIR:=..}" \
     -DCMAKE_INSTALL_PREFIX=/tmp/install-gdal \
     -DGDAL_USE_TIFF_INTERNAL=OFF \
     -DGDAL_USE_GEOTIFF_INTERNAL=OFF \
+    -DGDAL_USE_EXPRTK=ON \
     -DECW_ROOT=/opt/libecwj2-3.3 \
     -DMRSID_ROOT=/usr/local \
     -DFileGDB_ROOT=/usr/local/FileGDB_API \

.github/workflows/ubuntu_24.04/Dockerfile.ci

@@ -37,6 +37,7 @@ RUN apt-get update && \
     liblcms2-2 \
     liblz4-dev \
     liblzma-dev \
+    libmuparser-dev \
     libmysqlclient-dev \
     libnetcdf-dev \
     libogdi-dev \

.pre-commit-config.yaml

@@ -42,7 +42,9 @@ repos:
                 autotest/cpp/data/|
                 autotest/gdrivers/data/|
                 swig/|
+                third_party/exprtk/|
                 third_party/fast_float/|
+                third_party/muparser/|
                 third_party/LercLib/|
                 autotest/ogr/data/|
                 alg/internal_libqhull/|

autotest/gdrivers/data/vrt/arraysource_derived_expression.vrt

@@ -0,0 +1,14 @@
+<VRTDataset rasterXSize="20" rasterYSize="20">
+  <VRTRasterBand dataType="Float64" band="1" subClass="VRTDerivedRasterBand">
+    <PixelFunctionType>expression</PixelFunctionType>
+    <PixelFunctionArguments expression="B1 + 5" dialect="muparser"/>
+    <ArraySource>
+        <Array name="test">
+            <DataType>Float64</DataType>
+            <Dimension name="Y" size="20"/>
+            <Dimension name="X" size="20"/>
+            <ConstantValue>10</ConstantValue>
+        </Array>
+    </ArraySource>
+  </VRTRasterBand>
+</VRTDataset>

autotest/gdrivers/vrtderived.py

@@ -1055,6 +1055,225 @@ def identity(in_ar, out_ar, *args, **kwargs):
             assert vrt_ds.GetRasterBand(1).DataType == dtype
 
 
+###############################################################################
+# Test arbitrary expression pixel functions
+
+
+def vrt_expression_xml(tmpdir, expression, dialect, sources):
+
+    drv = gdal.GetDriverByName("GTiff")
+
+    nx = 1
+    ny = 1
+
+    expression = expression.replace("<", "&lt;").replace(">", "&gt;")
+
+    xml = f"""<VRTDataset rasterXSize="{nx}" rasterYSize="{ny}">
+              <VRTRasterBand dataType="Float64" band="1" subClass="VRTDerivedRasterBand">
+                 <PixelFunctionType>expression</PixelFunctionType>
+                 <PixelFunctionArguments expression="{expression}" dialect="{dialect}" />"""
+
+    for i, source in enumerate(sources):
+        if type(source) is tuple:
+            source_name, source_value = source
+        else:
+            source_name = ""
+            source_value = source
+
+        src_fname = tmpdir / f"source_{i}.tif"
+
+        with drv.Create(src_fname, 1, 1, 1, gdal.GDT_Float64) as ds:
+            ds.GetRasterBand(1).Fill(source_value)
+
+        xml += f"""<SimpleSource name="{source_name}">
+                     <SourceFilename relativeToVRT="0">{src_fname}</SourceFilename>
+                     <SourceBand>1</SourceBand>
+                   </SimpleSource>"""
+
+    xml += "</VRTRasterBand></VRTDataset>"
+
+    return xml
+
+
+@pytest.mark.parametrize(
+    "expression,sources,result,dialects",
+    [
+        pytest.param("A", [("A", 77)], 77, None, id="identity"),
+        pytest.param(
+            "(NIR-R)/(NIR+R)",
+            [("NIR", 77), ("R", 63)],
+            (77 - 63) / (77 + 63),
+            None,
+            id="simple expression",
+        ),
+        pytest.param(
+            "if (A > B) 1.5*C ; else A",
+            [("A", 77), ("B", 63), ("C", 18)],
+            27,
+            ["exprtk"],
+            id="exprtk conditional (explicit)",
+        ),
+        pytest.param(
+            "(A > B) ? 1.5*C : A",
+            [("A", 77), ("B", 63), ("C", 18)],
+            27,
+            ["muparser"],
+            id="muparser conditional (explicit)",
+        ),
+        pytest.param(
+            "(A > B)*(1.5*C) + (A <= B)*(A)",
+            [("A", 77), ("B", 63), ("C", 18)],
+            27,
+            None,
+            id="conditional (implicit)",
+        ),
+        pytest.param(
+            "B2 * PopDensity",
+            [("PopDensity", 3), ("", 7)],
+            21,
+            None,
+            id="implicit source name",
+        ),
+        pytest.param(
+            "B1 / sum(BANDS)",
+            [("", 3), ("", 5), ("", 31)],
+            3 / (3 + 5 + 31),
+            None,
+            id="use of BANDS variable",
+        ),
+        pytest.param(
+            "B1 / sum(B2, B3) ",
+            [("", 3), ("", 5), ("", 31)],
+            3 / (5 + 31),
+            None,
+            id="aggregate specified inputs",
+        ),
+        pytest.param(
+            "var q[2] := {B2, B3}; B1 * q",
+            [("", 3), ("", 5), ("", 31)],
+            15,  # First value in returned vector. This behavior doesn't seem desirable
+            # but I haven't figured out how to detect a vector return.
+            ["exprtk"],
+            id="return vector",
+        ),
+        pytest.param(
+            "B1 + B2 + B3",
+            (5, 9, float("nan")),
+            float("nan"),
+            None,
+            id="nan propagated via arithmetic",
+        ),
+        pytest.param(
+            "if (B3) B1 ; else B2",
+            (5, 9, float("nan")),
+            5,
+            ["exprtk"],
+            id="exprtk nan = truth in conditional?",
+        ),
+        pytest.param(
+            "B3 ? B1 : B2",
+            (5, 9, float("nan")),
+            5,
+            ["muparser"],
+            id="muparser nan = truth in conditional?",
+        ),
+        pytest.param(
+            "if (B3 > 0) B1 ; else B2",
+            (5, 9, float("nan")),
+            9,
+            ["exprtk"],
+            id="exprtk nan comparison is false in conditional",
+        ),
+        pytest.param(
+            "(B3 > 0) ? B1 : B2",
+            (5, 9, float("nan")),
+            9,
+            ["muparser"],
+            id="muparser nan comparison is false in conditional",
+        ),
+        pytest.param(
+            "if (B1 > 5) B1",
+            (1,),
+            float("nan"),
+            ["exprtk"],
+            id="expression returns nodata",
+        ),
+    ],
+)
+@pytest.mark.parametrize("dialect", ("exprtk", "muparser"))
+def test_vrt_pixelfn_expression(
+    tmp_vsimem, expression, sources, result, dialect, dialects
+):
+    pytest.importorskip("numpy")
+
+    if not gdaltest.gdal_has_vrt_expression_dialect(dialect):
+        pytest.skip(f"Expression dialect {dialect} is not available")
+
+    if dialects and dialect not in dialects:
+        pytest.skip(f"Expression not supported for dialect {dialect}")
+
+    xml = vrt_expression_xml(tmp_vsimem, expression, dialect, sources)
+
+    with gdal.Open(xml) as ds:
+        assert pytest.approx(ds.ReadAsArray()[0][0], nan_ok=True) == result
+
+
+@pytest.mark.parametrize(
+    "expression,sources,dialect,exception",
+    [
+        pytest.param(
+            "A*B + C",
+            [("A", 77), ("B", 63)],
+            "exprtk",
+            "Undefined symbol",
+            id="exprtk undefined variable",
+        ),
+        pytest.param(
+            "A*B + C",
+            [("A", 77), ("B", 63)],
+            "muparser",
+            "Unexpected token",
+            id="muparser undefined variable",
+        ),
+        pytest.param(
+            "(".join(["asin", "sin", "acos", "cos"] * 100) + "(X" + 100 * 4 * ")",
+            [("X", 0.5)],
+            "exprtk",
+            "exceeds maximum allowed stack depth",
+            id="expression is too complex",
+        ),
+        pytest.param(
+            " ".join(["sin(x) + cos(x)"] * 10000),
+            [("x", 0.5)],
+            "exprtk",
+            "exceeds maximum of 100000 set by GDAL_EXPRTK_MAX_EXPRESSION_LENGTH",
+            id="expression is too long",
+        ),
+    ],
+)
+def test_vrt_pixelfn_expression_invalid(
+    tmp_vsimem, expression, sources, dialect, exception
+):
+    pytest.importorskip("numpy")
+
+    if not gdaltest.gdal_has_vrt_expression_dialect(dialect):
+        pytest.skip(f"Expression dialect {dialect} is not available")
+
+    messages = []
+
+    def handle(ecls, ecode, emsg):
+        messages.append(emsg)
+
+    xml = vrt_expression_xml(tmp_vsimem, expression, dialect, sources)
+
+    with gdaltest.error_handler(handle):
+        ds = gdal.Open(xml)
+        if ds:
+            assert ds.ReadAsArray() is None
+
+    assert exception in "".join(messages)
+
+
 ###############################################################################
 # Cleanup.