LCOV - code coverage report
Current view: top level - frmts/vrt - vrtderivedrasterband.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 667 716 93.2 %
Date: 2026-01-31 22:56:34 Functions: 30 32 93.8 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  Virtual GDAL Datasets
       4             :  * Purpose:  Implementation of a sourced raster band that derives its raster
       5             :  *           by applying an algorithm (GDALDerivedPixelFunc) to the sources.
       6             :  * Author:   Pete Nagy
       7             :  *
       8             :  ******************************************************************************
       9             :  * Copyright (c) 2005 Vexcel Corp.
      10             :  * Copyright (c) 2008-2011, Even Rouault <even dot rouault at spatialys.com>
      11             :  *
      12             :  * SPDX-License-Identifier: MIT
      13             :  *****************************************************************************/
      14             : 
      15             : #include "cpl_minixml.h"
      16             : #include "cpl_string.h"
      17             : #include "gdal_priv.h"
      18             : #include "vrtdataset.h"
      19             : #include "cpl_multiproc.h"
      20             : #include "gdalpython.h"
      21             : #include "gdalantirecursion.h"
      22             : 
      23             : #include <algorithm>
      24             : #include <array>
      25             : #include <map>
      26             : #include <vector>
      27             : #include <utility>
      28             : 
      29             : /*! @cond Doxygen_Suppress */
      30             : 
      31             : using namespace GDALPy;
      32             : 
      33             : // #define GDAL_VRT_DISABLE_PYTHON
      34             : 
      35             : #ifndef GDAL_VRT_ENABLE_PYTHON_DEFAULT
      36             : // Can be YES, NO or TRUSTED_MODULES
      37             : #define GDAL_VRT_ENABLE_PYTHON_DEFAULT "TRUSTED_MODULES"
      38             : #endif
      39             : 
      40             : /* Flags for getting buffers */
      41             : #define PyBUF_WRITABLE 0x0001
      42             : #define PyBUF_FORMAT 0x0004
      43             : #define PyBUF_ND 0x0008
      44             : #define PyBUF_STRIDES (0x0010 | PyBUF_ND)
      45             : #define PyBUF_INDIRECT (0x0100 | PyBUF_STRIDES)
      46             : #define PyBUF_FULL (PyBUF_INDIRECT | PyBUF_WRITABLE | PyBUF_FORMAT)
      47             : 
      48             : /************************************************************************/
      49             : /*                        GDALCreateNumpyArray()                        */
      50             : /************************************************************************/
      51             : 
      52         557 : static PyObject *GDALCreateNumpyArray(PyObject *pCreateArray, void *pBuffer,
      53             :                                       GDALDataType eType, int nHeight,
      54             :                                       int nWidth)
      55             : {
      56             :     PyObject *poPyBuffer;
      57             :     const size_t nSize =
      58         557 :         static_cast<size_t>(nHeight) * nWidth * GDALGetDataTypeSizeBytes(eType);
      59             :     Py_buffer pybuffer;
      60         557 :     if (PyBuffer_FillInfo(&pybuffer, nullptr, static_cast<char *>(pBuffer),
      61         557 :                           nSize, 0, PyBUF_FULL) != 0)
      62             :     {
      63           0 :         return nullptr;
      64             :     }
      65         557 :     poPyBuffer = PyMemoryView_FromBuffer(&pybuffer);
      66         557 :     PyObject *pArgsCreateArray = PyTuple_New(4);
      67         557 :     PyTuple_SetItem(pArgsCreateArray, 0, poPyBuffer);
      68         557 :     const char *pszDataType = nullptr;
      69         557 :     switch (eType)
      70             :     {
      71         436 :         case GDT_UInt8:
      72         436 :             pszDataType = "uint8";
      73         436 :             break;
      74           3 :         case GDT_Int8:
      75           3 :             pszDataType = "int8";
      76           3 :             break;
      77           2 :         case GDT_UInt16:
      78           2 :             pszDataType = "uint16";
      79           2 :             break;
      80          88 :         case GDT_Int16:
      81          88 :             pszDataType = "int16";
      82          88 :             break;
      83           2 :         case GDT_UInt32:
      84           2 :             pszDataType = "uint32";
      85           2 :             break;
      86           4 :         case GDT_Int32:
      87           4 :             pszDataType = "int32";
      88           4 :             break;
      89           4 :         case GDT_Int64:
      90           4 :             pszDataType = "int64";
      91           4 :             break;
      92           2 :         case GDT_UInt64:
      93           2 :             pszDataType = "uint64";
      94           2 :             break;
      95           2 :         case GDT_Float16:
      96           2 :             pszDataType = "float16";
      97           2 :             break;
      98           4 :         case GDT_Float32:
      99           4 :             pszDataType = "float32";
     100           4 :             break;
     101           4 :         case GDT_Float64:
     102           4 :             pszDataType = "float64";
     103           4 :             break;
     104           0 :         case GDT_CInt16:
     105             :         case GDT_CInt32:
     106           0 :             CPLAssert(FALSE);
     107             :             break;
     108           0 :         case GDT_CFloat16:
     109           0 :             CPLAssert(FALSE);
     110             :             break;
     111           3 :         case GDT_CFloat32:
     112           3 :             pszDataType = "complex64";
     113           3 :             break;
     114           3 :         case GDT_CFloat64:
     115           3 :             pszDataType = "complex128";
     116           3 :             break;
     117           0 :         case GDT_Unknown:
     118             :         case GDT_TypeCount:
     119           0 :             CPLAssert(FALSE);
     120             :             break;
     121             :     }
     122         557 :     PyTuple_SetItem(
     123             :         pArgsCreateArray, 1,
     124             :         PyBytes_FromStringAndSize(pszDataType, strlen(pszDataType)));
     125         557 :     PyTuple_SetItem(pArgsCreateArray, 2, PyLong_FromLong(nHeight));
     126         557 :     PyTuple_SetItem(pArgsCreateArray, 3, PyLong_FromLong(nWidth));
     127             :     PyObject *poNumpyArray =
     128         557 :         PyObject_Call(pCreateArray, pArgsCreateArray, nullptr);
     129         557 :     Py_DecRef(pArgsCreateArray);
     130         557 :     if (PyErr_Occurred())
     131           0 :         PyErr_Print();
     132         557 :     return poNumpyArray;
     133             : }
     134             : 
     135             : /************************************************************************/
     136             : /* ==================================================================== */
     137             : /*                     VRTDerivedRasterBandPrivateData                  */
     138             : /* ==================================================================== */
     139             : /************************************************************************/
     140             : 
     141             : class VRTDerivedRasterBandPrivateData
     142             : {
     143             :     VRTDerivedRasterBandPrivateData(const VRTDerivedRasterBandPrivateData &) =
     144             :         delete;
     145             :     VRTDerivedRasterBandPrivateData &
     146             :     operator=(const VRTDerivedRasterBandPrivateData &) = delete;
     147             : 
     148             :   public:
     149             :     CPLString m_osCode{};
     150             :     CPLString m_osLanguage = "C";
     151             :     int m_nBufferRadius = 0;
     152             :     PyObject *m_poGDALCreateNumpyArray = nullptr;
     153             :     PyObject *m_poUserFunction = nullptr;
     154             :     bool m_bPythonInitializationDone = false;
     155             :     bool m_bPythonInitializationSuccess = false;
     156             :     bool m_bExclusiveLock = false;
     157             :     bool m_bFirstTime = true;
     158             :     std::vector<std::pair<CPLString, CPLString>> m_oFunctionArgs{};
     159             :     bool m_bSkipNonContributingSourcesSpecified = false;
     160             :     bool m_bSkipNonContributingSources = false;
     161             :     GIntBig m_nAllowedRAMUsage = 0;
     162             : 
     163        2252 :     VRTDerivedRasterBandPrivateData()
     164        2252 :         : m_nAllowedRAMUsage(CPLGetUsablePhysicalRAM() / 10 * 4)
     165             :     {
     166             :         // Use only up to 40% of RAM to acquire source bands and generate the
     167             :         // output buffer.
     168             :         // Only for tests now
     169        2252 :         const char *pszMAX_RAM = "VRT_DERIVED_DATASET_ALLOWED_RAM_USAGE";
     170        2252 :         if (const char *pszVal = CPLGetConfigOption(pszMAX_RAM, nullptr))
     171             :         {
     172           1 :             CPL_IGNORE_RET_VAL(
     173           1 :                 CPLParseMemorySize(pszVal, &m_nAllowedRAMUsage, nullptr));
     174             :         }
     175        2252 :     }
     176             : 
     177             :     ~VRTDerivedRasterBandPrivateData();
     178             : };
     179             : 
     180        2252 : VRTDerivedRasterBandPrivateData::~VRTDerivedRasterBandPrivateData()
     181             : {
     182        2252 :     if (m_poGDALCreateNumpyArray)
     183          51 :         Py_DecRef(m_poGDALCreateNumpyArray);
     184        2252 :     if (m_poUserFunction)
     185          52 :         Py_DecRef(m_poUserFunction);
     186        2252 : }
     187             : 
     188             : /************************************************************************/
     189             : /* ==================================================================== */
     190             : /*                          VRTDerivedRasterBand                        */
     191             : /* ==================================================================== */
     192             : /************************************************************************/
     193             : 
     194             : /************************************************************************/
     195             : /*                        VRTDerivedRasterBand()                        */
     196             : /************************************************************************/
     197             : 
     198        1490 : VRTDerivedRasterBand::VRTDerivedRasterBand(GDALDataset *poDSIn, int nBandIn)
     199             :     : VRTSourcedRasterBand(poDSIn, nBandIn), m_poPrivate(nullptr),
     200        1490 :       eSourceTransferType(GDT_Unknown)
     201             : {
     202        1490 :     m_poPrivate = new VRTDerivedRasterBandPrivateData;
     203        1490 : }
     204             : 
     205             : /************************************************************************/
     206             : /*                        VRTDerivedRasterBand()                        */
     207             : /************************************************************************/
     208             : 
     209         762 : VRTDerivedRasterBand::VRTDerivedRasterBand(GDALDataset *poDSIn, int nBandIn,
     210             :                                            GDALDataType eType, int nXSize,
     211             :                                            int nYSize, int nBlockXSizeIn,
     212         762 :                                            int nBlockYSizeIn)
     213             :     : VRTSourcedRasterBand(poDSIn, nBandIn, eType, nXSize, nYSize,
     214             :                            nBlockXSizeIn, nBlockYSizeIn),
     215         762 :       m_poPrivate(nullptr), eSourceTransferType(GDT_Unknown)
     216             : {
     217         762 :     m_poPrivate = new VRTDerivedRasterBandPrivateData;
     218         762 : }
     219             : 
     220             : /************************************************************************/
     221             : /*                       ~VRTDerivedRasterBand()                        */
     222             : /************************************************************************/
     223             : 
     224        4504 : VRTDerivedRasterBand::~VRTDerivedRasterBand()
     225             : 
     226             : {
     227        2252 :     delete m_poPrivate;
     228        4504 : }
     229             : 
     230             : /************************************************************************/
     231             : /*                              Cleanup()                               */
     232             : /************************************************************************/
     233             : 
     234        1126 : void VRTDerivedRasterBand::Cleanup()
     235             : {
     236        1126 : }
     237             : 
     238             : /************************************************************************/
     239             : /*                     GetGlobalMapPixelFunction()                      */
     240             : /************************************************************************/
     241             : 
     242             : static std::map<std::string,
     243             :                 std::pair<VRTDerivedRasterBand::PixelFunc, std::string>> &
     244       65963 : GetGlobalMapPixelFunction()
     245             : {
     246             :     static std::map<std::string,
     247             :                     std::pair<VRTDerivedRasterBand::PixelFunc, std::string>>
     248       65963 :         gosMapPixelFunction;
     249       65963 :     return gosMapPixelFunction;
     250             : }
     251             : 
     252             : /************************************************************************/
     253             : /*                          AddPixelFunction()                          */
     254             : /************************************************************************/
     255             : 
     256             : /*! @endcond */
     257             : 
     258             : /**
     259             :  * This adds a pixel function to the global list of available pixel
     260             :  * functions for derived bands.  Pixel functions must be registered
     261             :  * in this way before a derived band tries to access data.
     262             :  *
     263             :  * Derived bands are stored with only the name of the pixel function
     264             :  * that it will apply, and if a pixel function matching the name is not
     265             :  * found the IRasterIO() call will do nothing.
     266             :  *
     267             :  * @param pszName Name used to access pixel function
     268             :  * @param pfnNewFunction Pixel function associated with name.  An
     269             :  *  existing pixel function registered with the same name will be
     270             :  *  replaced with the new one.
     271             :  *
     272             :  * @return CE_None, invalid (NULL) parameters are currently ignored.
     273             :  */
     274       15752 : CPLErr CPL_STDCALL GDALAddDerivedBandPixelFunc(
     275             :     const char *pszName, GDALDerivedPixelFunc pfnNewFunction)
     276             : {
     277       15752 :     if (pszName == nullptr || pszName[0] == '\0' || pfnNewFunction == nullptr)
     278             :     {
     279           0 :         return CE_None;
     280             :     }
     281             : 
     282       31504 :     GetGlobalMapPixelFunction()[pszName] = {
     283          54 :         [pfnNewFunction](void **papoSources, int nSources, void *pData,
     284             :                          int nBufXSize, int nBufYSize, GDALDataType eSrcType,
     285             :                          GDALDataType eBufType, int nPixelSpace, int nLineSpace,
     286          54 :                          CSLConstList papszFunctionArgs)
     287             :         {
     288             :             (void)papszFunctionArgs;
     289          54 :             return pfnNewFunction(papoSources, nSources, pData, nBufXSize,
     290             :                                   nBufYSize, eSrcType, eBufType, nPixelSpace,
     291          54 :                                   nLineSpace);
     292             :         },
     293       47256 :         ""};
     294             : 
     295       15752 :     return CE_None;
     296             : }
     297             : 
     298             : /**
     299             :  * This adds a pixel function to the global list of available pixel
     300             :  * functions for derived bands.  Pixel functions must be registered
     301             :  * in this way before a derived band tries to access data.
     302             :  *
     303             :  * Derived bands are stored with only the name of the pixel function
     304             :  * that it will apply, and if a pixel function matching the name is not
     305             :  * found the IRasterIO() call will do nothing.
     306             :  *
     307             :  * @param pszName Name used to access pixel function
     308             :  * @param pfnNewFunction Pixel function associated with name.  An
     309             :  *  existing pixel function registered with the same name will be
     310             :  *  replaced with the new one.
     311             :  * @param pszMetadata Pixel function metadata (not currently implemented)
     312             :  *
     313             :  * @return CE_None, invalid (NULL) parameters are currently ignored.
     314             :  * @since GDAL 3.4
     315             :  */
     316       47252 : CPLErr CPL_STDCALL GDALAddDerivedBandPixelFuncWithArgs(
     317             :     const char *pszName, GDALDerivedPixelFuncWithArgs pfnNewFunction,
     318             :     const char *pszMetadata)
     319             : {
     320       47252 :     if (!pszName || pszName[0] == '\0' || !pfnNewFunction)
     321             :     {
     322           0 :         return CE_None;
     323             :     }
     324             : 
     325       94504 :     GetGlobalMapPixelFunction()[pszName] = {pfnNewFunction,
     326      141756 :                                             pszMetadata ? pszMetadata : ""};
     327             : 
     328       47252 :     return CE_None;
     329             : }
     330             : 
     331             : /*! @cond Doxygen_Suppress */
     332             : 
     333             : /**
     334             :  * This adds a pixel function to the global list of available pixel
     335             :  * functions for derived bands.
     336             :  *
     337             :  * This is the same as the C function GDALAddDerivedBandPixelFunc()
     338             :  *
     339             :  * @param pszFuncNameIn Name used to access pixel function
     340             :  * @param pfnNewFunction Pixel function associated with name.  An
     341             :  *  existing pixel function registered with the same name will be
     342             :  *  replaced with the new one.
     343             :  *
     344             :  * @return CE_None, invalid (NULL) parameters are currently ignored.
     345             :  */
     346             : CPLErr
     347           0 : VRTDerivedRasterBand::AddPixelFunction(const char *pszFuncNameIn,
     348             :                                        GDALDerivedPixelFunc pfnNewFunction)
     349             : {
     350           0 :     return GDALAddDerivedBandPixelFunc(pszFuncNameIn, pfnNewFunction);
     351             : }
     352             : 
     353           0 : CPLErr VRTDerivedRasterBand::AddPixelFunction(
     354             :     const char *pszFuncNameIn, GDALDerivedPixelFuncWithArgs pfnNewFunction,
     355             :     const char *pszMetadata)
     356             : {
     357           0 :     return GDALAddDerivedBandPixelFuncWithArgs(pszFuncNameIn, pfnNewFunction,
     358           0 :                                                pszMetadata);
     359             : }
     360             : 
     361             : /************************************************************************/
     362             : /*                          GetPixelFunction()                          */
     363             : /************************************************************************/
     364             : 
     365             : /**
     366             :  * Get a pixel function previously registered using the global
     367             :  * AddPixelFunction.
     368             :  *
     369             :  * @param pszFuncNameIn The name associated with the pixel function.
     370             :  *
     371             :  * @return A pointer to a std::pair whose first element is the pixel
     372             :  *         function pointer and second element is the pixel function
     373             :  *         metadata string. If no pixel function has been registered
     374             :  *         for pszFuncNameIn, nullptr will be returned.
     375             :  */
     376             : /* static */
     377             : const std::pair<VRTDerivedRasterBand::PixelFunc, std::string> *
     378        2856 : VRTDerivedRasterBand::GetPixelFunction(const char *pszFuncNameIn)
     379             : {
     380        2856 :     if (pszFuncNameIn == nullptr || pszFuncNameIn[0] == '\0')
     381             :     {
     382           0 :         return nullptr;
     383             :     }
     384             : 
     385        2856 :     const auto &oMapPixelFunction = GetGlobalMapPixelFunction();
     386        2856 :     const auto oIter = oMapPixelFunction.find(pszFuncNameIn);
     387             : 
     388        2856 :     if (oIter == oMapPixelFunction.end())
     389           3 :         return nullptr;
     390             : 
     391        2853 :     return &(oIter->second);
     392             : }
     393             : 
     394             : /************************************************************************/
     395             : /*                       GetPixelFunctionNames()                        */
     396             : /************************************************************************/
     397             : 
     398             : /**
     399             :  * Return the list of available pixel function names.
     400             :  */
     401             : /* static */
     402         103 : std::vector<std::string> VRTDerivedRasterBand::GetPixelFunctionNames()
     403             : {
     404         103 :     std::vector<std::string> res;
     405        4223 :     for (const auto &iter : GetGlobalMapPixelFunction())
     406             :     {
     407        4120 :         res.push_back(iter.first);
     408             :     }
     409         103 :     return res;
     410             : }
     411             : 
     412             : /************************************************************************/
     413             : /*                        SetPixelFunctionName()                        */
     414             : /************************************************************************/
     415             : 
     416             : /**
     417             :  * Set the pixel function name to be applied to this derived band.  The
     418             :  * name should match a pixel function registered using AddPixelFunction.
     419             :  *
     420             :  * @param pszFuncNameIn Name of pixel function to be applied to this derived
     421             :  * band.
     422             :  */
     423        2253 : void VRTDerivedRasterBand::SetPixelFunctionName(const char *pszFuncNameIn)
     424             : {
     425        2253 :     osFuncName = (pszFuncNameIn == nullptr) ? "" : pszFuncNameIn;
     426        2253 : }
     427             : 
     428             : /************************************************************************/
     429             : /*                      AddPixelFunctionArgument()                      */
     430             : /************************************************************************/
     431             : 
     432             : /**
     433             :  *  Set a pixel function argument to a specified value.
     434             :  * @param pszArg the argument name
     435             :  * @param pszValue the argument value
     436             :  *
     437             :  * @since 3.12
     438             :  */
     439        2664 : void VRTDerivedRasterBand::AddPixelFunctionArgument(const char *pszArg,
     440             :                                                     const char *pszValue)
     441             : {
     442        2664 :     m_poPrivate->m_oFunctionArgs.emplace_back(pszArg, pszValue);
     443        2664 : }
     444             : 
     445             : /************************************************************************/
     446             : /*                      SetPixelFunctionLanguage()                      */
     447             : /************************************************************************/
     448             : 
     449             : /**
     450             :  * Set the language of the pixel function.
     451             :  *
     452             :  * @param pszLanguage Language of the pixel function (only "C" and "Python"
     453             :  * are supported currently)
     454             :  * @since GDAL 2.3
     455             :  */
     456           1 : void VRTDerivedRasterBand::SetPixelFunctionLanguage(const char *pszLanguage)
     457             : {
     458           1 :     m_poPrivate->m_osLanguage = pszLanguage;
     459           1 : }
     460             : 
     461             : /************************************************************************/
     462             : /*                   SetSkipNonContributingSources()                    */
     463             : /************************************************************************/
     464             : 
     465             : /** Whether sources that do not intersect the VRTRasterBand RasterIO() requested
     466             :  * region should be omitted. By default, data for all sources, including ones
     467             :  * that do not intersect it, are passed to the pixel function. By setting this
     468             :  * parameter to true, only sources that intersect the requested region will be
     469             :  * passed.
     470             :  *
     471             :  * @param bSkip whether to skip non-contributing sources
     472             :  *
     473             :  * @since 3.12
     474             :  */
     475           5 : void VRTDerivedRasterBand::SetSkipNonContributingSources(bool bSkip)
     476             : {
     477           5 :     m_poPrivate->m_bSkipNonContributingSources = bSkip;
     478           5 :     m_poPrivate->m_bSkipNonContributingSourcesSpecified = true;
     479           5 : }
     480             : 
     481             : /************************************************************************/
     482             : /*                       SetSourceTransferType()                        */
     483             : /************************************************************************/
     484             : 
     485             : /**
     486             :  * Set the transfer type to be used to obtain pixel information from
     487             :  * all of the sources.  If unset, the transfer type used will be the
     488             :  * same as the derived band data type.  This makes it possible, for
     489             :  * example, to pass CFloat32 source pixels to the pixel function, even
     490             :  * if the pixel function generates a raster for a derived band that
     491             :  * is of type Byte.
     492             :  *
     493             :  * @param eDataTypeIn Data type to use to obtain pixel information from
     494             :  * the sources to be passed to the derived band pixel function.
     495             :  */
     496          21 : void VRTDerivedRasterBand::SetSourceTransferType(GDALDataType eDataTypeIn)
     497             : {
     498          21 :     eSourceTransferType = eDataTypeIn;
     499          21 : }
     500             : 
     501             : /************************************************************************/
     502             : /*                          InitializePython()                          */
     503             : /************************************************************************/
     504             : 
     505         463 : bool VRTDerivedRasterBand::InitializePython()
     506             : {
     507         463 :     if (m_poPrivate->m_bPythonInitializationDone)
     508         397 :         return m_poPrivate->m_bPythonInitializationSuccess;
     509             : 
     510          66 :     m_poPrivate->m_bPythonInitializationDone = true;
     511          66 :     m_poPrivate->m_bPythonInitializationSuccess = false;
     512             : 
     513          66 :     const size_t nIdxDot = osFuncName.rfind(".");
     514         132 :     CPLString osPythonModule;
     515         132 :     CPLString osPythonFunction;
     516          66 :     if (nIdxDot != std::string::npos)
     517             :     {
     518          29 :         osPythonModule = osFuncName.substr(0, nIdxDot);
     519          29 :         osPythonFunction = osFuncName.substr(nIdxDot + 1);
     520             :     }
     521             :     else
     522             :     {
     523          37 :         osPythonFunction = osFuncName;
     524             :     }
     525             : 
     526             : #ifndef GDAL_VRT_DISABLE_PYTHON
     527             :     const char *pszPythonEnabled =
     528          66 :         CPLGetConfigOption("GDAL_VRT_ENABLE_PYTHON", nullptr);
     529             : #else
     530             :     const char *pszPythonEnabled = "NO";
     531             : #endif
     532             :     const CPLString osPythonEnabled(
     533         132 :         pszPythonEnabled ? pszPythonEnabled : GDAL_VRT_ENABLE_PYTHON_DEFAULT);
     534             : 
     535          66 :     if (EQUAL(osPythonEnabled, "TRUSTED_MODULES"))
     536             :     {
     537          12 :         bool bIsTrustedModule = false;
     538             :         const CPLString osVRTTrustedModules(
     539          12 :             CPLGetConfigOption("GDAL_VRT_PYTHON_TRUSTED_MODULES", ""));
     540          12 :         if (!osPythonModule.empty())
     541             :         {
     542             :             char **papszTrustedModules =
     543          10 :                 CSLTokenizeString2(osVRTTrustedModules, ",", 0);
     544          23 :             for (char **papszIter = papszTrustedModules;
     545          23 :                  !bIsTrustedModule && papszIter && *papszIter; ++papszIter)
     546             :             {
     547          13 :                 const char *pszIterModule = *papszIter;
     548          13 :                 size_t nIterModuleLen = strlen(pszIterModule);
     549          13 :                 if (nIterModuleLen > 2 &&
     550          12 :                     strncmp(pszIterModule + nIterModuleLen - 2, ".*", 2) == 0)
     551             :                 {
     552           2 :                     bIsTrustedModule =
     553           2 :                         (strncmp(osPythonModule, pszIterModule,
     554           3 :                                  nIterModuleLen - 2) == 0) &&
     555           1 :                         (osPythonModule.size() == nIterModuleLen - 2 ||
     556           0 :                          (osPythonModule.size() >= nIterModuleLen &&
     557           0 :                           osPythonModule[nIterModuleLen - 1] == '.'));
     558             :                 }
     559          11 :                 else if (nIterModuleLen >= 1 &&
     560          11 :                          pszIterModule[nIterModuleLen - 1] == '*')
     561             :                 {
     562           4 :                     bIsTrustedModule = (strncmp(osPythonModule, pszIterModule,
     563             :                                                 nIterModuleLen - 1) == 0);
     564             :                 }
     565             :                 else
     566             :                 {
     567           7 :                     bIsTrustedModule =
     568           7 :                         (strcmp(osPythonModule, pszIterModule) == 0);
     569             :                 }
     570             :             }
     571          10 :             CSLDestroy(papszTrustedModules);
     572             :         }
     573             : 
     574          12 :         if (!bIsTrustedModule)
     575             :         {
     576           7 :             if (osPythonModule.empty())
     577             :             {
     578           2 :                 CPLError(
     579             :                     CE_Failure, CPLE_AppDefined,
     580             :                     "Python code needs to be executed, but it uses inline code "
     581             :                     "in the VRT whereas the current policy is to trust only "
     582             :                     "code from external trusted modules (defined in the "
     583             :                     "GDAL_VRT_PYTHON_TRUSTED_MODULES configuration option). "
     584             :                     "If you trust the code in %s, you can set the "
     585             :                     "GDAL_VRT_ENABLE_PYTHON configuration option to YES.",
     586           2 :                     GetDataset() ? GetDataset()->GetDescription()
     587             :                                  : "(unknown VRT)");
     588             :             }
     589           5 :             else if (osVRTTrustedModules.empty())
     590             :             {
     591           2 :                 CPLError(
     592             :                     CE_Failure, CPLE_AppDefined,
     593             :                     "Python code needs to be executed, but it uses code "
     594             :                     "from module '%s', whereas the current policy is to "
     595             :                     "trust only code from modules defined in the "
     596             :                     "GDAL_VRT_PYTHON_TRUSTED_MODULES configuration option, "
     597             :                     "which is currently unset. "
     598             :                     "If you trust the code in '%s', you can add module '%s' "
     599             :                     "to GDAL_VRT_PYTHON_TRUSTED_MODULES (or set the "
     600             :                     "GDAL_VRT_ENABLE_PYTHON configuration option to YES).",
     601             :                     osPythonModule.c_str(),
     602           1 :                     GetDataset() ? GetDataset()->GetDescription()
     603             :                                  : "(unknown VRT)",
     604             :                     osPythonModule.c_str());
     605             :             }
     606             :             else
     607             :             {
     608           8 :                 CPLError(
     609             :                     CE_Failure, CPLE_AppDefined,
     610             :                     "Python code needs to be executed, but it uses code "
     611             :                     "from module '%s', whereas the current policy is to "
     612             :                     "trust only code from modules '%s' (defined in the "
     613             :                     "GDAL_VRT_PYTHON_TRUSTED_MODULES configuration option). "
     614             :                     "If you trust the code in '%s', you can add module '%s' "
     615             :                     "to GDAL_VRT_PYTHON_TRUSTED_MODULES (or set the "
     616             :                     "GDAL_VRT_ENABLE_PYTHON configuration option to YES).",
     617             :                     osPythonModule.c_str(), osVRTTrustedModules.c_str(),
     618           4 :                     GetDataset() ? GetDataset()->GetDescription()
     619             :                                  : "(unknown VRT)",
     620             :                     osPythonModule.c_str());
     621             :             }
     622           7 :             return false;
     623             :         }
     624             :     }
     625             : 
     626             : #ifdef disabled_because_this_is_probably_broken_by_design
     627             :     // See https://lwn.net/Articles/574215/
     628             :     // and http://nedbatchelder.com/blog/201206/eval_really_is_dangerous.html
     629             :     else if (EQUAL(osPythonEnabled, "IF_SAFE"))
     630             :     {
     631             :         bool bSafe = true;
     632             :         // If the function comes from another module, then we don't know
     633             :         if (!osPythonModule.empty())
     634             :         {
     635             :             CPLDebug("VRT", "Python function is from another module");
     636             :             bSafe = false;
     637             :         }
     638             : 
     639             :         CPLString osCode(m_poPrivate->m_osCode);
     640             : 
     641             :         // Reject all imports except a few trusted modules
     642             :         const char *const apszTrustedImports[] = {
     643             :             "import math",
     644             :             "from math import",
     645             :             "import numpy",  // caution: numpy has lots of I/O functions !
     646             :             "from numpy import",
     647             :             // TODO: not sure if importing arbitrary stuff from numba is OK
     648             :             // so let's just restrict to jit.
     649             :             "from numba import jit",
     650             : 
     651             :             // Not imports but still whitelisted, whereas other __ is banned
     652             :             "__init__",
     653             :             "__call__",
     654             :         };
     655             :         for (size_t i = 0; i < CPL_ARRAYSIZE(apszTrustedImports); ++i)
     656             :         {
     657             :             osCode.replaceAll(CPLString(apszTrustedImports[i]), "");
     658             :         }
     659             : 
     660             :         // Some dangerous built-in functions or numpy functions
     661             :         const char *const apszUntrusted[] = {
     662             :             "import",  // and __import__
     663             :             "eval",       "compile", "open",
     664             :             "load",        // reload, numpy.load
     665             :             "file",        // and exec_file, numpy.fromfile, numpy.tofile
     666             :             "input",       // and raw_input
     667             :             "save",        // numpy.save
     668             :             "memmap",      // numpy.memmap
     669             :             "DataSource",  // numpy.DataSource
     670             :             "genfromtxt",  // numpy.genfromtxt
     671             :             "getattr",
     672             :             "ctypeslib",  // numpy.ctypeslib
     673             :             "testing",    // numpy.testing
     674             :             "dump",       // numpy.ndarray.dump
     675             :             "fromregex",  // numpy.fromregex
     676             :             "__"};
     677             :         for (size_t i = 0; i < CPL_ARRAYSIZE(apszUntrusted); ++i)
     678             :         {
     679             :             if (osCode.find(apszUntrusted[i]) != std::string::npos)
     680             :             {
     681             :                 CPLDebug("VRT", "Found '%s' word in Python code",
     682             :                          apszUntrusted[i]);
     683             :                 bSafe = false;
     684             :             }
     685             :         }
     686             : 
     687             :         if (!bSafe)
     688             :         {
     689             :             CPLError(CE_Failure, CPLE_AppDefined,
     690             :                      "Python code needs to be executed, but we cannot verify "
     691             :                      "if it is safe, so this is disabled by default. "
     692             :                      "If you trust the code in %s, you can set the "
     693             :                      "GDAL_VRT_ENABLE_PYTHON configuration option to YES.",
     694             :                      GetDataset() ? GetDataset()->GetDescription()
     695             :                                   : "(unknown VRT)");
     696             :             return false;
     697             :         }
     698             :     }
     699             : #endif  // disabled_because_this_is_probably_broken_by_design
     700             : 
     701          55 :     else if (!EQUAL(osPythonEnabled, "YES") && !EQUAL(osPythonEnabled, "ON") &&
     702           1 :              !EQUAL(osPythonEnabled, "TRUE"))
     703             :     {
     704           1 :         if (pszPythonEnabled == nullptr)
     705             :         {
     706             :             // Note: this is dead code with our current default policy
     707             :             // GDAL_VRT_ENABLE_PYTHON == "TRUSTED_MODULES"
     708           0 :             CPLError(CE_Failure, CPLE_AppDefined,
     709             :                      "Python code needs to be executed, but this is "
     710             :                      "disabled by default. If you trust the code in %s, "
     711             :                      "you can set the GDAL_VRT_ENABLE_PYTHON configuration "
     712             :                      "option to YES.",
     713           0 :                      GetDataset() ? GetDataset()->GetDescription()
     714             :                                   : "(unknown VRT)");
     715             :         }
     716             :         else
     717             :         {
     718           1 :             CPLError(
     719             :                 CE_Failure, CPLE_AppDefined,
     720             :                 "Python code in %s needs to be executed, but this has been "
     721             :                 "explicitly disabled.",
     722           1 :                 GetDataset() ? GetDataset()->GetDescription()
     723             :                              : "(unknown VRT)");
     724             :         }
     725           1 :         return false;
     726             :     }
     727             : 
     728          58 :     if (!GDALPythonInitialize())
     729           2 :         return false;
     730             : 
     731             :     // Whether we should just use our own global mutex, in addition to Python
     732             :     // GIL locking.
     733         112 :     m_poPrivate->m_bExclusiveLock =
     734          56 :         CPLTestBool(CPLGetConfigOption("GDAL_VRT_PYTHON_EXCLUSIVE_LOCK", "NO"));
     735             : 
     736             :     // numba jit'ification doesn't seem to be thread-safe, so force use of
     737             :     // lock now and at first execution of function. Later executions seem to
     738             :     // be thread-safe. This problem doesn't seem to appear for code in
     739             :     // regular files
     740             :     const bool bUseExclusiveLock =
     741         112 :         m_poPrivate->m_bExclusiveLock ||
     742          56 :         m_poPrivate->m_osCode.find("@jit") != std::string::npos;
     743         112 :     GIL_Holder oHolder(bUseExclusiveLock);
     744             : 
     745             :     // As we don't want to depend on numpy C API/ABI, we use a trick to build
     746             :     // a numpy array object. We define a Python function to which we pass a
     747             :     // Python buffer object.
     748             : 
     749             :     // We need to build a unique module name, otherwise this will crash in
     750             :     // multithreaded use cases.
     751         112 :     CPLString osModuleName(CPLSPrintf("gdal_vrt_module_%p", this));
     752         112 :     PyObject *poCompiledString = Py_CompileString(
     753             :         ("import numpy\n"
     754             :          "def GDALCreateNumpyArray(buffer, dtype, height, width):\n"
     755             :          "    return numpy.frombuffer(buffer, str(dtype.decode('ascii')))."
     756             :          "reshape([height, width])\n"
     757          56 :          "\n" +
     758          56 :          m_poPrivate->m_osCode)
     759             :             .c_str(),
     760             :         osModuleName, Py_file_input);
     761          56 :     if (poCompiledString == nullptr || PyErr_Occurred())
     762             :     {
     763           1 :         CPLError(CE_Failure, CPLE_AppDefined, "Couldn't compile code:\n%s",
     764           2 :                  GetPyExceptionString().c_str());
     765           1 :         return false;
     766             :     }
     767             :     PyObject *poModule =
     768          55 :         PyImport_ExecCodeModule(osModuleName, poCompiledString);
     769          55 :     Py_DecRef(poCompiledString);
     770             : 
     771          55 :     if (poModule == nullptr || PyErr_Occurred())
     772             :     {
     773           1 :         CPLError(CE_Failure, CPLE_AppDefined, "%s",
     774           2 :                  GetPyExceptionString().c_str());
     775           1 :         return false;
     776             :     }
     777             : 
     778             :     // Fetch user computation function
     779          54 :     if (!osPythonModule.empty())
     780             :     {
     781          24 :         PyObject *poUserModule = PyImport_ImportModule(osPythonModule);
     782          24 :         if (poUserModule == nullptr || PyErr_Occurred())
     783             :         {
     784           1 :             CPLString osException = GetPyExceptionString();
     785           1 :             if (!osException.empty() && osException.back() == '\n')
     786             :             {
     787           1 :                 osException.pop_back();
     788             :             }
     789           1 :             if (osException.find("ModuleNotFoundError") == 0)
     790             :             {
     791           1 :                 osException += ". You may need to define PYTHONPATH";
     792             :             }
     793           1 :             CPLError(CE_Failure, CPLE_AppDefined, "%s", osException.c_str());
     794           1 :             Py_DecRef(poModule);
     795           1 :             return false;
     796             :         }
     797          46 :         m_poPrivate->m_poUserFunction =
     798          23 :             PyObject_GetAttrString(poUserModule, osPythonFunction);
     799          23 :         Py_DecRef(poUserModule);
     800             :     }
     801             :     else
     802             :     {
     803          60 :         m_poPrivate->m_poUserFunction =
     804          30 :             PyObject_GetAttrString(poModule, osPythonFunction);
     805             :     }
     806          53 :     if (m_poPrivate->m_poUserFunction == nullptr || PyErr_Occurred())
     807             :     {
     808           1 :         CPLError(CE_Failure, CPLE_AppDefined, "%s",
     809           2 :                  GetPyExceptionString().c_str());
     810           1 :         Py_DecRef(poModule);
     811           1 :         return false;
     812             :     }
     813          52 :     if (!PyCallable_Check(m_poPrivate->m_poUserFunction))
     814             :     {
     815           1 :         CPLError(CE_Failure, CPLE_AppDefined, "Object '%s' is not callable",
     816             :                  osPythonFunction.c_str());
     817           1 :         Py_DecRef(poModule);
     818           1 :         return false;
     819             :     }
     820             : 
     821             :     // Fetch our GDALCreateNumpyArray python function
     822         102 :     m_poPrivate->m_poGDALCreateNumpyArray =
     823          51 :         PyObject_GetAttrString(poModule, "GDALCreateNumpyArray");
     824          51 :     if (m_poPrivate->m_poGDALCreateNumpyArray == nullptr || PyErr_Occurred())
     825             :     {
     826             :         // Shouldn't happen normally...
     827           0 :         CPLError(CE_Failure, CPLE_AppDefined, "%s",
     828           0 :                  GetPyExceptionString().c_str());
     829           0 :         Py_DecRef(poModule);
     830           0 :         return false;
     831             :     }
     832          51 :     Py_DecRef(poModule);
     833             : 
     834          51 :     m_poPrivate->m_bPythonInitializationSuccess = true;
     835          51 :     return true;
     836             : }
     837             : 
     838        2776 : CPLErr VRTDerivedRasterBand::GetPixelFunctionArguments(
     839             :     const CPLString &osMetadata,
     840             :     const std::vector<int> &anMapBufferIdxToSourceIdx, int nXOff, int nYOff,
     841             :     std::vector<std::pair<CPLString, CPLString>> &oAdditionalArgs)
     842             : {
     843             : 
     844        5552 :     auto poArgs = CPLXMLTreeCloser(CPLParseXMLString(osMetadata));
     845        5552 :     if (poArgs != nullptr && poArgs->eType == CXT_Element &&
     846        2776 :         !strcmp(poArgs->pszValue, "PixelFunctionArgumentsList"))
     847             :     {
     848       12701 :         for (CPLXMLNode *psIter = poArgs->psChild; psIter != nullptr;
     849        9925 :              psIter = psIter->psNext)
     850             :         {
     851        9927 :             if (psIter->eType == CXT_Element &&
     852        9927 :                 !strcmp(psIter->pszValue, "Argument"))
     853             :             {
     854        9927 :                 CPLString osName, osType, osValue;
     855        9927 :                 auto pszName = CPLGetXMLValue(psIter, "name", nullptr);
     856        9927 :                 if (pszName != nullptr)
     857        5647 :                     osName = pszName;
     858        9927 :                 auto pszType = CPLGetXMLValue(psIter, "type", nullptr);
     859        9927 :                 if (pszType != nullptr)
     860        9927 :                     osType = pszType;
     861        9927 :                 auto pszValue = CPLGetXMLValue(psIter, "value", nullptr);
     862        9927 :                 if (pszValue != nullptr)
     863        4658 :                     osValue = pszValue;
     864        9927 :                 if (osType == "constant" && osValue != "" && osName != "")
     865           1 :                     oAdditionalArgs.push_back(
     866           2 :                         std::pair<CPLString, CPLString>(osName, osValue));
     867        9927 :                 if (osType == "builtin")
     868             :                 {
     869        4280 :                     const CPLString &osArgName = osValue;
     870        4280 :                     CPLString osVal;
     871        4280 :                     double dfVal = 0;
     872             : 
     873        4280 :                     int success(FALSE);
     874        4280 :                     if (osArgName == "NoData")
     875        2766 :                         dfVal = this->GetNoDataValue(&success);
     876        1514 :                     else if (osArgName == "scale")
     877           3 :                         dfVal = this->GetScale(&success);
     878        1511 :                     else if (osArgName == "offset")
     879           2 :                         dfVal = this->GetOffset(&success);
     880        1509 :                     else if (osArgName == "xoff")
     881             :                     {
     882         377 :                         dfVal = static_cast<double>(nXOff);
     883         377 :                         success = true;
     884             :                     }
     885        1132 :                     else if (osArgName == "yoff")
     886             :                     {
     887         377 :                         dfVal = static_cast<double>(nYOff);
     888         377 :                         success = true;
     889             :                     }
     890         755 :                     else if (osArgName == "crs")
     891             :                     {
     892             :                         const auto *crs =
     893           5 :                             GetDataset()->GetSpatialRefRasterOnly();
     894           5 :                         if (crs)
     895             :                         {
     896             :                             osVal =
     897           4 :                                 std::to_string(reinterpret_cast<size_t>(crs));
     898           4 :                             success = true;
     899             :                         }
     900             :                         else
     901             :                         {
     902           1 :                             CPLError(CE_Failure, CPLE_AppDefined,
     903             :                                      "VRTDataset has no <SRS>");
     904             :                         }
     905             :                     }
     906         750 :                     else if (osArgName == "geotransform")
     907             :                     {
     908         377 :                         GDALGeoTransform gt;
     909         377 :                         if (GetDataset()->GetGeoTransform(gt) != CE_None)
     910             :                         {
     911             :                             // Do not fail here because the argument is most
     912             :                             // likely not needed by the pixel function. If it
     913             :                             // is needed, the pixel function can emit the error.
     914         143 :                             continue;
     915             :                         }
     916             :                         osVal = CPLSPrintf(
     917         468 :                             "%.17g,%.17g,%.17g,%.17g,%.17g,%.17g", gt[0], gt[1],
     918         234 :                             gt[2], gt[3], gt[4], gt[5]);
     919         234 :                         success = true;
     920             :                     }
     921         373 :                     else if (osArgName == "source_names")
     922             :                     {
     923         984 :                         for (size_t iBuffer = 0;
     924         984 :                              iBuffer < anMapBufferIdxToSourceIdx.size();
     925             :                              iBuffer++)
     926             :                         {
     927             :                             const int iSource =
     928         611 :                                 anMapBufferIdxToSourceIdx[iBuffer];
     929             :                             const VRTSource *poSource =
     930         611 :                                 m_papoSources[iSource].get();
     931             : 
     932         611 :                             if (iBuffer > 0)
     933             :                             {
     934         245 :                                 osVal += "|";
     935             :                             }
     936             : 
     937         611 :                             const auto &osSourceName = poSource->GetName();
     938         611 :                             if (osSourceName.empty())
     939             :                             {
     940          42 :                                 osVal += "B" + std::to_string(iBuffer + 1);
     941             :                             }
     942             :                             else
     943             :                             {
     944         569 :                                 osVal += osSourceName;
     945             :                             }
     946             :                         }
     947             : 
     948         373 :                         success = true;
     949             :                     }
     950             :                     else
     951             :                     {
     952           0 :                         CPLError(
     953             :                             CE_Failure, CPLE_NotSupported,
     954             :                             "PixelFunction builtin argument %s not supported",
     955             :                             osArgName.c_str());
     956           0 :                         return CE_Failure;
     957             :                     }
     958        4137 :                     if (!success)
     959             :                     {
     960        2645 :                         if (CPLTestBool(
     961             :                                 CPLGetXMLValue(psIter, "optional", "false")))
     962        2643 :                             continue;
     963             : 
     964           2 :                         CPLError(CE_Failure, CPLE_AppDefined,
     965             :                                  "Raster has no %s", osValue.c_str());
     966           2 :                         return CE_Failure;
     967             :                     }
     968             : 
     969        1492 :                     if (osVal.empty())
     970             :                     {
     971         888 :                         osVal = CPLSPrintf("%.17g", dfVal);
     972             :                     }
     973             : 
     974        1492 :                     oAdditionalArgs.push_back(
     975        2984 :                         std::pair<CPLString, CPLString>(osArgName, osVal));
     976        1492 :                     CPLDebug("VRT",
     977             :                              "Added builtin pixel function argument %s = %s",
     978             :                              osArgName.c_str(), osVal.c_str());
     979             :                 }
     980             :             }
     981             :         }
     982             :     }
     983             : 
     984        2774 :     return CE_None;
     985             : }
     986             : 
     987             : /************************************************************************/
     988             : /*                             IRasterIO()                              */
     989             : /************************************************************************/
     990             : 
     991             : /**
     992             :  * Read/write a region of image data for this band.
     993             :  *
     994             :  * Each of the sources for this derived band will be read and passed to
     995             :  * the derived band pixel function.  The pixel function is responsible
     996             :  * for applying whatever algorithm is necessary to generate this band's
     997             :  * pixels from the sources.
     998             :  *
     999             :  * The sources will be read using the transfer type specified for sources
    1000             :  * using SetSourceTransferType().  If no transfer type has been set for
    1001             :  * this derived band, the band's data type will be used as the transfer type.
    1002             :  *
    1003             :  * @see gdalrasterband
    1004             :  *
    1005             :  * @param eRWFlag Either GF_Read to read a region of data, or GT_Write to
    1006             :  * write a region of data.
    1007             :  *
    1008             :  * @param nXOff The pixel offset to the top left corner of the region
    1009             :  * of the band to be accessed.  This would be zero to start from the left side.
    1010             :  *
    1011             :  * @param nYOff The line offset to the top left corner of the region
    1012             :  * of the band to be accessed.  This would be zero to start from the top.
    1013             :  *
    1014             :  * @param nXSize The width of the region of the band to be accessed in pixels.
    1015             :  *
    1016             :  * @param nYSize The height of the region of the band to be accessed in lines.
    1017             :  *
    1018             :  * @param pData The buffer into which the data should be read, or from which
    1019             :  * it should be written.  This buffer must contain at least nBufXSize *
    1020             :  * nBufYSize words of type eBufType.  It is organized in left to right,
    1021             :  * top to bottom pixel order.  Spacing is controlled by the nPixelSpace,
    1022             :  * and nLineSpace parameters.
    1023             :  *
    1024             :  * @param nBufXSize The width of the buffer image into which the desired
    1025             :  * region is to be read, or from which it is to be written.
    1026             :  *
    1027             :  * @param nBufYSize The height of the buffer image into which the desired
    1028             :  * region is to be read, or from which it is to be written.
    1029             :  *
    1030             :  * @param eBufType The type of the pixel values in the pData data buffer.  The
    1031             :  * pixel values will automatically be translated to/from the GDALRasterBand
    1032             :  * data type as needed.
    1033             :  *
    1034             :  * @param nPixelSpace The byte offset from the start of one pixel value in
    1035             :  * pData to the start of the next pixel value within a scanline.  If defaulted
    1036             :  * (0) the size of the datatype eBufType is used.
    1037             :  *
    1038             :  * @param nLineSpace The byte offset from the start of one scanline in
    1039             :  * pData to the start of the next.  If defaulted the size of the datatype
    1040             :  * eBufType * nBufXSize is used.
    1041             :  *
    1042             :  * @return CE_Failure if the access fails, otherwise CE_None.
    1043             :  */
    1044        4305 : CPLErr VRTDerivedRasterBand::IRasterIO(
    1045             :     GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
    1046             :     void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
    1047             :     GSpacing nPixelSpace, GSpacing nLineSpace, GDALRasterIOExtraArg *psExtraArg)
    1048             : {
    1049        4305 :     if (eRWFlag == GF_Write)
    1050             :     {
    1051           1 :         CPLError(CE_Failure, CPLE_AppDefined,
    1052             :                  "Writing through VRTSourcedRasterBand is not supported.");
    1053           1 :         return CE_Failure;
    1054             :     }
    1055             : 
    1056        8608 :     const std::string osFctId("VRTDerivedRasterBand::IRasterIO");
    1057        8608 :     GDALAntiRecursionGuard oGuard(osFctId);
    1058        4304 :     if (oGuard.GetCallDepth() >= 32)
    1059             :     {
    1060           0 :         CPLError(
    1061             :             CE_Failure, CPLE_AppDefined,
    1062             :             "VRTDerivedRasterBand::IRasterIO(): Recursion detected (case 1)");
    1063           0 :         return CE_Failure;
    1064             :     }
    1065             : 
    1066       12912 :     GDALAntiRecursionGuard oGuard2(oGuard, poDS->GetDescription());
    1067             :     // Allow multiple recursion depths on the same dataset in case the split strategy is applied
    1068        4304 :     if (oGuard2.GetCallDepth() > 15)
    1069             :     {
    1070           0 :         CPLError(
    1071             :             CE_Failure, CPLE_AppDefined,
    1072             :             "VRTDerivedRasterBand::IRasterIO(): Recursion detected (case 2)");
    1073           0 :         return CE_Failure;
    1074             :     }
    1075             : 
    1076             :     if constexpr (sizeof(GSpacing) > sizeof(int))
    1077             :     {
    1078        4304 :         if (nLineSpace > INT_MAX)
    1079             :         {
    1080           0 :             if (nBufYSize == 1)
    1081             :             {
    1082           0 :                 nLineSpace = 0;
    1083             :             }
    1084             :             else
    1085             :             {
    1086           0 :                 CPLError(CE_Failure, CPLE_NotSupported,
    1087             :                          "VRTDerivedRasterBand::IRasterIO(): nLineSpace > "
    1088             :                          "INT_MAX not supported");
    1089           0 :                 return CE_Failure;
    1090             :             }
    1091             :         }
    1092             :     }
    1093             : 
    1094             :     /* -------------------------------------------------------------------- */
    1095             :     /*      Do we have overviews that would be appropriate to satisfy       */
    1096             :     /*      this request?                                                   */
    1097             :     /* -------------------------------------------------------------------- */
    1098        4304 :     auto l_poDS = dynamic_cast<VRTDataset *>(poDS);
    1099        4304 :     if (l_poDS &&
    1100        8607 :         l_poDS->m_apoOverviews.empty() &&  // do not use virtual overviews
    1101        8608 :         (nBufXSize < nXSize || nBufYSize < nYSize) && GetOverviewCount() > 0)
    1102             :     {
    1103           0 :         if (OverviewRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize, pData,
    1104             :                              nBufXSize, nBufYSize, eBufType, nPixelSpace,
    1105           0 :                              nLineSpace, psExtraArg) == CE_None)
    1106           0 :             return CE_None;
    1107             :     }
    1108             : 
    1109        4304 :     const int nBufTypeSize = GDALGetDataTypeSizeBytes(eBufType);
    1110        4304 :     GDALDataType eSrcType = eSourceTransferType;
    1111        4304 :     if (eSrcType == GDT_Unknown || eSrcType >= GDT_TypeCount)
    1112             :     {
    1113             :         // Check the largest data type for all sources
    1114        3327 :         GDALDataType eAllSrcType = GDT_Unknown;
    1115      107461 :         for (auto &poSource : m_papoSources)
    1116             :         {
    1117      104135 :             if (poSource->IsSimpleSource())
    1118             :             {
    1119             :                 const auto poSS =
    1120      104134 :                     static_cast<VRTSimpleSource *>(poSource.get());
    1121      104134 :                 auto l_poBand = poSS->GetRasterBand();
    1122      104134 :                 if (l_poBand)
    1123             :                 {
    1124      104134 :                     eAllSrcType = GDALDataTypeUnion(
    1125             :                         eAllSrcType, l_poBand->GetRasterDataType());
    1126             :                 }
    1127             :                 else
    1128             :                 {
    1129           0 :                     eAllSrcType = GDT_Unknown;
    1130           0 :                     break;
    1131             :                 }
    1132             :             }
    1133             :             else
    1134             :             {
    1135           1 :                 eAllSrcType = GDT_Unknown;
    1136           1 :                 break;
    1137             :             }
    1138             :         }
    1139             : 
    1140        3327 :         if (eAllSrcType != GDT_Unknown)
    1141        2941 :             eSrcType = GDALDataTypeUnion(eAllSrcType, eDataType);
    1142             :         else
    1143         386 :             eSrcType = GDALDataTypeUnion(GDT_Float64, eDataType);
    1144             :     }
    1145        4304 :     const int nSrcTypeSize = GDALGetDataTypeSizeBytes(eSrcType);
    1146             : 
    1147             :     // If acquiring the region of interest in a single time is going
    1148             :     // to consume too much RAM, split in halves, and that recursively
    1149             :     // until we get below m_nAllowedRAMUsage.
    1150        4304 :     if (m_poPrivate->m_nAllowedRAMUsage > 0 && !m_papoSources.empty() &&
    1151       12521 :         nSrcTypeSize > 0 && nBufXSize == nXSize && nBufYSize == nYSize &&
    1152        3913 :         static_cast<GIntBig>(nBufXSize) * nBufYSize >
    1153        7826 :             m_poPrivate->m_nAllowedRAMUsage /
    1154        3913 :                 (static_cast<int>(m_papoSources.size()) * nSrcTypeSize))
    1155             :     {
    1156         999 :         CPLErr eErr = SplitRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
    1157             :                                     pData, nBufXSize, nBufYSize, eBufType,
    1158             :                                     nPixelSpace, nLineSpace, psExtraArg);
    1159         999 :         if (eErr != CE_Warning)
    1160         999 :             return eErr;
    1161             :     }
    1162             : 
    1163             :     /* ---- Get pixel function for band ---- */
    1164        3305 :     const std::pair<PixelFunc, std::string> *poPixelFunc = nullptr;
    1165        6610 :     std::vector<std::pair<CPLString, CPLString>> oAdditionalArgs;
    1166             : 
    1167        3305 :     if (EQUAL(m_poPrivate->m_osLanguage, "C"))
    1168             :     {
    1169             :         poPixelFunc =
    1170        2832 :             VRTDerivedRasterBand::GetPixelFunction(osFuncName.c_str());
    1171        2832 :         if (poPixelFunc == nullptr)
    1172             :         {
    1173           1 :             CPLError(CE_Failure, CPLE_IllegalArg,
    1174             :                      "VRTDerivedRasterBand::IRasterIO:"
    1175             :                      "Derived band pixel function '%s' not registered.",
    1176             :                      osFuncName.c_str());
    1177           1 :             return CE_Failure;
    1178             :         }
    1179             :     }
    1180             : 
    1181             :     /* TODO: It would be nice to use a MallocBlock function for each
    1182             :        individual buffer that would recycle blocks of memory from a
    1183             :        cache by reassigning blocks that are nearly the same size.
    1184             :        A corresponding FreeBlock might only truly free if the total size
    1185             :        of freed blocks gets to be too great of a percentage of the size
    1186             :        of the allocated blocks. */
    1187             : 
    1188             :     // Get buffers for each source.
    1189        3304 :     const int nBufferRadius = m_poPrivate->m_nBufferRadius;
    1190        3304 :     if (nBufferRadius > (INT_MAX - nBufXSize) / 2 ||
    1191        3304 :         nBufferRadius > (INT_MAX - nBufYSize) / 2)
    1192             :     {
    1193           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    1194             :                  "Integer overflow: "
    1195             :                  "nBufferRadius > (INT_MAX - nBufXSize) / 2 || "
    1196             :                  "nBufferRadius > (INT_MAX - nBufYSize) / 2)");
    1197           0 :         return CE_Failure;
    1198             :     }
    1199        3304 :     const int nExtBufXSize = nBufXSize + 2 * nBufferRadius;
    1200        3304 :     const int nExtBufYSize = nBufYSize + 2 * nBufferRadius;
    1201        3304 :     int nBufferCount = 0;
    1202             : 
    1203             :     std::vector<std::unique_ptr<void, VSIFreeReleaser>> apBuffers(
    1204        6608 :         m_papoSources.size());
    1205        6608 :     std::vector<int> anMapBufferIdxToSourceIdx(m_papoSources.size());
    1206        3304 :     bool bSkipOutputBufferInitialization = !m_papoSources.empty();
    1207      107247 :     for (int iSource = 0; iSource < static_cast<int>(m_papoSources.size());
    1208             :          iSource++)
    1209             :     {
    1210      103957 :         if (m_poPrivate->m_bSkipNonContributingSources &&
    1211          14 :             m_papoSources[iSource]->IsSimpleSource())
    1212             :         {
    1213          14 :             bool bError = false;
    1214             :             double dfReqXOff, dfReqYOff, dfReqXSize, dfReqYSize;
    1215             :             int nReqXOff, nReqYOff, nReqXSize, nReqYSize;
    1216             :             int nOutXOff, nOutYOff, nOutXSize, nOutYSize;
    1217             :             auto poSource =
    1218          14 :                 static_cast<VRTSimpleSource *>(m_papoSources[iSource].get());
    1219          14 :             if (!poSource->GetSrcDstWindow(
    1220             :                     nXOff, nYOff, nXSize, nYSize, nBufXSize, nBufYSize,
    1221             :                     psExtraArg->eResampleAlg, &dfReqXOff, &dfReqYOff,
    1222             :                     &dfReqXSize, &dfReqYSize, &nReqXOff, &nReqYOff, &nReqXSize,
    1223             :                     &nReqYSize, &nOutXOff, &nOutYOff, &nOutXSize, &nOutYSize,
    1224             :                     bError))
    1225             :             {
    1226           4 :                 if (bError)
    1227             :                 {
    1228           0 :                     return CE_Failure;
    1229             :                 }
    1230             : 
    1231             :                 // Skip non contributing source
    1232           4 :                 bSkipOutputBufferInitialization = false;
    1233           4 :                 continue;
    1234             :             }
    1235             :         }
    1236             : 
    1237      103939 :         anMapBufferIdxToSourceIdx[nBufferCount] = iSource;
    1238      103939 :         apBuffers[nBufferCount].reset(
    1239             :             VSI_MALLOC3_VERBOSE(nSrcTypeSize, nExtBufXSize, nExtBufYSize));
    1240      103939 :         if (apBuffers[nBufferCount] == nullptr)
    1241             :         {
    1242           0 :             return CE_Failure;
    1243             :         }
    1244             : 
    1245      103939 :         bool bBufferInit = true;
    1246      103939 :         if (m_papoSources[iSource]->IsSimpleSource())
    1247             :         {
    1248             :             const auto poSS =
    1249      103938 :                 static_cast<VRTSimpleSource *>(m_papoSources[iSource].get());
    1250      103938 :             auto l_poBand = poSS->GetRasterBand();
    1251      103938 :             if (l_poBand != nullptr && poSS->m_dfSrcXOff == 0.0 &&
    1252        1046 :                 poSS->m_dfSrcYOff == 0.0 &&
    1253        2092 :                 poSS->m_dfSrcXOff + poSS->m_dfSrcXSize ==
    1254        1046 :                     l_poBand->GetXSize() &&
    1255        2056 :                 poSS->m_dfSrcYOff + poSS->m_dfSrcYSize ==
    1256        1028 :                     l_poBand->GetYSize() &&
    1257        1028 :                 poSS->m_dfDstXOff == 0.0 && poSS->m_dfDstYOff == 0.0 &&
    1258      208898 :                 poSS->m_dfDstXOff + poSS->m_dfDstXSize == nRasterXSize &&
    1259        1022 :                 poSS->m_dfDstYOff + poSS->m_dfDstYSize == nRasterYSize)
    1260             :             {
    1261        1022 :                 if (m_papoSources[iSource]->GetType() ==
    1262        1022 :                     VRTSimpleSource::GetTypeStatic())
    1263         995 :                     bBufferInit = false;
    1264             :             }
    1265             :             else
    1266             :             {
    1267      102916 :                 bSkipOutputBufferInitialization = false;
    1268             :             }
    1269             :         }
    1270             :         else
    1271             :         {
    1272           1 :             bSkipOutputBufferInitialization = false;
    1273             :         }
    1274      103939 :         if (bBufferInit)
    1275             :         {
    1276             :             /* ------------------------------------------------------------ */
    1277             :             /* #4045: Initialize the newly allocated buffers before handing */
    1278             :             /* them off to the sources. These buffers are packed, so we     */
    1279             :             /* don't need any special line-by-line handling when a nonzero  */
    1280             :             /* nodata value is set.                                         */
    1281             :             /* ------------------------------------------------------------ */
    1282      102944 :             if (!m_bNoDataValueSet || m_dfNoDataValue == 0)
    1283             :             {
    1284      102671 :                 memset(apBuffers[nBufferCount].get(), 0,
    1285      102671 :                        static_cast<size_t>(nSrcTypeSize) * nExtBufXSize *
    1286      102671 :                            nExtBufYSize);
    1287             :             }
    1288             :             else
    1289             :             {
    1290         546 :                 GDALCopyWords64(
    1291         273 :                     &m_dfNoDataValue, GDT_Float64, 0,
    1292         273 :                     static_cast<GByte *>(apBuffers[nBufferCount].get()),
    1293             :                     eSrcType, nSrcTypeSize,
    1294         273 :                     static_cast<GPtrDiff_t>(nExtBufXSize) * nExtBufYSize);
    1295             :             }
    1296             :         }
    1297             : 
    1298      103939 :         ++nBufferCount;
    1299             :     }
    1300             : 
    1301             :     /* -------------------------------------------------------------------- */
    1302             :     /*      Initialize the buffer to some background value. Use the         */
    1303             :     /*      nodata value if available.                                      */
    1304             :     /* -------------------------------------------------------------------- */
    1305        3304 :     if (bSkipOutputBufferInitialization)
    1306             :     {
    1307             :         // Do nothing
    1308             :     }
    1309        2602 :     else if (nPixelSpace == nBufTypeSize &&
    1310        2584 :              (!m_bNoDataValueSet || m_dfNoDataValue == 0))
    1311             :     {
    1312        2486 :         memset(pData, 0,
    1313        2486 :                static_cast<size_t>(nBufXSize) * nBufYSize * nBufTypeSize);
    1314             :     }
    1315         116 :     else if (m_bNoDataValueSet)
    1316             :     {
    1317          98 :         double dfWriteValue = m_dfNoDataValue;
    1318             : 
    1319         261 :         for (int iLine = 0; iLine < nBufYSize; iLine++)
    1320             :         {
    1321         163 :             GDALCopyWords64(&dfWriteValue, GDT_Float64, 0,
    1322         163 :                             static_cast<GByte *>(pData) + nLineSpace * iLine,
    1323             :                             eBufType, static_cast<int>(nPixelSpace), nBufXSize);
    1324             :         }
    1325             :     }
    1326             : 
    1327             :     // No contributing sources and SkipNonContributingSources mode ?
    1328             :     // Do not call the pixel function and just return the 0/nodata initialized
    1329             :     // output buffer.
    1330        3304 :     if (nBufferCount == 0 && m_poPrivate->m_bSkipNonContributingSources)
    1331             :     {
    1332           1 :         return CE_None;
    1333             :     }
    1334             : 
    1335             :     GDALRasterIOExtraArg sExtraArg;
    1336        3303 :     GDALCopyRasterIOExtraArg(&sExtraArg, psExtraArg);
    1337             : 
    1338        3303 :     int nXShiftInBuffer = 0;
    1339        3303 :     int nYShiftInBuffer = 0;
    1340        3303 :     int nExtBufXSizeReq = nExtBufXSize;
    1341        3303 :     int nExtBufYSizeReq = nExtBufYSize;
    1342             : 
    1343        3303 :     int nXOffExt = nXOff;
    1344        3303 :     int nYOffExt = nYOff;
    1345        3303 :     int nXSizeExt = nXSize;
    1346        3303 :     int nYSizeExt = nYSize;
    1347             : 
    1348        3303 :     if (nBufferRadius)
    1349             :     {
    1350          88 :         double dfXRatio = static_cast<double>(nXSize) / nBufXSize;
    1351          88 :         double dfYRatio = static_cast<double>(nYSize) / nBufYSize;
    1352             : 
    1353          88 :         if (!sExtraArg.bFloatingPointWindowValidity)
    1354             :         {
    1355          88 :             sExtraArg.dfXOff = nXOff;
    1356          88 :             sExtraArg.dfYOff = nYOff;
    1357          88 :             sExtraArg.dfXSize = nXSize;
    1358          88 :             sExtraArg.dfYSize = nYSize;
    1359             :         }
    1360             : 
    1361          88 :         sExtraArg.dfXOff -= dfXRatio * nBufferRadius;
    1362          88 :         sExtraArg.dfYOff -= dfYRatio * nBufferRadius;
    1363          88 :         sExtraArg.dfXSize += 2 * dfXRatio * nBufferRadius;
    1364          88 :         sExtraArg.dfYSize += 2 * dfYRatio * nBufferRadius;
    1365          88 :         if (sExtraArg.dfXOff < 0)
    1366             :         {
    1367          88 :             nXShiftInBuffer = -static_cast<int>(sExtraArg.dfXOff / dfXRatio);
    1368          88 :             nExtBufXSizeReq -= nXShiftInBuffer;
    1369          88 :             sExtraArg.dfXSize += sExtraArg.dfXOff;
    1370          88 :             sExtraArg.dfXOff = 0;
    1371             :         }
    1372          88 :         if (sExtraArg.dfYOff < 0)
    1373             :         {
    1374          88 :             nYShiftInBuffer = -static_cast<int>(sExtraArg.dfYOff / dfYRatio);
    1375          88 :             nExtBufYSizeReq -= nYShiftInBuffer;
    1376          88 :             sExtraArg.dfYSize += sExtraArg.dfYOff;
    1377          88 :             sExtraArg.dfYOff = 0;
    1378             :         }
    1379          88 :         if (sExtraArg.dfXOff + sExtraArg.dfXSize > nRasterXSize)
    1380             :         {
    1381          88 :             nExtBufXSizeReq -= static_cast<int>(
    1382          88 :                 (sExtraArg.dfXOff + sExtraArg.dfXSize - nRasterXSize) /
    1383             :                 dfXRatio);
    1384          88 :             sExtraArg.dfXSize = nRasterXSize - sExtraArg.dfXOff;
    1385             :         }
    1386          88 :         if (sExtraArg.dfYOff + sExtraArg.dfYSize > nRasterYSize)
    1387             :         {
    1388          88 :             nExtBufYSizeReq -= static_cast<int>(
    1389          88 :                 (sExtraArg.dfYOff + sExtraArg.dfYSize - nRasterYSize) /
    1390             :                 dfYRatio);
    1391          88 :             sExtraArg.dfYSize = nRasterYSize - sExtraArg.dfYOff;
    1392             :         }
    1393             : 
    1394          88 :         nXOffExt = static_cast<int>(sExtraArg.dfXOff);
    1395          88 :         nYOffExt = static_cast<int>(sExtraArg.dfYOff);
    1396         176 :         nXSizeExt = std::min(static_cast<int>(sExtraArg.dfXSize + 0.5),
    1397          88 :                              nRasterXSize - nXOffExt);
    1398         176 :         nYSizeExt = std::min(static_cast<int>(sExtraArg.dfYSize + 0.5),
    1399          88 :                              nRasterYSize - nYOffExt);
    1400             :     }
    1401             : 
    1402             :     // Load values for sources into packed buffers.
    1403        3303 :     CPLErr eErr = CE_None;
    1404        6606 :     VRTSource::WorkingState oWorkingState;
    1405      107242 :     for (int iBuffer = 0; iBuffer < nBufferCount && eErr == CE_None; iBuffer++)
    1406             :     {
    1407      103939 :         const int iSource = anMapBufferIdxToSourceIdx[iBuffer];
    1408      103939 :         GByte *pabyBuffer = static_cast<GByte *>(apBuffers[iBuffer].get());
    1409      103939 :         eErr = static_cast<VRTSource *>(m_papoSources[iSource].get())
    1410      207878 :                    ->RasterIO(
    1411             :                        eSrcType, nXOffExt, nYOffExt, nXSizeExt, nYSizeExt,
    1412      103939 :                        pabyBuffer + (static_cast<size_t>(nYShiftInBuffer) *
    1413      103939 :                                          nExtBufXSize +
    1414      103939 :                                      nXShiftInBuffer) *
    1415      103939 :                                         nSrcTypeSize,
    1416             :                        nExtBufXSizeReq, nExtBufYSizeReq, eSrcType, nSrcTypeSize,
    1417      103939 :                        static_cast<GSpacing>(nSrcTypeSize) * nExtBufXSize,
    1418      103939 :                        &sExtraArg, oWorkingState);
    1419             : 
    1420             :         // Extend first lines
    1421      104027 :         for (int iY = 0; iY < nYShiftInBuffer; iY++)
    1422             :         {
    1423          88 :             memcpy(pabyBuffer +
    1424          88 :                        static_cast<size_t>(iY) * nExtBufXSize * nSrcTypeSize,
    1425          88 :                    pabyBuffer + static_cast<size_t>(nYShiftInBuffer) *
    1426          88 :                                     nExtBufXSize * nSrcTypeSize,
    1427          88 :                    static_cast<size_t>(nExtBufXSize) * nSrcTypeSize);
    1428             :         }
    1429             :         // Extend last lines
    1430      104027 :         for (int iY = nYShiftInBuffer + nExtBufYSizeReq; iY < nExtBufYSize;
    1431             :              iY++)
    1432             :         {
    1433          88 :             memcpy(pabyBuffer +
    1434          88 :                        static_cast<size_t>(iY) * nExtBufXSize * nSrcTypeSize,
    1435          88 :                    pabyBuffer + static_cast<size_t>(nYShiftInBuffer +
    1436          88 :                                                     nExtBufYSizeReq - 1) *
    1437          88 :                                     nExtBufXSize * nSrcTypeSize,
    1438          88 :                    static_cast<size_t>(nExtBufXSize) * nSrcTypeSize);
    1439             :         }
    1440             :         // Extend first cols
    1441      103939 :         if (nXShiftInBuffer)
    1442             :         {
    1443       10912 :             for (int iY = 0; iY < nExtBufYSize; iY++)
    1444             :             {
    1445       21648 :                 for (int iX = 0; iX < nXShiftInBuffer; iX++)
    1446             :                 {
    1447       10824 :                     memcpy(pabyBuffer +
    1448       10824 :                                static_cast<size_t>(iY * nExtBufXSize + iX) *
    1449       10824 :                                    nSrcTypeSize,
    1450       10824 :                            pabyBuffer +
    1451       10824 :                                (static_cast<size_t>(iY) * nExtBufXSize +
    1452       10824 :                                 nXShiftInBuffer) *
    1453       10824 :                                    nSrcTypeSize,
    1454             :                            nSrcTypeSize);
    1455             :                 }
    1456             :             }
    1457             :         }
    1458             :         // Extent last cols
    1459      103939 :         if (nXShiftInBuffer + nExtBufXSizeReq < nExtBufXSize)
    1460             :         {
    1461       10912 :             for (int iY = 0; iY < nExtBufYSize; iY++)
    1462             :             {
    1463       10824 :                 for (int iX = nXShiftInBuffer + nExtBufXSizeReq;
    1464       21648 :                      iX < nExtBufXSize; iX++)
    1465             :                 {
    1466       10824 :                     memcpy(pabyBuffer +
    1467       10824 :                                (static_cast<size_t>(iY) * nExtBufXSize + iX) *
    1468       10824 :                                    nSrcTypeSize,
    1469       10824 :                            pabyBuffer +
    1470       10824 :                                (static_cast<size_t>(iY) * nExtBufXSize +
    1471       10824 :                                 nXShiftInBuffer + nExtBufXSizeReq - 1) *
    1472       10824 :                                    nSrcTypeSize,
    1473             :                            nSrcTypeSize);
    1474             :                 }
    1475             :             }
    1476             :         }
    1477             :     }
    1478             : 
    1479             :     // Collect any pixel function arguments into oAdditionalArgs
    1480        3303 :     if (poPixelFunc != nullptr && !poPixelFunc->second.empty())
    1481             :     {
    1482        5552 :         if (GetPixelFunctionArguments(poPixelFunc->second,
    1483             :                                       anMapBufferIdxToSourceIdx, nXOff, nYOff,
    1484        2776 :                                       oAdditionalArgs) != CE_None)
    1485             :         {
    1486           2 :             eErr = CE_Failure;
    1487             :         }
    1488             :     }
    1489             : 
    1490             :     // Apply pixel function.
    1491        3303 :     if (eErr == CE_None && EQUAL(m_poPrivate->m_osLanguage, "Python"))
    1492             :     {
    1493             :         // numpy doesn't have native cint16/cint32/cfloat16
    1494         472 :         if (eSrcType == GDT_CInt16 || eSrcType == GDT_CInt32 ||
    1495             :             eSrcType == GDT_CFloat16)
    1496             :         {
    1497           2 :             CPLError(CE_Failure, CPLE_AppDefined,
    1498             :                      "CInt16/CInt32/CFloat16 data type not supported for "
    1499             :                      "SourceTransferType");
    1500          24 :             return CE_Failure;
    1501             :         }
    1502         470 :         if (eDataType == GDT_CInt16 || eDataType == GDT_CInt32 ||
    1503         466 :             eDataType == GDT_CFloat16)
    1504             :         {
    1505           7 :             CPLError(
    1506             :                 CE_Failure, CPLE_AppDefined,
    1507             :                 "CInt16/CInt32/CFloat16 data type not supported for data type");
    1508           7 :             return CE_Failure;
    1509             :         }
    1510             : 
    1511         463 :         if (!InitializePython())
    1512          15 :             return CE_Failure;
    1513             : 
    1514           0 :         std::unique_ptr<GByte, VSIFreeReleaser> pabyTmpBuffer;
    1515             :         // Do we need a temporary buffer or can we use directly the output
    1516             :         // buffer ?
    1517         448 :         if (nBufferRadius != 0 || eDataType != eBufType ||
    1518          25 :             nPixelSpace != nBufTypeSize ||
    1519          25 :             nLineSpace != static_cast<GSpacing>(nBufTypeSize) * nBufXSize)
    1520             :         {
    1521         423 :             pabyTmpBuffer.reset(static_cast<GByte *>(VSI_CALLOC_VERBOSE(
    1522             :                 static_cast<size_t>(nExtBufXSize) * nExtBufYSize,
    1523             :                 GDALGetDataTypeSizeBytes(eDataType))));
    1524         423 :             if (!pabyTmpBuffer)
    1525           0 :                 return CE_Failure;
    1526             :         }
    1527             : 
    1528             :         {
    1529             :             const bool bUseExclusiveLock =
    1530         896 :                 m_poPrivate->m_bExclusiveLock ||
    1531         499 :                 (m_poPrivate->m_bFirstTime &&
    1532          51 :                  m_poPrivate->m_osCode.find("@jit") != std::string::npos);
    1533         448 :             m_poPrivate->m_bFirstTime = false;
    1534         448 :             GIL_Holder oHolder(bUseExclusiveLock);
    1535             : 
    1536             :             // Prepare target numpy array
    1537         448 :             PyObject *poPyDstArray = GDALCreateNumpyArray(
    1538         448 :                 m_poPrivate->m_poGDALCreateNumpyArray,
    1539         871 :                 pabyTmpBuffer ? pabyTmpBuffer.get() : pData, eDataType,
    1540             :                 nExtBufYSize, nExtBufXSize);
    1541         448 :             if (!poPyDstArray)
    1542             :             {
    1543           0 :                 return CE_Failure;
    1544             :             }
    1545             : 
    1546             :             // Wrap source buffers as input numpy arrays
    1547         448 :             PyObject *pyArgInputArray = PyTuple_New(nBufferCount);
    1548         557 :             for (int i = 0; i < nBufferCount; i++)
    1549             :             {
    1550         109 :                 GByte *pabyBuffer = static_cast<GByte *>(apBuffers[i].get());
    1551         218 :                 PyObject *poPySrcArray = GDALCreateNumpyArray(
    1552         109 :                     m_poPrivate->m_poGDALCreateNumpyArray, pabyBuffer, eSrcType,
    1553             :                     nExtBufYSize, nExtBufXSize);
    1554         109 :                 CPLAssert(poPySrcArray);
    1555         109 :                 PyTuple_SetItem(pyArgInputArray, i, poPySrcArray);
    1556             :             }
    1557             : 
    1558             :             // Create arguments
    1559         448 :             PyObject *pyArgs = PyTuple_New(10);
    1560         448 :             PyTuple_SetItem(pyArgs, 0, pyArgInputArray);
    1561         448 :             PyTuple_SetItem(pyArgs, 1, poPyDstArray);
    1562         448 :             PyTuple_SetItem(pyArgs, 2, PyLong_FromLong(nXOff));
    1563         448 :             PyTuple_SetItem(pyArgs, 3, PyLong_FromLong(nYOff));
    1564         448 :             PyTuple_SetItem(pyArgs, 4, PyLong_FromLong(nXSize));
    1565         448 :             PyTuple_SetItem(pyArgs, 5, PyLong_FromLong(nYSize));
    1566         448 :             PyTuple_SetItem(pyArgs, 6, PyLong_FromLong(nRasterXSize));
    1567         448 :             PyTuple_SetItem(pyArgs, 7, PyLong_FromLong(nRasterYSize));
    1568         448 :             PyTuple_SetItem(pyArgs, 8, PyLong_FromLong(nBufferRadius));
    1569             : 
    1570         448 :             GDALGeoTransform gt;
    1571         448 :             if (GetDataset())
    1572         448 :                 GetDataset()->GetGeoTransform(gt);
    1573         448 :             PyObject *pyGT = PyTuple_New(6);
    1574        3136 :             for (int i = 0; i < 6; i++)
    1575        2688 :                 PyTuple_SetItem(pyGT, i, PyFloat_FromDouble(gt[i]));
    1576         448 :             PyTuple_SetItem(pyArgs, 9, pyGT);
    1577             : 
    1578             :             // Prepare kwargs
    1579         448 :             PyObject *pyKwargs = PyDict_New();
    1580         616 :             for (size_t i = 0; i < m_poPrivate->m_oFunctionArgs.size(); ++i)
    1581             :             {
    1582             :                 const char *pszKey =
    1583         168 :                     m_poPrivate->m_oFunctionArgs[i].first.c_str();
    1584             :                 const char *pszValue =
    1585         168 :                     m_poPrivate->m_oFunctionArgs[i].second.c_str();
    1586         168 :                 PyDict_SetItemString(
    1587             :                     pyKwargs, pszKey,
    1588             :                     PyBytes_FromStringAndSize(pszValue, strlen(pszValue)));
    1589             :             }
    1590             : 
    1591             :             // Call user function
    1592             :             PyObject *pRetValue =
    1593         448 :                 PyObject_Call(m_poPrivate->m_poUserFunction, pyArgs, pyKwargs);
    1594             : 
    1595         448 :             Py_DecRef(pyArgs);
    1596         448 :             Py_DecRef(pyKwargs);
    1597             : 
    1598         448 :             if (ErrOccurredEmitCPLError())
    1599             :             {
    1600           2 :                 eErr = CE_Failure;
    1601             :             }
    1602         448 :             if (pRetValue)
    1603         446 :                 Py_DecRef(pRetValue);
    1604             :         }  // End of GIL section
    1605             : 
    1606         448 :         if (pabyTmpBuffer)
    1607             :         {
    1608             :             // Copy numpy destination array to user buffer
    1609       50867 :             for (int iY = 0; iY < nBufYSize; iY++)
    1610             :             {
    1611             :                 size_t nSrcOffset =
    1612       50444 :                     (static_cast<size_t>(iY + nBufferRadius) * nExtBufXSize +
    1613       50444 :                      nBufferRadius) *
    1614       50444 :                     GDALGetDataTypeSizeBytes(eDataType);
    1615       50444 :                 GDALCopyWords64(pabyTmpBuffer.get() + nSrcOffset, eDataType,
    1616             :                                 GDALGetDataTypeSizeBytes(eDataType),
    1617       50444 :                                 static_cast<GByte *>(pData) + iY * nLineSpace,
    1618             :                                 eBufType, static_cast<int>(nPixelSpace),
    1619             :                                 nBufXSize);
    1620             :             }
    1621             :         }
    1622             :     }
    1623        2831 :     else if (eErr == CE_None && poPixelFunc != nullptr)
    1624             :     {
    1625        2829 :         CPLStringList aosArgs;
    1626             : 
    1627             :         // Apply arguments specified using <PixelFunctionArguments>
    1628        4939 :         for (const auto &[pszKey, pszValue] : m_poPrivate->m_oFunctionArgs)
    1629             :         {
    1630        2110 :             aosArgs.SetNameValue(pszKey, pszValue);
    1631             :         }
    1632             : 
    1633             :         // Apply built-in arguments, potentially overwriting those in <PixelFunctionArguments>
    1634             :         // This is important because some pixel functions rely on built-in arguments being
    1635             :         // properly formatted, or even being a valid pointer. If a user can override these, we could have a crash.
    1636        4322 :         for (const auto &[pszKey, pszValue] : oAdditionalArgs)
    1637             :         {
    1638        1493 :             aosArgs.SetNameValue(pszKey, pszValue);
    1639             :         }
    1640             : 
    1641             :         static_assert(sizeof(apBuffers[0]) == sizeof(void *));
    1642        2829 :         eErr = (poPixelFunc->first)(
    1643             :             // We cast vector<unique_ptr<void>>.data() as void**. This is OK
    1644             :             // given above static_assert
    1645        2829 :             reinterpret_cast<void **>(apBuffers.data()), nBufferCount, pData,
    1646             :             nBufXSize, nBufYSize, eSrcType, eBufType,
    1647             :             static_cast<int>(nPixelSpace), static_cast<int>(nLineSpace),
    1648        2829 :             aosArgs.List());
    1649             :     }
    1650             : 
    1651        3279 :     return eErr;
    1652             : }
    1653             : 
    1654             : /************************************************************************/
    1655             : /*                       IGetDataCoverageStatus()                       */
    1656             : /************************************************************************/
    1657             : 
    1658          57 : int VRTDerivedRasterBand::IGetDataCoverageStatus(
    1659             :     int /* nXOff */, int /* nYOff */, int /* nXSize */, int /* nYSize */,
    1660             :     int /* nMaskFlagStop */, double *pdfDataPct)
    1661             : {
    1662          57 :     if (pdfDataPct != nullptr)
    1663           0 :         *pdfDataPct = -1.0;
    1664             :     return GDAL_DATA_COVERAGE_STATUS_UNIMPLEMENTED |
    1665          57 :            GDAL_DATA_COVERAGE_STATUS_DATA;
    1666             : }
    1667             : 
    1668             : /************************************************************************/
    1669             : /*                              XMLInit()                               */
    1670             : /************************************************************************/
    1671             : 
    1672        1490 : CPLErr VRTDerivedRasterBand::XMLInit(const CPLXMLNode *psTree,
    1673             :                                      const char *pszVRTPath,
    1674             :                                      VRTMapSharedResources &oMapSharedSources)
    1675             : 
    1676             : {
    1677             :     const CPLErr eErr =
    1678        1490 :         VRTSourcedRasterBand::XMLInit(psTree, pszVRTPath, oMapSharedSources);
    1679        1490 :     if (eErr != CE_None)
    1680           0 :         return eErr;
    1681             : 
    1682             :     // Read derived pixel function type.
    1683        1490 :     SetPixelFunctionName(CPLGetXMLValue(psTree, "PixelFunctionType", nullptr));
    1684        1490 :     if (osFuncName.empty())
    1685             :     {
    1686           1 :         CPLError(CE_Failure, CPLE_AppDefined, "PixelFunctionType missing");
    1687           1 :         return CE_Failure;
    1688             :     }
    1689             : 
    1690        1489 :     m_poPrivate->m_osLanguage =
    1691        1489 :         CPLGetXMLValue(psTree, "PixelFunctionLanguage", "C");
    1692        1567 :     if (!EQUAL(m_poPrivate->m_osLanguage, "C") &&
    1693          78 :         !EQUAL(m_poPrivate->m_osLanguage, "Python"))
    1694             :     {
    1695           1 :         CPLError(CE_Failure, CPLE_NotSupported,
    1696             :                  "Unsupported PixelFunctionLanguage");
    1697           1 :         return CE_Failure;
    1698             :     }
    1699             : 
    1700        1488 :     m_poPrivate->m_osCode = CPLGetXMLValue(psTree, "PixelFunctionCode", "");
    1701        1530 :     if (!m_poPrivate->m_osCode.empty() &&
    1702          42 :         !EQUAL(m_poPrivate->m_osLanguage, "Python"))
    1703             :     {
    1704           1 :         CPLError(CE_Failure, CPLE_NotSupported,
    1705             :                  "PixelFunctionCode can only be used with Python");
    1706           1 :         return CE_Failure;
    1707             :     }
    1708             : 
    1709        1487 :     m_poPrivate->m_nBufferRadius =
    1710        1487 :         atoi(CPLGetXMLValue(psTree, "BufferRadius", "0"));
    1711        1487 :     if (m_poPrivate->m_nBufferRadius < 0 || m_poPrivate->m_nBufferRadius > 1024)
    1712             :     {
    1713           1 :         CPLError(CE_Failure, CPLE_AppDefined, "Invalid value for BufferRadius");
    1714           1 :         return CE_Failure;
    1715             :     }
    1716        1500 :     if (m_poPrivate->m_nBufferRadius != 0 &&
    1717          14 :         !EQUAL(m_poPrivate->m_osLanguage, "Python"))
    1718             :     {
    1719           1 :         CPLError(CE_Failure, CPLE_NotSupported,
    1720             :                  "BufferRadius can only be used with Python");
    1721           1 :         return CE_Failure;
    1722             :     }
    1723             : 
    1724             :     const CPLXMLNode *const psArgs =
    1725        1485 :         CPLGetXMLNode(psTree, "PixelFunctionArguments");
    1726        1485 :     if (psArgs != nullptr)
    1727             :     {
    1728        2663 :         for (const CPLXMLNode *psIter = psArgs->psChild; psIter;
    1729        1410 :              psIter = psIter->psNext)
    1730             :         {
    1731        1410 :             if (psIter->eType == CXT_Attribute)
    1732             :             {
    1733        1410 :                 AddPixelFunctionArgument(psIter->pszValue,
    1734        1410 :                                          psIter->psChild->pszValue);
    1735             :             }
    1736             :         }
    1737             :     }
    1738             : 
    1739             :     // Read optional source transfer data type.
    1740             :     const char *pszTypeName =
    1741        1485 :         CPLGetXMLValue(psTree, "SourceTransferType", nullptr);
    1742        1485 :     if (pszTypeName != nullptr)
    1743             :     {
    1744         891 :         eSourceTransferType = GDALGetDataTypeByName(pszTypeName);
    1745             :     }
    1746             : 
    1747             :     // Whether to skip non contributing sources
    1748             :     const char *pszSkipNonContributingSources =
    1749        1485 :         CPLGetXMLValue(psTree, "SkipNonContributingSources", nullptr);
    1750        1485 :     if (pszSkipNonContributingSources)
    1751             :     {
    1752           2 :         SetSkipNonContributingSources(
    1753           2 :             CPLTestBool(pszSkipNonContributingSources));
    1754             :     }
    1755             : 
    1756        1485 :     return CE_None;
    1757             : }
    1758             : 
    1759             : /************************************************************************/
    1760             : /*                           SerializeToXML()                           */
    1761             : /************************************************************************/
    1762             : 
    1763          66 : CPLXMLNode *VRTDerivedRasterBand::SerializeToXML(const char *pszVRTPath,
    1764             :                                                  bool &bHasWarnedAboutRAMUsage,
    1765             :                                                  size_t &nAccRAMUsage)
    1766             : {
    1767          66 :     CPLXMLNode *psTree = VRTSourcedRasterBand::SerializeToXML(
    1768             :         pszVRTPath, bHasWarnedAboutRAMUsage, nAccRAMUsage);
    1769             : 
    1770             :     /* -------------------------------------------------------------------- */
    1771             :     /*      Set subclass.                                                   */
    1772             :     /* -------------------------------------------------------------------- */
    1773          66 :     CPLCreateXMLNode(CPLCreateXMLNode(psTree, CXT_Attribute, "subClass"),
    1774             :                      CXT_Text, "VRTDerivedRasterBand");
    1775             : 
    1776             :     /* ---- Encode DerivedBand-specific fields ---- */
    1777          66 :     if (!EQUAL(m_poPrivate->m_osLanguage, "C"))
    1778             :     {
    1779           5 :         CPLSetXMLValue(psTree, "PixelFunctionLanguage",
    1780           5 :                        m_poPrivate->m_osLanguage);
    1781             :     }
    1782          66 :     if (!osFuncName.empty())
    1783          65 :         CPLSetXMLValue(psTree, "PixelFunctionType", osFuncName.c_str());
    1784          66 :     if (!m_poPrivate->m_oFunctionArgs.empty())
    1785             :     {
    1786             :         CPLXMLNode *psArgs =
    1787          59 :             CPLCreateXMLNode(psTree, CXT_Element, "PixelFunctionArguments");
    1788         154 :         for (size_t i = 0; i < m_poPrivate->m_oFunctionArgs.size(); ++i)
    1789             :         {
    1790          95 :             const char *pszKey = m_poPrivate->m_oFunctionArgs[i].first.c_str();
    1791             :             const char *pszValue =
    1792          95 :                 m_poPrivate->m_oFunctionArgs[i].second.c_str();
    1793          95 :             CPLCreateXMLNode(CPLCreateXMLNode(psArgs, CXT_Attribute, pszKey),
    1794             :                              CXT_Text, pszValue);
    1795             :         }
    1796             :     }
    1797          66 :     if (!m_poPrivate->m_osCode.empty())
    1798             :     {
    1799           4 :         if (m_poPrivate->m_osCode.find("<![CDATA[") == std::string::npos)
    1800             :         {
    1801           4 :             CPLCreateXMLNode(
    1802             :                 CPLCreateXMLNode(psTree, CXT_Element, "PixelFunctionCode"),
    1803             :                 CXT_Literal,
    1804           8 :                 ("<![CDATA[" + m_poPrivate->m_osCode + "]]>").c_str());
    1805             :         }
    1806             :         else
    1807             :         {
    1808           0 :             CPLSetXMLValue(psTree, "PixelFunctionCode", m_poPrivate->m_osCode);
    1809             :         }
    1810             :     }
    1811          66 :     if (m_poPrivate->m_nBufferRadius != 0)
    1812           1 :         CPLSetXMLValue(psTree, "BufferRadius",
    1813           1 :                        CPLSPrintf("%d", m_poPrivate->m_nBufferRadius));
    1814          66 :     if (this->eSourceTransferType != GDT_Unknown)
    1815           4 :         CPLSetXMLValue(psTree, "SourceTransferType",
    1816             :                        GDALGetDataTypeName(eSourceTransferType));
    1817             : 
    1818          66 :     if (m_poPrivate->m_bSkipNonContributingSourcesSpecified)
    1819             :     {
    1820           1 :         CPLSetXMLValue(psTree, "SkipNonContributingSources",
    1821           1 :                        m_poPrivate->m_bSkipNonContributingSources ? "true"
    1822             :                                                                   : "false");
    1823             :     }
    1824             : 
    1825          66 :     return psTree;
    1826             : }
    1827             : 
    1828             : /************************************************************************/
    1829             : /*                             GetMinimum()                             */
    1830             : /************************************************************************/
    1831             : 
    1832           5 : double VRTDerivedRasterBand::GetMinimum(int *pbSuccess)
    1833             : {
    1834           5 :     return GDALRasterBand::GetMinimum(pbSuccess);
    1835             : }
    1836             : 
    1837             : /************************************************************************/
    1838             : /*                             GetMaximum()                             */
    1839             : /************************************************************************/
    1840             : 
    1841           5 : double VRTDerivedRasterBand::GetMaximum(int *pbSuccess)
    1842             : {
    1843           5 :     return GDALRasterBand::GetMaximum(pbSuccess);
    1844             : }
    1845             : 
    1846             : /************************************************************************/
    1847             : /*                        ComputeRasterMinMax()                         */
    1848             : /************************************************************************/
    1849             : 
    1850          15 : CPLErr VRTDerivedRasterBand::ComputeRasterMinMax(int bApproxOK,
    1851             :                                                  double *adfMinMax)
    1852             : {
    1853          15 :     return GDALRasterBand::ComputeRasterMinMax(bApproxOK, adfMinMax);
    1854             : }
    1855             : 
    1856             : /************************************************************************/
    1857             : /*                         ComputeStatistics()                          */
    1858             : /************************************************************************/
    1859             : 
    1860           1 : CPLErr VRTDerivedRasterBand::ComputeStatistics(int bApproxOK, double *pdfMin,
    1861             :                                                double *pdfMax, double *pdfMean,
    1862             :                                                double *pdfStdDev,
    1863             :                                                GDALProgressFunc pfnProgress,
    1864             :                                                void *pProgressData)
    1865             : 
    1866             : {
    1867           1 :     return GDALRasterBand::ComputeStatistics(bApproxOK, pdfMin, pdfMax, pdfMean,
    1868             :                                              pdfStdDev, pfnProgress,
    1869           1 :                                              pProgressData);
    1870             : }
    1871             : 
    1872             : /************************************************************************/
    1873             : /*                            GetHistogram()                            */
    1874             : /************************************************************************/
    1875             : 
    1876           1 : CPLErr VRTDerivedRasterBand::GetHistogram(double dfMin, double dfMax,
    1877             :                                           int nBuckets, GUIntBig *panHistogram,
    1878             :                                           int bIncludeOutOfRange, int bApproxOK,
    1879             :                                           GDALProgressFunc pfnProgress,
    1880             :                                           void *pProgressData)
    1881             : 
    1882             : {
    1883           1 :     return VRTRasterBand::GetHistogram(dfMin, dfMax, nBuckets, panHistogram,
    1884             :                                        bIncludeOutOfRange, bApproxOK,
    1885           1 :                                        pfnProgress, pProgressData);
    1886             : }
    1887             : 
    1888             : /*! @endcond */

Generated by: LCOV version 1.14