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

Generated by: LCOV version 1.14