LCOV - code coverage report
Current view: top level - frmts/vrt - vrtderivedrasterband.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 664 713 93.1 %
Date: 2025-10-25 23:36:32 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         399 : 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         399 :         static_cast<size_t>(nHeight) * nWidth * GDALGetDataTypeSizeBytes(eType);
      59             :     Py_buffer pybuffer;
      60         399 :     if (PyBuffer_FillInfo(&pybuffer, nullptr, static_cast<char *>(pBuffer),
      61         399 :                           nSize, 0, PyBUF_FULL) != 0)
      62             :     {
      63           0 :         return nullptr;
      64             :     }
      65         399 :     poPyBuffer = PyMemoryView_FromBuffer(&pybuffer);
      66         399 :     PyObject *pArgsCreateArray = PyTuple_New(4);
      67         399 :     PyTuple_SetItem(pArgsCreateArray, 0, poPyBuffer);
      68         399 :     const char *pszDataType = nullptr;
      69         399 :     switch (eType)
      70             :     {
      71         357 :         case GDT_Byte:
      72         357 :             pszDataType = "uint8";
      73         357 :             break;
      74          15 :         case GDT_Int8:
      75          15 :             pszDataType = "int8";
      76          15 :             break;
      77           2 :         case GDT_UInt16:
      78           2 :             pszDataType = "uint16";
      79           2 :             break;
      80           7 :         case GDT_Int16:
      81           7 :             pszDataType = "int16";
      82           7 :             break;
      83           2 :         case GDT_UInt32:
      84           2 :             pszDataType = "uint32";
      85           2 :             break;
      86           2 :         case GDT_Int32:
      87           2 :             pszDataType = "int32";
      88           2 :             break;
      89           2 :         case GDT_Int64:
      90           2 :             pszDataType = "int64";
      91           2 :             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           2 :         case GDT_Float32:
      99           2 :             pszDataType = "float32";
     100           2 :             break;
     101           2 :         case GDT_Float64:
     102           2 :             pszDataType = "float64";
     103           2 :             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           2 :         case GDT_CFloat32:
     112           2 :             pszDataType = "complex64";
     113           2 :             break;
     114           2 :         case GDT_CFloat64:
     115           2 :             pszDataType = "complex128";
     116           2 :             break;
     117           0 :         case GDT_Unknown:
     118             :         case GDT_TypeCount:
     119           0 :             CPLAssert(FALSE);
     120             :             break;
     121             :     }
     122         399 :     PyTuple_SetItem(
     123             :         pArgsCreateArray, 1,
     124             :         PyBytes_FromStringAndSize(pszDataType, strlen(pszDataType)));
     125         399 :     PyTuple_SetItem(pArgsCreateArray, 2, PyLong_FromLong(nHeight));
     126         399 :     PyTuple_SetItem(pArgsCreateArray, 3, PyLong_FromLong(nWidth));
     127             :     PyObject *poNumpyArray =
     128         399 :         PyObject_Call(pCreateArray, pArgsCreateArray, nullptr);
     129         399 :     Py_DecRef(pArgsCreateArray);
     130         399 :     if (PyErr_Occurred())
     131           0 :         PyErr_Print();
     132         399 :     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        2218 :     VRTDerivedRasterBandPrivateData()
     164        2218 :         : 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        2218 :         const char *pszMAX_RAM = "VRT_DERIVED_DATASET_ALLOWED_RAM_USAGE";
     170        2218 :         if (const char *pszVal = CPLGetConfigOption(pszMAX_RAM, nullptr))
     171             :         {
     172           1 :             CPL_IGNORE_RET_VAL(
     173           1 :                 CPLParseMemorySize(pszVal, &m_nAllowedRAMUsage, nullptr));
     174             :         }
     175        2218 :     }
     176             : 
     177             :     ~VRTDerivedRasterBandPrivateData();
     178             : };
     179             : 
     180        2218 : VRTDerivedRasterBandPrivateData::~VRTDerivedRasterBandPrivateData()
     181             : {
     182        2218 :     if (m_poGDALCreateNumpyArray)
     183          48 :         Py_DecRef(m_poGDALCreateNumpyArray);
     184        2218 :     if (m_poUserFunction)
     185          49 :         Py_DecRef(m_poUserFunction);
     186        2218 : }
     187             : 
     188             : /************************************************************************/
     189             : /* ==================================================================== */
     190             : /*                          VRTDerivedRasterBand                        */
     191             : /* ==================================================================== */
     192             : /************************************************************************/
     193             : 
     194             : /************************************************************************/
     195             : /*                        VRTDerivedRasterBand()                        */
     196             : /************************************************************************/
     197             : 
     198        1456 : VRTDerivedRasterBand::VRTDerivedRasterBand(GDALDataset *poDSIn, int nBandIn)
     199             :     : VRTSourcedRasterBand(poDSIn, nBandIn), m_poPrivate(nullptr),
     200        1456 :       eSourceTransferType(GDT_Unknown)
     201             : {
     202        1456 :     m_poPrivate = new VRTDerivedRasterBandPrivateData;
     203        1456 : }
     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        4436 : VRTDerivedRasterBand::~VRTDerivedRasterBand()
     225             : 
     226             : {
     227        2218 :     delete m_poPrivate;
     228        4436 : }
     229             : 
     230             : /************************************************************************/
     231             : /*                               Cleanup()                              */
     232             : /************************************************************************/
     233             : 
     234        1123 : void VRTDerivedRasterBand::Cleanup()
     235             : {
     236        1123 : }
     237             : 
     238             : /************************************************************************/
     239             : /*                      GetGlobalMapPixelFunction()                     */
     240             : /************************************************************************/
     241             : 
     242             : static std::map<std::string,
     243             :                 std::pair<VRTDerivedRasterBand::PixelFunc, std::string>> &
     244       60460 : GetGlobalMapPixelFunction()
     245             : {
     246             :     static std::map<std::string,
     247             :                     std::pair<VRTDerivedRasterBand::PixelFunc, std::string>>
     248       60460 :         gosMapPixelFunction;
     249       60460 :     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       15552 : CPLErr CPL_STDCALL GDALAddDerivedBandPixelFunc(
     275             :     const char *pszName, GDALDerivedPixelFunc pfnNewFunction)
     276             : {
     277       15552 :     if (pszName == nullptr || pszName[0] == '\0' || pfnNewFunction == nullptr)
     278             :     {
     279           0 :         return CE_None;
     280             :     }
     281             : 
     282       31104 :     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       46656 :         ""};
     294             : 
     295       15552 :     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       41987 : CPLErr CPL_STDCALL GDALAddDerivedBandPixelFuncWithArgs(
     317             :     const char *pszName, GDALDerivedPixelFuncWithArgs pfnNewFunction,
     318             :     const char *pszMetadata)
     319             : {
     320       41987 :     if (!pszName || pszName[0] == '\0' || !pfnNewFunction)
     321             :     {
     322           0 :         return CE_None;
     323             :     }
     324             : 
     325       83974 :     GetGlobalMapPixelFunction()[pszName] = {pfnNewFunction,
     326      125961 :                                             pszMetadata ? pszMetadata : ""};
     327             : 
     328       41987 :     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        2825 : VRTDerivedRasterBand::GetPixelFunction(const char *pszFuncNameIn)
     379             : {
     380        2825 :     if (pszFuncNameIn == nullptr || pszFuncNameIn[0] == '\0')
     381             :     {
     382           0 :         return nullptr;
     383             :     }
     384             : 
     385        2825 :     const auto &oMapPixelFunction = GetGlobalMapPixelFunction();
     386        2825 :     const auto oIter = oMapPixelFunction.find(pszFuncNameIn);
     387             : 
     388        2825 :     if (oIter == oMapPixelFunction.end())
     389           3 :         return nullptr;
     390             : 
     391        2822 :     return &(oIter->second);
     392             : }
     393             : 
     394             : /************************************************************************/
     395             : /*                        GetPixelFunctionNames()                       */
     396             : /************************************************************************/
     397             : 
     398             : /**
     399             :  * Return the list of available pixel function names.
     400             :  */
     401             : /* static */
     402          96 : std::vector<std::string> VRTDerivedRasterBand::GetPixelFunctionNames()
     403             : {
     404          96 :     std::vector<std::string> res;
     405        3648 :     for (const auto &iter : GetGlobalMapPixelFunction())
     406             :     {
     407        3552 :         res.push_back(iter.first);
     408             :     }
     409          96 :     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        2219 : void VRTDerivedRasterBand::SetPixelFunctionName(const char *pszFuncNameIn)
     424             : {
     425        2219 :     osFuncName = (pszFuncNameIn == nullptr) ? "" : pszFuncNameIn;
     426        2219 : }
     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        2632 : void VRTDerivedRasterBand::AddPixelFunctionArgument(const char *pszArg,
     440             :                                                     const char *pszValue)
     441             : {
     442        2632 :     m_poPrivate->m_oFunctionArgs.emplace_back(pszArg, pszValue);
     443        2632 : }
     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         384 : bool VRTDerivedRasterBand::InitializePython()
     506             : {
     507         384 :     if (m_poPrivate->m_bPythonInitializationDone)
     508         321 :         return m_poPrivate->m_bPythonInitializationSuccess;
     509             : 
     510          63 :     m_poPrivate->m_bPythonInitializationDone = true;
     511          63 :     m_poPrivate->m_bPythonInitializationSuccess = false;
     512             : 
     513          63 :     const size_t nIdxDot = osFuncName.rfind(".");
     514         126 :     CPLString osPythonModule;
     515         126 :     CPLString osPythonFunction;
     516          63 :     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          34 :         osPythonFunction = osFuncName;
     524             :     }
     525             : 
     526             : #ifndef GDAL_VRT_DISABLE_PYTHON
     527             :     const char *pszPythonEnabled =
     528          63 :         CPLGetConfigOption("GDAL_VRT_ENABLE_PYTHON", nullptr);
     529             : #else
     530             :     const char *pszPythonEnabled = "NO";
     531             : #endif
     532             :     const CPLString osPythonEnabled(
     533         126 :         pszPythonEnabled ? pszPythonEnabled : GDAL_VRT_ENABLE_PYTHON_DEFAULT);
     534             : 
     535          63 :     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          52 :     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          55 :     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         106 :     m_poPrivate->m_bExclusiveLock =
     734          53 :         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         106 :         m_poPrivate->m_bExclusiveLock ||
     742          53 :         m_poPrivate->m_osCode.find("@jit") != std::string::npos;
     743         106 :     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         106 :     CPLString osModuleName(CPLSPrintf("gdal_vrt_module_%p", this));
     752         106 :     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          53 :          "\n" +
     758          53 :          m_poPrivate->m_osCode)
     759             :             .c_str(),
     760             :         osModuleName, Py_file_input);
     761          53 :     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          52 :         PyImport_ExecCodeModule(osModuleName, poCompiledString);
     769          52 :     Py_DecRef(poCompiledString);
     770             : 
     771          52 :     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          51 :     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          54 :         m_poPrivate->m_poUserFunction =
     804          27 :             PyObject_GetAttrString(poModule, osPythonFunction);
     805             :     }
     806          50 :     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          49 :     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          96 :     m_poPrivate->m_poGDALCreateNumpyArray =
     823          48 :         PyObject_GetAttrString(poModule, "GDALCreateNumpyArray");
     824          48 :     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          48 :     Py_DecRef(poModule);
     833             : 
     834          48 :     m_poPrivate->m_bPythonInitializationSuccess = true;
     835          48 :     return true;
     836             : }
     837             : 
     838        2745 : 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        5490 :     auto poArgs = CPLXMLTreeCloser(CPLParseXMLString(osMetadata));
     845        5490 :     if (poArgs != nullptr && poArgs->eType == CXT_Element &&
     846        2745 :         !strcmp(poArgs->pszValue, "PixelFunctionArgumentsList"))
     847             :     {
     848       12572 :         for (CPLXMLNode *psIter = poArgs->psChild; psIter != nullptr;
     849        9827 :              psIter = psIter->psNext)
     850             :         {
     851        9828 :             if (psIter->eType == CXT_Element &&
     852        9828 :                 !strcmp(psIter->pszValue, "Argument"))
     853             :             {
     854        9828 :                 CPLString osName, osType, osValue;
     855        9828 :                 auto pszName = CPLGetXMLValue(psIter, "name", nullptr);
     856        9828 :                 if (pszName != nullptr)
     857        5599 :                     osName = pszName;
     858        9828 :                 auto pszType = CPLGetXMLValue(psIter, "type", nullptr);
     859        9828 :                 if (pszType != nullptr)
     860        9828 :                     osType = pszType;
     861        9828 :                 auto pszValue = CPLGetXMLValue(psIter, "value", nullptr);
     862        9828 :                 if (pszValue != nullptr)
     863        4605 :                     osValue = pszValue;
     864        9828 :                 if (osType == "constant" && osValue != "" && osName != "")
     865           1 :                     oAdditionalArgs.push_back(
     866           2 :                         std::pair<CPLString, CPLString>(osName, osValue));
     867        9828 :                 if (osType == "builtin")
     868             :                 {
     869        4229 :                     const CPLString &osArgName = osValue;
     870        4229 :                     CPLString osVal;
     871        4229 :                     double dfVal = 0;
     872             : 
     873        4229 :                     int success(FALSE);
     874        4229 :                     if (osArgName == "NoData")
     875        2740 :                         dfVal = this->GetNoDataValue(&success);
     876        1489 :                     else if (osArgName == "scale")
     877           3 :                         dfVal = this->GetScale(&success);
     878        1486 :                     else if (osArgName == "offset")
     879           2 :                         dfVal = this->GetOffset(&success);
     880        1484 :                     else if (osArgName == "xoff")
     881             :                     {
     882         371 :                         dfVal = static_cast<double>(nXOff);
     883         371 :                         success = true;
     884             :                     }
     885        1113 :                     else if (osArgName == "yoff")
     886             :                     {
     887         371 :                         dfVal = static_cast<double>(nYOff);
     888         371 :                         success = true;
     889             :                     }
     890         742 :                     else if (osArgName == "geotransform")
     891             :                     {
     892         371 :                         GDALGeoTransform gt;
     893         371 :                         if (GetDataset()->GetGeoTransform(gt) != CE_None)
     894             :                         {
     895             :                             // Do not fail here because the argument is most
     896             :                             // likely not needed by the pixel function. If it
     897             :                             // is needed, the pixel function can emit the error.
     898         142 :                             continue;
     899             :                         }
     900             :                         osVal = CPLSPrintf(
     901         458 :                             "%.17g,%.17g,%.17g,%.17g,%.17g,%.17g", gt[0], gt[1],
     902         229 :                             gt[2], gt[3], gt[4], gt[5]);
     903         229 :                         success = true;
     904             :                     }
     905         371 :                     else if (osArgName == "source_names")
     906             :                     {
     907         980 :                         for (size_t iBuffer = 0;
     908         980 :                              iBuffer < anMapBufferIdxToSourceIdx.size();
     909             :                              iBuffer++)
     910             :                         {
     911             :                             const int iSource =
     912         609 :                                 anMapBufferIdxToSourceIdx[iBuffer];
     913             :                             const VRTSource *poSource =
     914         609 :                                 m_papoSources[iSource].get();
     915             : 
     916         609 :                             if (iBuffer > 0)
     917             :                             {
     918         245 :                                 osVal += "|";
     919             :                             }
     920             : 
     921         609 :                             const auto &osSourceName = poSource->GetName();
     922         609 :                             if (osSourceName.empty())
     923             :                             {
     924          42 :                                 osVal += "B" + std::to_string(iBuffer + 1);
     925             :                             }
     926             :                             else
     927             :                             {
     928         567 :                                 osVal += osSourceName;
     929             :                             }
     930             :                         }
     931             : 
     932         371 :                         success = true;
     933             :                     }
     934             :                     else
     935             :                     {
     936           0 :                         CPLError(
     937             :                             CE_Failure, CPLE_NotSupported,
     938             :                             "PixelFunction builtin argument %s not supported",
     939             :                             osArgName.c_str());
     940           0 :                         return CE_Failure;
     941             :                     }
     942        4087 :                     if (!success)
     943             :                     {
     944        2632 :                         if (CPLTestBool(
     945             :                                 CPLGetXMLValue(psIter, "optional", "false")))
     946        2631 :                             continue;
     947             : 
     948           1 :                         CPLError(CE_Failure, CPLE_AppDefined,
     949             :                                  "Raster has no %s", osValue.c_str());
     950           1 :                         return CE_Failure;
     951             :                     }
     952             : 
     953        1455 :                     if (osVal.empty())
     954             :                     {
     955         862 :                         osVal = CPLSPrintf("%.17g", dfVal);
     956             :                     }
     957             : 
     958        1455 :                     oAdditionalArgs.push_back(
     959        2910 :                         std::pair<CPLString, CPLString>(osArgName, osVal));
     960        1455 :                     CPLDebug("VRT",
     961             :                              "Added builtin pixel function argument %s = %s",
     962             :                              osArgName.c_str(), osVal.c_str());
     963             :                 }
     964             :             }
     965             :         }
     966             :     }
     967             : 
     968        2744 :     return CE_None;
     969             : }
     970             : 
     971             : /************************************************************************/
     972             : /*                             IRasterIO()                              */
     973             : /************************************************************************/
     974             : 
     975             : /**
     976             :  * Read/write a region of image data for this band.
     977             :  *
     978             :  * Each of the sources for this derived band will be read and passed to
     979             :  * the derived band pixel function.  The pixel function is responsible
     980             :  * for applying whatever algorithm is necessary to generate this band's
     981             :  * pixels from the sources.
     982             :  *
     983             :  * The sources will be read using the transfer type specified for sources
     984             :  * using SetSourceTransferType().  If no transfer type has been set for
     985             :  * this derived band, the band's data type will be used as the transfer type.
     986             :  *
     987             :  * @see gdalrasterband
     988             :  *
     989             :  * @param eRWFlag Either GF_Read to read a region of data, or GT_Write to
     990             :  * write a region of data.
     991             :  *
     992             :  * @param nXOff The pixel offset to the top left corner of the region
     993             :  * of the band to be accessed.  This would be zero to start from the left side.
     994             :  *
     995             :  * @param nYOff The line offset to the top left corner of the region
     996             :  * of the band to be accessed.  This would be zero to start from the top.
     997             :  *
     998             :  * @param nXSize The width of the region of the band to be accessed in pixels.
     999             :  *
    1000             :  * @param nYSize The height of the region of the band to be accessed in lines.
    1001             :  *
    1002             :  * @param pData The buffer into which the data should be read, or from which
    1003             :  * it should be written.  This buffer must contain at least nBufXSize *
    1004             :  * nBufYSize words of type eBufType.  It is organized in left to right,
    1005             :  * top to bottom pixel order.  Spacing is controlled by the nPixelSpace,
    1006             :  * and nLineSpace parameters.
    1007             :  *
    1008             :  * @param nBufXSize The width of the buffer image into which the desired
    1009             :  * region is to be read, or from which it is to be written.
    1010             :  *
    1011             :  * @param nBufYSize The height of the buffer image into which the desired
    1012             :  * region is to be read, or from which it is to be written.
    1013             :  *
    1014             :  * @param eBufType The type of the pixel values in the pData data buffer.  The
    1015             :  * pixel values will automatically be translated to/from the GDALRasterBand
    1016             :  * data type as needed.
    1017             :  *
    1018             :  * @param nPixelSpace The byte offset from the start of one pixel value in
    1019             :  * pData to the start of the next pixel value within a scanline.  If defaulted
    1020             :  * (0) the size of the datatype eBufType is used.
    1021             :  *
    1022             :  * @param nLineSpace The byte offset from the start of one scanline in
    1023             :  * pData to the start of the next.  If defaulted the size of the datatype
    1024             :  * eBufType * nBufXSize is used.
    1025             :  *
    1026             :  * @return CE_Failure if the access fails, otherwise CE_None.
    1027             :  */
    1028        4195 : CPLErr VRTDerivedRasterBand::IRasterIO(
    1029             :     GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
    1030             :     void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
    1031             :     GSpacing nPixelSpace, GSpacing nLineSpace, GDALRasterIOExtraArg *psExtraArg)
    1032             : {
    1033        4195 :     if (eRWFlag == GF_Write)
    1034             :     {
    1035           1 :         CPLError(CE_Failure, CPLE_AppDefined,
    1036             :                  "Writing through VRTSourcedRasterBand is not supported.");
    1037           1 :         return CE_Failure;
    1038             :     }
    1039             : 
    1040        8388 :     const std::string osFctId("VRTDerivedRasterBand::IRasterIO");
    1041        8388 :     GDALAntiRecursionGuard oGuard(osFctId);
    1042        4194 :     if (oGuard.GetCallDepth() >= 32)
    1043             :     {
    1044           0 :         CPLError(
    1045             :             CE_Failure, CPLE_AppDefined,
    1046             :             "VRTDerivedRasterBand::IRasterIO(): Recursion detected (case 1)");
    1047           0 :         return CE_Failure;
    1048             :     }
    1049             : 
    1050       12582 :     GDALAntiRecursionGuard oGuard2(oGuard, poDS->GetDescription());
    1051             :     // Allow multiple recursion depths on the same dataset in case the split strategy is applied
    1052        4194 :     if (oGuard2.GetCallDepth() > 15)
    1053             :     {
    1054           0 :         CPLError(
    1055             :             CE_Failure, CPLE_AppDefined,
    1056             :             "VRTDerivedRasterBand::IRasterIO(): Recursion detected (case 2)");
    1057           0 :         return CE_Failure;
    1058             :     }
    1059             : 
    1060             :     if constexpr (sizeof(GSpacing) > sizeof(int))
    1061             :     {
    1062        4194 :         if (nLineSpace > INT_MAX)
    1063             :         {
    1064           0 :             if (nBufYSize == 1)
    1065             :             {
    1066           0 :                 nLineSpace = 0;
    1067             :             }
    1068             :             else
    1069             :             {
    1070           0 :                 CPLError(CE_Failure, CPLE_NotSupported,
    1071             :                          "VRTDerivedRasterBand::IRasterIO(): nLineSpace > "
    1072             :                          "INT_MAX not supported");
    1073           0 :                 return CE_Failure;
    1074             :             }
    1075             :         }
    1076             :     }
    1077             : 
    1078             :     /* -------------------------------------------------------------------- */
    1079             :     /*      Do we have overviews that would be appropriate to satisfy       */
    1080             :     /*      this request?                                                   */
    1081             :     /* -------------------------------------------------------------------- */
    1082        4194 :     auto l_poDS = dynamic_cast<VRTDataset *>(poDS);
    1083        4194 :     if (l_poDS &&
    1084        8387 :         l_poDS->m_apoOverviews.empty() &&  // do not use virtual overviews
    1085        8388 :         (nBufXSize < nXSize || nBufYSize < nYSize) && GetOverviewCount() > 0)
    1086             :     {
    1087           0 :         if (OverviewRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize, pData,
    1088             :                              nBufXSize, nBufYSize, eBufType, nPixelSpace,
    1089           0 :                              nLineSpace, psExtraArg) == CE_None)
    1090           0 :             return CE_None;
    1091             :     }
    1092             : 
    1093        4194 :     const int nBufTypeSize = GDALGetDataTypeSizeBytes(eBufType);
    1094        4194 :     GDALDataType eSrcType = eSourceTransferType;
    1095        4194 :     if (eSrcType == GDT_Unknown || eSrcType >= GDT_TypeCount)
    1096             :     {
    1097             :         // Check the largest data type for all sources
    1098        3296 :         GDALDataType eAllSrcType = GDT_Unknown;
    1099      107358 :         for (auto &poSource : m_papoSources)
    1100             :         {
    1101      104079 :             if (poSource->GetType() == VRTSimpleSource::GetTypeStatic())
    1102             :             {
    1103             :                 const auto poSS =
    1104      104062 :                     static_cast<VRTSimpleSource *>(poSource.get());
    1105      104062 :                 auto l_poBand = poSS->GetRasterBand();
    1106      104062 :                 if (l_poBand)
    1107             :                 {
    1108      104062 :                     eAllSrcType = GDALDataTypeUnion(
    1109             :                         eAllSrcType, l_poBand->GetRasterDataType());
    1110             :                 }
    1111             :                 else
    1112             :                 {
    1113           0 :                     eAllSrcType = GDT_Unknown;
    1114           0 :                     break;
    1115             :                 }
    1116             :             }
    1117             :             else
    1118             :             {
    1119          17 :                 eAllSrcType = GDT_Unknown;
    1120          17 :                 break;
    1121             :             }
    1122             :         }
    1123             : 
    1124        3296 :         if (eAllSrcType != GDT_Unknown)
    1125        2898 :             eSrcType = eAllSrcType;
    1126             :         else
    1127         398 :             eSrcType = eBufType;
    1128             :     }
    1129        4194 :     const int nSrcTypeSize = GDALGetDataTypeSizeBytes(eSrcType);
    1130             : 
    1131             :     // If acquiring the region of interest in a single time is going
    1132             :     // to consume too much RAM, split in halves, and that recursively
    1133             :     // until we get below m_nAllowedRAMUsage.
    1134        4194 :     if (m_poPrivate->m_nAllowedRAMUsage > 0 && !m_papoSources.empty() &&
    1135       12195 :         nSrcTypeSize > 0 && nBufXSize == nXSize && nBufYSize == nYSize &&
    1136        3807 :         static_cast<GIntBig>(nBufXSize) * nBufYSize >
    1137        7614 :             m_poPrivate->m_nAllowedRAMUsage /
    1138        3807 :                 (static_cast<int>(m_papoSources.size()) * nSrcTypeSize))
    1139             :     {
    1140         999 :         CPLErr eErr = SplitRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
    1141             :                                     pData, nBufXSize, nBufYSize, eBufType,
    1142             :                                     nPixelSpace, nLineSpace, psExtraArg);
    1143         999 :         if (eErr != CE_Warning)
    1144         999 :             return eErr;
    1145             :     }
    1146             : 
    1147             :     /* ---- Get pixel function for band ---- */
    1148        3195 :     const std::pair<PixelFunc, std::string> *poPixelFunc = nullptr;
    1149        6390 :     std::vector<std::pair<CPLString, CPLString>> oAdditionalArgs;
    1150             : 
    1151        3195 :     if (EQUAL(m_poPrivate->m_osLanguage, "C"))
    1152             :     {
    1153             :         poPixelFunc =
    1154        2801 :             VRTDerivedRasterBand::GetPixelFunction(osFuncName.c_str());
    1155        2801 :         if (poPixelFunc == nullptr)
    1156             :         {
    1157           1 :             CPLError(CE_Failure, CPLE_IllegalArg,
    1158             :                      "VRTDerivedRasterBand::IRasterIO:"
    1159             :                      "Derived band pixel function '%s' not registered.",
    1160             :                      osFuncName.c_str());
    1161           1 :             return CE_Failure;
    1162             :         }
    1163             :     }
    1164             : 
    1165             :     /* TODO: It would be nice to use a MallocBlock function for each
    1166             :        individual buffer that would recycle blocks of memory from a
    1167             :        cache by reassigning blocks that are nearly the same size.
    1168             :        A corresponding FreeBlock might only truly free if the total size
    1169             :        of freed blocks gets to be too great of a percentage of the size
    1170             :        of the allocated blocks. */
    1171             : 
    1172             :     // Get buffers for each source.
    1173        3194 :     const int nBufferRadius = m_poPrivate->m_nBufferRadius;
    1174        3194 :     if (nBufferRadius > (INT_MAX - nBufXSize) / 2 ||
    1175        3194 :         nBufferRadius > (INT_MAX - nBufYSize) / 2)
    1176             :     {
    1177           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    1178             :                  "Integer overflow: "
    1179             :                  "nBufferRadius > (INT_MAX - nBufXSize) / 2 || "
    1180             :                  "nBufferRadius > (INT_MAX - nBufYSize) / 2)");
    1181           0 :         return CE_Failure;
    1182             :     }
    1183        3194 :     const int nExtBufXSize = nBufXSize + 2 * nBufferRadius;
    1184        3194 :     const int nExtBufYSize = nBufYSize + 2 * nBufferRadius;
    1185        3194 :     int nBufferCount = 0;
    1186             : 
    1187             :     std::vector<std::unique_ptr<void, VSIFreeReleaser>> apBuffers(
    1188        6388 :         m_papoSources.size());
    1189        6388 :     std::vector<int> anMapBufferIdxToSourceIdx(m_papoSources.size());
    1190        3194 :     bool bSkipOutputBufferInitialization = !m_papoSources.empty();
    1191      107010 :     for (int iSource = 0; iSource < static_cast<int>(m_papoSources.size());
    1192             :          iSource++)
    1193             :     {
    1194      103830 :         if (m_poPrivate->m_bSkipNonContributingSources &&
    1195          14 :             m_papoSources[iSource]->IsSimpleSource())
    1196             :         {
    1197          14 :             bool bError = false;
    1198             :             double dfReqXOff, dfReqYOff, dfReqXSize, dfReqYSize;
    1199             :             int nReqXOff, nReqYOff, nReqXSize, nReqYSize;
    1200             :             int nOutXOff, nOutYOff, nOutXSize, nOutYSize;
    1201             :             auto poSource =
    1202          14 :                 static_cast<VRTSimpleSource *>(m_papoSources[iSource].get());
    1203          14 :             if (!poSource->GetSrcDstWindow(
    1204             :                     nXOff, nYOff, nXSize, nYSize, nBufXSize, nBufYSize,
    1205             :                     &dfReqXOff, &dfReqYOff, &dfReqXSize, &dfReqYSize, &nReqXOff,
    1206             :                     &nReqYOff, &nReqXSize, &nReqYSize, &nOutXOff, &nOutYOff,
    1207             :                     &nOutXSize, &nOutYSize, bError))
    1208             :             {
    1209           4 :                 if (bError)
    1210             :                 {
    1211           0 :                     return CE_Failure;
    1212             :                 }
    1213             : 
    1214             :                 // Skip non contributing source
    1215           4 :                 bSkipOutputBufferInitialization = false;
    1216           4 :                 continue;
    1217             :             }
    1218             :         }
    1219             : 
    1220      103812 :         anMapBufferIdxToSourceIdx[nBufferCount] = iSource;
    1221      103812 :         apBuffers[nBufferCount].reset(
    1222             :             VSI_MALLOC3_VERBOSE(nSrcTypeSize, nExtBufXSize, nExtBufYSize));
    1223      103812 :         if (apBuffers[nBufferCount] == nullptr)
    1224             :         {
    1225           0 :             return CE_Failure;
    1226             :         }
    1227             : 
    1228      103812 :         bool bBufferInit = true;
    1229      103812 :         if (m_papoSources[iSource]->IsSimpleSource())
    1230             :         {
    1231             :             const auto poSS =
    1232      103811 :                 static_cast<VRTSimpleSource *>(m_papoSources[iSource].get());
    1233      103811 :             auto l_poBand = poSS->GetRasterBand();
    1234      103811 :             if (l_poBand != nullptr && poSS->m_dfSrcXOff == 0.0 &&
    1235         965 :                 poSS->m_dfSrcYOff == 0.0 &&
    1236        1930 :                 poSS->m_dfSrcXOff + poSS->m_dfSrcXSize ==
    1237         965 :                     l_poBand->GetXSize() &&
    1238        1894 :                 poSS->m_dfSrcYOff + poSS->m_dfSrcYSize ==
    1239         947 :                     l_poBand->GetYSize() &&
    1240         947 :                 poSS->m_dfDstXOff == 0.0 && poSS->m_dfDstYOff == 0.0 &&
    1241      208563 :                 poSS->m_dfDstXOff + poSS->m_dfDstXSize == nRasterXSize &&
    1242         941 :                 poSS->m_dfDstYOff + poSS->m_dfDstYSize == nRasterYSize)
    1243             :             {
    1244         941 :                 if (m_papoSources[iSource]->GetType() ==
    1245         941 :                     VRTSimpleSource::GetTypeStatic())
    1246         916 :                     bBufferInit = false;
    1247             :             }
    1248             :             else
    1249             :             {
    1250      102870 :                 bSkipOutputBufferInitialization = false;
    1251             :             }
    1252             :         }
    1253             :         else
    1254             :         {
    1255           1 :             bSkipOutputBufferInitialization = false;
    1256             :         }
    1257      103812 :         if (bBufferInit)
    1258             :         {
    1259             :             /* ------------------------------------------------------------ */
    1260             :             /* #4045: Initialize the newly allocated buffers before handing */
    1261             :             /* them off to the sources. These buffers are packed, so we     */
    1262             :             /* don't need any special line-by-line handling when a nonzero  */
    1263             :             /* nodata value is set.                                         */
    1264             :             /* ------------------------------------------------------------ */
    1265      102896 :             if (!m_bNoDataValueSet || m_dfNoDataValue == 0)
    1266             :             {
    1267      102651 :                 memset(apBuffers[nBufferCount].get(), 0,
    1268      102651 :                        static_cast<size_t>(nSrcTypeSize) * nExtBufXSize *
    1269      102651 :                            nExtBufYSize);
    1270             :             }
    1271             :             else
    1272             :             {
    1273         490 :                 GDALCopyWords64(
    1274         245 :                     &m_dfNoDataValue, GDT_Float64, 0,
    1275         245 :                     static_cast<GByte *>(apBuffers[nBufferCount].get()),
    1276             :                     eSrcType, nSrcTypeSize,
    1277         245 :                     static_cast<GPtrDiff_t>(nExtBufXSize) * nExtBufYSize);
    1278             :             }
    1279             :         }
    1280             : 
    1281      103812 :         ++nBufferCount;
    1282             :     }
    1283             : 
    1284             :     /* -------------------------------------------------------------------- */
    1285             :     /*      Initialize the buffer to some background value. Use the         */
    1286             :     /*      nodata value if available.                                      */
    1287             :     /* -------------------------------------------------------------------- */
    1288        3194 :     if (bSkipOutputBufferInitialization)
    1289             :     {
    1290             :         // Do nothing
    1291             :     }
    1292        2573 :     else if (nPixelSpace == nBufTypeSize &&
    1293        2555 :              (!m_bNoDataValueSet || m_dfNoDataValue == 0))
    1294             :     {
    1295        2469 :         memset(pData, 0,
    1296        2469 :                static_cast<size_t>(nBufXSize) * nBufYSize * nBufTypeSize);
    1297             :     }
    1298         104 :     else if (m_bNoDataValueSet)
    1299             :     {
    1300          86 :         double dfWriteValue = m_dfNoDataValue;
    1301             : 
    1302         237 :         for (int iLine = 0; iLine < nBufYSize; iLine++)
    1303             :         {
    1304         151 :             GDALCopyWords64(&dfWriteValue, GDT_Float64, 0,
    1305         151 :                             static_cast<GByte *>(pData) + nLineSpace * iLine,
    1306             :                             eBufType, static_cast<int>(nPixelSpace), nBufXSize);
    1307             :         }
    1308             :     }
    1309             : 
    1310             :     // No contributing sources and SkipNonContributingSources mode ?
    1311             :     // Do not call the pixel function and just return the 0/nodata initialized
    1312             :     // output buffer.
    1313        3194 :     if (nBufferCount == 0 && m_poPrivate->m_bSkipNonContributingSources)
    1314             :     {
    1315           1 :         return CE_None;
    1316             :     }
    1317             : 
    1318             :     GDALRasterIOExtraArg sExtraArg;
    1319        3193 :     GDALCopyRasterIOExtraArg(&sExtraArg, psExtraArg);
    1320             : 
    1321        3193 :     int nXShiftInBuffer = 0;
    1322        3193 :     int nYShiftInBuffer = 0;
    1323        3193 :     int nExtBufXSizeReq = nExtBufXSize;
    1324        3193 :     int nExtBufYSizeReq = nExtBufYSize;
    1325             : 
    1326        3193 :     int nXOffExt = nXOff;
    1327        3193 :     int nYOffExt = nYOff;
    1328        3193 :     int nXSizeExt = nXSize;
    1329        3193 :     int nYSizeExt = nYSize;
    1330             : 
    1331        3193 :     if (nBufferRadius)
    1332             :     {
    1333           9 :         double dfXRatio = static_cast<double>(nXSize) / nBufXSize;
    1334           9 :         double dfYRatio = static_cast<double>(nYSize) / nBufYSize;
    1335             : 
    1336           9 :         if (!sExtraArg.bFloatingPointWindowValidity)
    1337             :         {
    1338           9 :             sExtraArg.dfXOff = nXOff;
    1339           9 :             sExtraArg.dfYOff = nYOff;
    1340           9 :             sExtraArg.dfXSize = nXSize;
    1341           9 :             sExtraArg.dfYSize = nYSize;
    1342             :         }
    1343             : 
    1344           9 :         sExtraArg.dfXOff -= dfXRatio * nBufferRadius;
    1345           9 :         sExtraArg.dfYOff -= dfYRatio * nBufferRadius;
    1346           9 :         sExtraArg.dfXSize += 2 * dfXRatio * nBufferRadius;
    1347           9 :         sExtraArg.dfYSize += 2 * dfYRatio * nBufferRadius;
    1348           9 :         if (sExtraArg.dfXOff < 0)
    1349             :         {
    1350           9 :             nXShiftInBuffer = -static_cast<int>(sExtraArg.dfXOff / dfXRatio);
    1351           9 :             nExtBufXSizeReq -= nXShiftInBuffer;
    1352           9 :             sExtraArg.dfXSize += sExtraArg.dfXOff;
    1353           9 :             sExtraArg.dfXOff = 0;
    1354             :         }
    1355           9 :         if (sExtraArg.dfYOff < 0)
    1356             :         {
    1357           9 :             nYShiftInBuffer = -static_cast<int>(sExtraArg.dfYOff / dfYRatio);
    1358           9 :             nExtBufYSizeReq -= nYShiftInBuffer;
    1359           9 :             sExtraArg.dfYSize += sExtraArg.dfYOff;
    1360           9 :             sExtraArg.dfYOff = 0;
    1361             :         }
    1362           9 :         if (sExtraArg.dfXOff + sExtraArg.dfXSize > nRasterXSize)
    1363             :         {
    1364           9 :             nExtBufXSizeReq -= static_cast<int>(
    1365           9 :                 (sExtraArg.dfXOff + sExtraArg.dfXSize - nRasterXSize) /
    1366             :                 dfXRatio);
    1367           9 :             sExtraArg.dfXSize = nRasterXSize - sExtraArg.dfXOff;
    1368             :         }
    1369           9 :         if (sExtraArg.dfYOff + sExtraArg.dfYSize > nRasterYSize)
    1370             :         {
    1371           9 :             nExtBufYSizeReq -= static_cast<int>(
    1372           9 :                 (sExtraArg.dfYOff + sExtraArg.dfYSize - nRasterYSize) /
    1373             :                 dfYRatio);
    1374           9 :             sExtraArg.dfYSize = nRasterYSize - sExtraArg.dfYOff;
    1375             :         }
    1376             : 
    1377           9 :         nXOffExt = static_cast<int>(sExtraArg.dfXOff);
    1378           9 :         nYOffExt = static_cast<int>(sExtraArg.dfYOff);
    1379          18 :         nXSizeExt = std::min(static_cast<int>(sExtraArg.dfXSize + 0.5),
    1380           9 :                              nRasterXSize - nXOffExt);
    1381          18 :         nYSizeExt = std::min(static_cast<int>(sExtraArg.dfYSize + 0.5),
    1382           9 :                              nRasterYSize - nYOffExt);
    1383             :     }
    1384             : 
    1385             :     // Load values for sources into packed buffers.
    1386        3193 :     CPLErr eErr = CE_None;
    1387        6386 :     VRTSource::WorkingState oWorkingState;
    1388      107005 :     for (int iBuffer = 0; iBuffer < nBufferCount && eErr == CE_None; iBuffer++)
    1389             :     {
    1390      103812 :         const int iSource = anMapBufferIdxToSourceIdx[iBuffer];
    1391      103812 :         GByte *pabyBuffer = static_cast<GByte *>(apBuffers[iBuffer].get());
    1392      103812 :         eErr = static_cast<VRTSource *>(m_papoSources[iSource].get())
    1393      207624 :                    ->RasterIO(
    1394             :                        eSrcType, nXOffExt, nYOffExt, nXSizeExt, nYSizeExt,
    1395      103812 :                        pabyBuffer + (static_cast<size_t>(nYShiftInBuffer) *
    1396      103812 :                                          nExtBufXSize +
    1397      103812 :                                      nXShiftInBuffer) *
    1398      103812 :                                         nSrcTypeSize,
    1399             :                        nExtBufXSizeReq, nExtBufYSizeReq, eSrcType, nSrcTypeSize,
    1400      103812 :                        static_cast<GSpacing>(nSrcTypeSize) * nExtBufXSize,
    1401      103812 :                        &sExtraArg, oWorkingState);
    1402             : 
    1403             :         // Extend first lines
    1404      103821 :         for (int iY = 0; iY < nYShiftInBuffer; iY++)
    1405             :         {
    1406           9 :             memcpy(pabyBuffer +
    1407           9 :                        static_cast<size_t>(iY) * nExtBufXSize * nSrcTypeSize,
    1408           9 :                    pabyBuffer + static_cast<size_t>(nYShiftInBuffer) *
    1409           9 :                                     nExtBufXSize * nSrcTypeSize,
    1410           9 :                    static_cast<size_t>(nExtBufXSize) * nSrcTypeSize);
    1411             :         }
    1412             :         // Extend last lines
    1413      103821 :         for (int iY = nYShiftInBuffer + nExtBufYSizeReq; iY < nExtBufYSize;
    1414             :              iY++)
    1415             :         {
    1416           9 :             memcpy(pabyBuffer +
    1417           9 :                        static_cast<size_t>(iY) * nExtBufXSize * nSrcTypeSize,
    1418           9 :                    pabyBuffer + static_cast<size_t>(nYShiftInBuffer +
    1419           9 :                                                     nExtBufYSizeReq - 1) *
    1420           9 :                                     nExtBufXSize * nSrcTypeSize,
    1421           9 :                    static_cast<size_t>(nExtBufXSize) * nSrcTypeSize);
    1422             :         }
    1423             :         // Extend first cols
    1424      103812 :         if (nXShiftInBuffer)
    1425             :         {
    1426        1116 :             for (int iY = 0; iY < nExtBufYSize; iY++)
    1427             :             {
    1428        2214 :                 for (int iX = 0; iX < nXShiftInBuffer; iX++)
    1429             :                 {
    1430        1107 :                     memcpy(pabyBuffer +
    1431        1107 :                                static_cast<size_t>(iY * nExtBufXSize + iX) *
    1432        1107 :                                    nSrcTypeSize,
    1433        1107 :                            pabyBuffer +
    1434        1107 :                                (static_cast<size_t>(iY) * nExtBufXSize +
    1435        1107 :                                 nXShiftInBuffer) *
    1436        1107 :                                    nSrcTypeSize,
    1437             :                            nSrcTypeSize);
    1438             :                 }
    1439             :             }
    1440             :         }
    1441             :         // Extent last cols
    1442      103812 :         if (nXShiftInBuffer + nExtBufXSizeReq < nExtBufXSize)
    1443             :         {
    1444        1116 :             for (int iY = 0; iY < nExtBufYSize; iY++)
    1445             :             {
    1446        1107 :                 for (int iX = nXShiftInBuffer + nExtBufXSizeReq;
    1447        2214 :                      iX < nExtBufXSize; iX++)
    1448             :                 {
    1449        1107 :                     memcpy(pabyBuffer +
    1450        1107 :                                (static_cast<size_t>(iY) * nExtBufXSize + iX) *
    1451        1107 :                                    nSrcTypeSize,
    1452        1107 :                            pabyBuffer +
    1453        1107 :                                (static_cast<size_t>(iY) * nExtBufXSize +
    1454        1107 :                                 nXShiftInBuffer + nExtBufXSizeReq - 1) *
    1455        1107 :                                    nSrcTypeSize,
    1456             :                            nSrcTypeSize);
    1457             :                 }
    1458             :             }
    1459             :         }
    1460             :     }
    1461             : 
    1462             :     // Collect any pixel function arguments
    1463        3193 :     if (poPixelFunc != nullptr && !poPixelFunc->second.empty())
    1464             :     {
    1465        5490 :         if (GetPixelFunctionArguments(poPixelFunc->second,
    1466             :                                       anMapBufferIdxToSourceIdx, nXOff, nYOff,
    1467        2745 :                                       oAdditionalArgs) != CE_None)
    1468             :         {
    1469           1 :             eErr = CE_Failure;
    1470             :         }
    1471             :     }
    1472             : 
    1473             :     // Apply pixel function.
    1474        3193 :     if (eErr == CE_None && EQUAL(m_poPrivate->m_osLanguage, "Python"))
    1475             :     {
    1476             :         // numpy doesn't have native cint16/cint32/cfloat16
    1477         393 :         if (eSrcType == GDT_CInt16 || eSrcType == GDT_CInt32 ||
    1478             :             eSrcType == GDT_CFloat16)
    1479             :         {
    1480           2 :             CPLError(CE_Failure, CPLE_AppDefined,
    1481             :                      "CInt16/CInt32/CFloat16 data type not supported for "
    1482             :                      "SourceTransferType");
    1483          24 :             return CE_Failure;
    1484             :         }
    1485         391 :         if (eDataType == GDT_CInt16 || eDataType == GDT_CInt32 ||
    1486         387 :             eDataType == GDT_CFloat16)
    1487             :         {
    1488           7 :             CPLError(
    1489             :                 CE_Failure, CPLE_AppDefined,
    1490             :                 "CInt16/CInt32/CFloat16 data type not supported for data type");
    1491           7 :             return CE_Failure;
    1492             :         }
    1493             : 
    1494         384 :         if (!InitializePython())
    1495          15 :             return CE_Failure;
    1496             : 
    1497           0 :         std::unique_ptr<GByte, VSIFreeReleaser> pabyTmpBuffer;
    1498             :         // Do we need a temporary buffer or can we use directly the output
    1499             :         // buffer ?
    1500         369 :         if (nBufferRadius != 0 || eDataType != eBufType ||
    1501          25 :             nPixelSpace != nBufTypeSize ||
    1502          25 :             nLineSpace != static_cast<GSpacing>(nBufTypeSize) * nBufXSize)
    1503             :         {
    1504         344 :             pabyTmpBuffer.reset(static_cast<GByte *>(VSI_CALLOC_VERBOSE(
    1505             :                 static_cast<size_t>(nExtBufXSize) * nExtBufYSize,
    1506             :                 GDALGetDataTypeSizeBytes(eDataType))));
    1507         344 :             if (!pabyTmpBuffer)
    1508           0 :                 return CE_Failure;
    1509             :         }
    1510             : 
    1511             :         {
    1512             :             const bool bUseExclusiveLock =
    1513         738 :                 m_poPrivate->m_bExclusiveLock ||
    1514         417 :                 (m_poPrivate->m_bFirstTime &&
    1515          48 :                  m_poPrivate->m_osCode.find("@jit") != std::string::npos);
    1516         369 :             m_poPrivate->m_bFirstTime = false;
    1517         369 :             GIL_Holder oHolder(bUseExclusiveLock);
    1518             : 
    1519             :             // Prepare target numpy array
    1520         369 :             PyObject *poPyDstArray = GDALCreateNumpyArray(
    1521         369 :                 m_poPrivate->m_poGDALCreateNumpyArray,
    1522         713 :                 pabyTmpBuffer ? pabyTmpBuffer.get() : pData, eDataType,
    1523             :                 nExtBufYSize, nExtBufXSize);
    1524         369 :             if (!poPyDstArray)
    1525             :             {
    1526           0 :                 return CE_Failure;
    1527             :             }
    1528             : 
    1529             :             // Wrap source buffers as input numpy arrays
    1530         369 :             PyObject *pyArgInputArray = PyTuple_New(nBufferCount);
    1531         399 :             for (int i = 0; i < nBufferCount; i++)
    1532             :             {
    1533          30 :                 GByte *pabyBuffer = static_cast<GByte *>(apBuffers[i].get());
    1534          60 :                 PyObject *poPySrcArray = GDALCreateNumpyArray(
    1535          30 :                     m_poPrivate->m_poGDALCreateNumpyArray, pabyBuffer, eSrcType,
    1536             :                     nExtBufYSize, nExtBufXSize);
    1537          30 :                 CPLAssert(poPySrcArray);
    1538          30 :                 PyTuple_SetItem(pyArgInputArray, i, poPySrcArray);
    1539             :             }
    1540             : 
    1541             :             // Create arguments
    1542         369 :             PyObject *pyArgs = PyTuple_New(10);
    1543         369 :             PyTuple_SetItem(pyArgs, 0, pyArgInputArray);
    1544         369 :             PyTuple_SetItem(pyArgs, 1, poPyDstArray);
    1545         369 :             PyTuple_SetItem(pyArgs, 2, PyLong_FromLong(nXOff));
    1546         369 :             PyTuple_SetItem(pyArgs, 3, PyLong_FromLong(nYOff));
    1547         369 :             PyTuple_SetItem(pyArgs, 4, PyLong_FromLong(nXSize));
    1548         369 :             PyTuple_SetItem(pyArgs, 5, PyLong_FromLong(nYSize));
    1549         369 :             PyTuple_SetItem(pyArgs, 6, PyLong_FromLong(nRasterXSize));
    1550         369 :             PyTuple_SetItem(pyArgs, 7, PyLong_FromLong(nRasterYSize));
    1551         369 :             PyTuple_SetItem(pyArgs, 8, PyLong_FromLong(nBufferRadius));
    1552             : 
    1553         369 :             GDALGeoTransform gt;
    1554         369 :             if (GetDataset())
    1555         369 :                 GetDataset()->GetGeoTransform(gt);
    1556         369 :             PyObject *pyGT = PyTuple_New(6);
    1557        2583 :             for (int i = 0; i < 6; i++)
    1558        2214 :                 PyTuple_SetItem(pyGT, i, PyFloat_FromDouble(gt[i]));
    1559         369 :             PyTuple_SetItem(pyArgs, 9, pyGT);
    1560             : 
    1561             :             // Prepare kwargs
    1562         369 :             PyObject *pyKwargs = PyDict_New();
    1563         379 :             for (size_t i = 0; i < m_poPrivate->m_oFunctionArgs.size(); ++i)
    1564             :             {
    1565             :                 const char *pszKey =
    1566          10 :                     m_poPrivate->m_oFunctionArgs[i].first.c_str();
    1567             :                 const char *pszValue =
    1568          10 :                     m_poPrivate->m_oFunctionArgs[i].second.c_str();
    1569          10 :                 PyDict_SetItemString(
    1570             :                     pyKwargs, pszKey,
    1571             :                     PyBytes_FromStringAndSize(pszValue, strlen(pszValue)));
    1572             :             }
    1573             : 
    1574             :             // Call user function
    1575             :             PyObject *pRetValue =
    1576         369 :                 PyObject_Call(m_poPrivate->m_poUserFunction, pyArgs, pyKwargs);
    1577             : 
    1578         369 :             Py_DecRef(pyArgs);
    1579         369 :             Py_DecRef(pyKwargs);
    1580             : 
    1581         369 :             if (ErrOccurredEmitCPLError())
    1582             :             {
    1583           2 :                 eErr = CE_Failure;
    1584             :             }
    1585         369 :             if (pRetValue)
    1586         367 :                 Py_DecRef(pRetValue);
    1587             :         }  // End of GIL section
    1588             : 
    1589         369 :         if (pabyTmpBuffer)
    1590             :         {
    1591             :             // Copy numpy destination array to user buffer
    1592       41219 :             for (int iY = 0; iY < nBufYSize; iY++)
    1593             :             {
    1594             :                 size_t nSrcOffset =
    1595       40873 :                     (static_cast<size_t>(iY + nBufferRadius) * nExtBufXSize +
    1596       40873 :                      nBufferRadius) *
    1597       40873 :                     GDALGetDataTypeSizeBytes(eDataType);
    1598       40869 :                 GDALCopyWords64(pabyTmpBuffer.get() + nSrcOffset, eDataType,
    1599             :                                 GDALGetDataTypeSizeBytes(eDataType),
    1600       40869 :                                 static_cast<GByte *>(pData) + iY * nLineSpace,
    1601             :                                 eBufType, static_cast<int>(nPixelSpace),
    1602             :                                 nBufXSize);
    1603             :             }
    1604             :         }
    1605             :     }
    1606        2800 :     else if (eErr == CE_None && poPixelFunc != nullptr)
    1607             :     {
    1608        2799 :         CPLStringList aosArgs;
    1609             : 
    1610        2799 :         oAdditionalArgs.insert(oAdditionalArgs.end(),
    1611        2799 :                                m_poPrivate->m_oFunctionArgs.begin(),
    1612        5598 :                                m_poPrivate->m_oFunctionArgs.end());
    1613        6339 :         for (const auto &oArg : oAdditionalArgs)
    1614             :         {
    1615        3540 :             const char *pszKey = oArg.first.c_str();
    1616        3540 :             const char *pszValue = oArg.second.c_str();
    1617        3540 :             aosArgs.SetNameValue(pszKey, pszValue);
    1618             :         }
    1619             : 
    1620             :         static_assert(sizeof(apBuffers[0]) == sizeof(void *));
    1621        2799 :         eErr = (poPixelFunc->first)(
    1622             :             // We cast vector<unique_ptr<void>>.data() as void**. This is OK
    1623             :             // given above static_assert
    1624        2799 :             reinterpret_cast<void **>(apBuffers.data()), nBufferCount, pData,
    1625             :             nBufXSize, nBufYSize, eSrcType, eBufType,
    1626             :             static_cast<int>(nPixelSpace), static_cast<int>(nLineSpace),
    1627        2799 :             aosArgs.List());
    1628             :     }
    1629             : 
    1630        3169 :     return eErr;
    1631             : }
    1632             : 
    1633             : /************************************************************************/
    1634             : /*                         IGetDataCoverageStatus()                     */
    1635             : /************************************************************************/
    1636             : 
    1637          57 : int VRTDerivedRasterBand::IGetDataCoverageStatus(
    1638             :     int /* nXOff */, int /* nYOff */, int /* nXSize */, int /* nYSize */,
    1639             :     int /* nMaskFlagStop */, double *pdfDataPct)
    1640             : {
    1641          57 :     if (pdfDataPct != nullptr)
    1642           0 :         *pdfDataPct = -1.0;
    1643             :     return GDAL_DATA_COVERAGE_STATUS_UNIMPLEMENTED |
    1644          57 :            GDAL_DATA_COVERAGE_STATUS_DATA;
    1645             : }
    1646             : 
    1647             : /************************************************************************/
    1648             : /*                              XMLInit()                               */
    1649             : /************************************************************************/
    1650             : 
    1651        1456 : CPLErr VRTDerivedRasterBand::XMLInit(const CPLXMLNode *psTree,
    1652             :                                      const char *pszVRTPath,
    1653             :                                      VRTMapSharedResources &oMapSharedSources)
    1654             : 
    1655             : {
    1656             :     const CPLErr eErr =
    1657        1456 :         VRTSourcedRasterBand::XMLInit(psTree, pszVRTPath, oMapSharedSources);
    1658        1456 :     if (eErr != CE_None)
    1659           0 :         return eErr;
    1660             : 
    1661             :     // Read derived pixel function type.
    1662        1456 :     SetPixelFunctionName(CPLGetXMLValue(psTree, "PixelFunctionType", nullptr));
    1663        1456 :     if (osFuncName.empty())
    1664             :     {
    1665           1 :         CPLError(CE_Failure, CPLE_AppDefined, "PixelFunctionType missing");
    1666           1 :         return CE_Failure;
    1667             :     }
    1668             : 
    1669        1455 :     m_poPrivate->m_osLanguage =
    1670        1455 :         CPLGetXMLValue(psTree, "PixelFunctionLanguage", "C");
    1671        1530 :     if (!EQUAL(m_poPrivate->m_osLanguage, "C") &&
    1672          75 :         !EQUAL(m_poPrivate->m_osLanguage, "Python"))
    1673             :     {
    1674           1 :         CPLError(CE_Failure, CPLE_NotSupported,
    1675             :                  "Unsupported PixelFunctionLanguage");
    1676           1 :         return CE_Failure;
    1677             :     }
    1678             : 
    1679        1454 :     m_poPrivate->m_osCode = CPLGetXMLValue(psTree, "PixelFunctionCode", "");
    1680        1493 :     if (!m_poPrivate->m_osCode.empty() &&
    1681          39 :         !EQUAL(m_poPrivate->m_osLanguage, "Python"))
    1682             :     {
    1683           1 :         CPLError(CE_Failure, CPLE_NotSupported,
    1684             :                  "PixelFunctionCode can only be used with Python");
    1685           1 :         return CE_Failure;
    1686             :     }
    1687             : 
    1688        1453 :     m_poPrivate->m_nBufferRadius =
    1689        1453 :         atoi(CPLGetXMLValue(psTree, "BufferRadius", "0"));
    1690        1453 :     if (m_poPrivate->m_nBufferRadius < 0 || m_poPrivate->m_nBufferRadius > 1024)
    1691             :     {
    1692           1 :         CPLError(CE_Failure, CPLE_AppDefined, "Invalid value for BufferRadius");
    1693           1 :         return CE_Failure;
    1694             :     }
    1695        1463 :     if (m_poPrivate->m_nBufferRadius != 0 &&
    1696          11 :         !EQUAL(m_poPrivate->m_osLanguage, "Python"))
    1697             :     {
    1698           1 :         CPLError(CE_Failure, CPLE_NotSupported,
    1699             :                  "BufferRadius can only be used with Python");
    1700           1 :         return CE_Failure;
    1701             :     }
    1702             : 
    1703             :     const CPLXMLNode *const psArgs =
    1704        1451 :         CPLGetXMLNode(psTree, "PixelFunctionArguments");
    1705        1451 :     if (psArgs != nullptr)
    1706             :     {
    1707        2598 :         for (const CPLXMLNode *psIter = psArgs->psChild; psIter;
    1708        1378 :              psIter = psIter->psNext)
    1709             :         {
    1710        1378 :             if (psIter->eType == CXT_Attribute)
    1711             :             {
    1712        1378 :                 AddPixelFunctionArgument(psIter->pszValue,
    1713        1378 :                                          psIter->psChild->pszValue);
    1714             :             }
    1715             :         }
    1716             :     }
    1717             : 
    1718             :     // Read optional source transfer data type.
    1719             :     const char *pszTypeName =
    1720        1451 :         CPLGetXMLValue(psTree, "SourceTransferType", nullptr);
    1721        1451 :     if (pszTypeName != nullptr)
    1722             :     {
    1723         888 :         eSourceTransferType = GDALGetDataTypeByName(pszTypeName);
    1724             :     }
    1725             : 
    1726             :     // Whether to skip non contributing sources
    1727             :     const char *pszSkipNonContributingSources =
    1728        1451 :         CPLGetXMLValue(psTree, "SkipNonContributingSources", nullptr);
    1729        1451 :     if (pszSkipNonContributingSources)
    1730             :     {
    1731           2 :         SetSkipNonContributingSources(
    1732           2 :             CPLTestBool(pszSkipNonContributingSources));
    1733             :     }
    1734             : 
    1735        1451 :     return CE_None;
    1736             : }
    1737             : 
    1738             : /************************************************************************/
    1739             : /*                           SerializeToXML()                           */
    1740             : /************************************************************************/
    1741             : 
    1742          66 : CPLXMLNode *VRTDerivedRasterBand::SerializeToXML(const char *pszVRTPath,
    1743             :                                                  bool &bHasWarnedAboutRAMUsage,
    1744             :                                                  size_t &nAccRAMUsage)
    1745             : {
    1746          66 :     CPLXMLNode *psTree = VRTSourcedRasterBand::SerializeToXML(
    1747             :         pszVRTPath, bHasWarnedAboutRAMUsage, nAccRAMUsage);
    1748             : 
    1749             :     /* -------------------------------------------------------------------- */
    1750             :     /*      Set subclass.                                                   */
    1751             :     /* -------------------------------------------------------------------- */
    1752          66 :     CPLCreateXMLNode(CPLCreateXMLNode(psTree, CXT_Attribute, "subClass"),
    1753             :                      CXT_Text, "VRTDerivedRasterBand");
    1754             : 
    1755             :     /* ---- Encode DerivedBand-specific fields ---- */
    1756          66 :     if (!EQUAL(m_poPrivate->m_osLanguage, "C"))
    1757             :     {
    1758           5 :         CPLSetXMLValue(psTree, "PixelFunctionLanguage",
    1759           5 :                        m_poPrivate->m_osLanguage);
    1760             :     }
    1761          66 :     if (!osFuncName.empty())
    1762          65 :         CPLSetXMLValue(psTree, "PixelFunctionType", osFuncName.c_str());
    1763          66 :     if (!m_poPrivate->m_oFunctionArgs.empty())
    1764             :     {
    1765             :         CPLXMLNode *psArgs =
    1766          59 :             CPLCreateXMLNode(psTree, CXT_Element, "PixelFunctionArguments");
    1767         154 :         for (size_t i = 0; i < m_poPrivate->m_oFunctionArgs.size(); ++i)
    1768             :         {
    1769          95 :             const char *pszKey = m_poPrivate->m_oFunctionArgs[i].first.c_str();
    1770             :             const char *pszValue =
    1771          95 :                 m_poPrivate->m_oFunctionArgs[i].second.c_str();
    1772          95 :             CPLCreateXMLNode(CPLCreateXMLNode(psArgs, CXT_Attribute, pszKey),
    1773             :                              CXT_Text, pszValue);
    1774             :         }
    1775             :     }
    1776          66 :     if (!m_poPrivate->m_osCode.empty())
    1777             :     {
    1778           4 :         if (m_poPrivate->m_osCode.find("<![CDATA[") == std::string::npos)
    1779             :         {
    1780           4 :             CPLCreateXMLNode(
    1781             :                 CPLCreateXMLNode(psTree, CXT_Element, "PixelFunctionCode"),
    1782             :                 CXT_Literal,
    1783           8 :                 ("<![CDATA[" + m_poPrivate->m_osCode + "]]>").c_str());
    1784             :         }
    1785             :         else
    1786             :         {
    1787           0 :             CPLSetXMLValue(psTree, "PixelFunctionCode", m_poPrivate->m_osCode);
    1788             :         }
    1789             :     }
    1790          66 :     if (m_poPrivate->m_nBufferRadius != 0)
    1791           1 :         CPLSetXMLValue(psTree, "BufferRadius",
    1792           1 :                        CPLSPrintf("%d", m_poPrivate->m_nBufferRadius));
    1793          66 :     if (this->eSourceTransferType != GDT_Unknown)
    1794           4 :         CPLSetXMLValue(psTree, "SourceTransferType",
    1795             :                        GDALGetDataTypeName(eSourceTransferType));
    1796             : 
    1797          66 :     if (m_poPrivate->m_bSkipNonContributingSourcesSpecified)
    1798             :     {
    1799           1 :         CPLSetXMLValue(psTree, "SkipNonContributingSources",
    1800           1 :                        m_poPrivate->m_bSkipNonContributingSources ? "true"
    1801             :                                                                   : "false");
    1802             :     }
    1803             : 
    1804          66 :     return psTree;
    1805             : }
    1806             : 
    1807             : /************************************************************************/
    1808             : /*                             GetMinimum()                             */
    1809             : /************************************************************************/
    1810             : 
    1811           5 : double VRTDerivedRasterBand::GetMinimum(int *pbSuccess)
    1812             : {
    1813           5 :     return GDALRasterBand::GetMinimum(pbSuccess);
    1814             : }
    1815             : 
    1816             : /************************************************************************/
    1817             : /*                             GetMaximum()                             */
    1818             : /************************************************************************/
    1819             : 
    1820           5 : double VRTDerivedRasterBand::GetMaximum(int *pbSuccess)
    1821             : {
    1822           5 :     return GDALRasterBand::GetMaximum(pbSuccess);
    1823             : }
    1824             : 
    1825             : /************************************************************************/
    1826             : /*                       ComputeRasterMinMax()                          */
    1827             : /************************************************************************/
    1828             : 
    1829          15 : CPLErr VRTDerivedRasterBand::ComputeRasterMinMax(int bApproxOK,
    1830             :                                                  double *adfMinMax)
    1831             : {
    1832          15 :     return GDALRasterBand::ComputeRasterMinMax(bApproxOK, adfMinMax);
    1833             : }
    1834             : 
    1835             : /************************************************************************/
    1836             : /*                         ComputeStatistics()                          */
    1837             : /************************************************************************/
    1838             : 
    1839           1 : CPLErr VRTDerivedRasterBand::ComputeStatistics(int bApproxOK, double *pdfMin,
    1840             :                                                double *pdfMax, double *pdfMean,
    1841             :                                                double *pdfStdDev,
    1842             :                                                GDALProgressFunc pfnProgress,
    1843             :                                                void *pProgressData)
    1844             : 
    1845             : {
    1846           1 :     return GDALRasterBand::ComputeStatistics(bApproxOK, pdfMin, pdfMax, pdfMean,
    1847             :                                              pdfStdDev, pfnProgress,
    1848           1 :                                              pProgressData);
    1849             : }
    1850             : 
    1851             : /************************************************************************/
    1852             : /*                            GetHistogram()                            */
    1853             : /************************************************************************/
    1854             : 
    1855           1 : CPLErr VRTDerivedRasterBand::GetHistogram(double dfMin, double dfMax,
    1856             :                                           int nBuckets, GUIntBig *panHistogram,
    1857             :                                           int bIncludeOutOfRange, int bApproxOK,
    1858             :                                           GDALProgressFunc pfnProgress,
    1859             :                                           void *pProgressData)
    1860             : 
    1861             : {
    1862           1 :     return VRTRasterBand::GetHistogram(dfMin, dfMax, nBuckets, panHistogram,
    1863             :                                        bIncludeOutOfRange, bApproxOK,
    1864           1 :                                        pfnProgress, pProgressData);
    1865             : }
    1866             : 
    1867             : /*! @endcond */

Generated by: LCOV version 1.14