/******************************************************************************
 *
 * Project:  Virtual GDAL Datasets
 * Purpose:  Implementation of a sourced raster band that derives its raster
 *           by applying an algorithm (GDALDerivedPixelFunc) to the sources.
 * Author:   Pete Nagy
 *
 ******************************************************************************
 * Copyright (c) 2005 Vexcel Corp.
 * Copyright (c) 2008-2011, Even Rouault <even dot rouault at mines-paris dot org>
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 *****************************************************************************/

#include "cpl_minixml.h"
#include "cpl_string.h"
#include "vrtdataset.h"
#include "cpl_multiproc.h"
#include "cpl_spawn.h"

#if !defined(WIN32)
  #include <sys/stat.h>
  #include <sys/types.h>
  #ifdef HAVE_UNISTD_H
    #include <unistd.h>
  #endif
#endif

#include <algorithm>
#include <map>
#include <vector>
#include <utility>

/*! @cond Doxygen_Suppress */

CPL_CVSID("$Id: vrtderivedrasterband.cpp b2155cc20097aa9ab786396ed9bfe45686193e93 2019-06-20 10:43:54 +0200 Even Rouault $")

// #define GDAL_VRT_DISABLE_PYTHON
// #define PYTHONSO_DEFAULT "libpython2.7.so"

#ifndef GDAL_VRT_ENABLE_PYTHON_DEFAULT
// Can be YES, NO or TRUSTED_MODULES
#define GDAL_VRT_ENABLE_PYTHON_DEFAULT "TRUSTED_MODULES"
#endif

static std::map<CPLString, GDALDerivedPixelFunc> osMapPixelFunction;
static bool gbHasInitializedPython = false;
static int gnPythonInstanceCounter = 0;
static CPLMutex* ghMutex = nullptr;

// Subset of Python API defined as function of pointers
typedef struct PyObject_t PyObject;
#define Py_file_input 257
static void (*Py_SetProgramName)(const char*) = nullptr;
static PyObject* (*PyBuffer_FromReadWriteMemory)(void*, size_t) = nullptr;
static PyObject* (*PyTuple_New)(size_t) = nullptr;
static PyObject* (*PyInt_FromLong)(long) = nullptr;
static PyObject* (*PyFloat_FromDouble)(double) = nullptr;
static PyObject* (*PyObject_Call)(PyObject*, PyObject*, PyObject*) = nullptr;
static void (*Py_IncRef)(PyObject*) = nullptr;
static void (*Py_DecRef)(PyObject*) = nullptr;
static PyObject* (*PyErr_Occurred)(void) = nullptr;
static void (*PyErr_Print)(void) = nullptr;
static int (*Py_IsInitialized)(void) = nullptr;
static void (*Py_InitializeEx)(int) = nullptr;
static void (*PyEval_InitThreads)(void) = nullptr;
typedef struct PyThreadState_t PyThreadState;
static PyThreadState* (*PyEval_SaveThread)(void) = nullptr;
static void (*PyEval_RestoreThread)(PyThreadState*) = nullptr;
static void (*Py_Finalize)(void) = nullptr;
static PyObject* (*Py_CompileString)(const char*, const char*, int) = nullptr;
static PyObject* (*PyImport_ExecCodeModule)(const char*, PyObject*) = nullptr;
static PyObject* (*PyObject_GetAttrString)(PyObject*, const char*) = nullptr;
static int (*PyTuple_SetItem)(PyObject *, size_t, PyObject *) = nullptr;
static void (*PyObject_Print)(PyObject*,FILE*,int) = nullptr;
static PyObject* (*PyString_FromStringAndSize)(const void*, size_t) = nullptr;
static PyObject* (*PyImport_ImportModule)(const char*) = nullptr;
static int (*PyCallable_Check)(PyObject*) = nullptr;
static PyObject* (*PyDict_New)(void) = nullptr;
static int (*PyDict_SetItemString)(PyObject *p, const char *key,
                                   PyObject *val) = nullptr;
static void (*PyErr_Fetch)(PyObject **poPyType, PyObject **poPyValue,
                           PyObject **poPyTraceback) = nullptr;
static void (*PyErr_Clear)(void) = nullptr;
static const char* (*PyString_AsString)(PyObject*) = nullptr;
static const char* (*Py_GetVersion)(void) = nullptr;

typedef int PyGILState_STATE;
static PyGILState_STATE (*PyGILState_Ensure)(void) = nullptr;
static void (*PyGILState_Release)(PyGILState_STATE) = nullptr;

/* Flags for getting buffers */
#define PyBUF_WRITABLE 0x0001
#define PyBUF_FORMAT 0x0004
#define PyBUF_ND 0x0008
#define PyBUF_STRIDES (0x0010 | PyBUF_ND)
#define PyBUF_INDIRECT (0x0100 | PyBUF_STRIDES)
#define PyBUF_FULL (PyBUF_INDIRECT | PyBUF_WRITABLE | PyBUF_FORMAT)

typedef struct
{
    //cppcheck-suppress unusedStructMember
    char big_enough[256];
} Py_buffer;
static int (*PyBuffer_FillInfo)(Py_buffer *view, PyObject *obj, void *buf,
                                size_t len, int readonly, int infoflags) = nullptr;
static PyObject* (*PyMemoryView_FromBuffer)(Py_buffer *view) = nullptr;

static PyThreadState* gphThreadState = nullptr;

/************************************************************************/
/*                        GDALCreateNumpyArray()                        */
/************************************************************************/

static PyObject* GDALCreateNumpyArray(PyObject* pCreateArray,
                                      void* pBuffer,
                                      GDALDataType eType,
                                      int nHeight,
                                      int nWidth )
{
    PyObject* poPyBuffer;
    const size_t nSize = static_cast<size_t>(nHeight) * nWidth *
                                    GDALGetDataTypeSizeBytes(eType);
    if( PyBuffer_FromReadWriteMemory )
    {
        // Python 2
        poPyBuffer = PyBuffer_FromReadWriteMemory(pBuffer, nSize);
    }
    else
    {
        // Python 3
        Py_buffer pybuffer;
        if( PyBuffer_FillInfo(&pybuffer, nullptr, (char*)pBuffer,
                              nSize,
                              0, PyBUF_FULL) != 0)
        {
            return nullptr;
        }
        poPyBuffer = PyMemoryView_FromBuffer(&pybuffer);
    }
    PyObject* pArgsCreateArray = PyTuple_New(4);
    PyTuple_SetItem(pArgsCreateArray, 0, poPyBuffer);
    const char* pszDataType = nullptr;
    switch( eType )
    {
        case GDT_Byte: pszDataType = "uint8"; break;
        case GDT_UInt16: pszDataType = "uint16"; break;
        case GDT_Int16: pszDataType = "int16"; break;
        case GDT_UInt32: pszDataType = "uint32"; break;
        case GDT_Int32: pszDataType = "int32"; break;
        case GDT_Float32: pszDataType = "float32"; break;
        case GDT_Float64: pszDataType = "float64"; break;
        case GDT_CInt16:
        case GDT_CInt32:
            CPLAssert(FALSE);
            break;
        case GDT_CFloat32: pszDataType = "complex64"; break;
        case GDT_CFloat64: pszDataType = "complex128"; break;
        default:
            CPLAssert(FALSE);
            break;
    }
    PyTuple_SetItem(pArgsCreateArray, 1,
                PyString_FromStringAndSize(pszDataType, strlen(pszDataType)));
    PyTuple_SetItem(pArgsCreateArray, 2, PyInt_FromLong(nHeight));
    PyTuple_SetItem(pArgsCreateArray, 3, PyInt_FromLong(nWidth));
    PyObject* poNumpyArray = PyObject_Call(pCreateArray, pArgsCreateArray, nullptr);
    Py_DecRef(pArgsCreateArray);
    if (PyErr_Occurred())
        PyErr_Print();
    return poNumpyArray;
}

/* MinGW32 might define HAVE_DLFCN_H, so skip the unix implementation */
#if defined(HAVE_DLFCN_H) && !defined(WIN32)

#include <dlfcn.h>

typedef void* LibraryHandle;

#define LOAD_NOCHECK_WITH_NAME(libHandle, x, name) \
    do { \
            void* ptr = dlsym(libHandle, name); \
            memcpy(&x, &ptr, sizeof(void*)); \
    } while(0)

#elif defined(WIN32)

#include <windows.h>
#include <psapi.h>

typedef HMODULE LibraryHandle;

#define LOAD_NOCHECK_WITH_NAME(libHandle, x, name) \
    do { \
            FARPROC ptr = GetProcAddress(libHandle, name); \
            memcpy(&x, &ptr, sizeof(void*)); \
    } while(0)

#endif

#define STRINGIFY(x) #x

#define LOAD_NOCHECK(libHandle, x) LOAD_NOCHECK_WITH_NAME(libHandle, x, STRINGIFY(x))
#define LOAD_WITH_NAME(libHandle, x, name) \
    do { \
            LOAD_NOCHECK_WITH_NAME(libHandle, x, name); \
            if (!x) \
            { \
                CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s", name); \
                return false; \
            } \
    } while(0)
#define LOAD(libHandle, x) LOAD_WITH_NAME(libHandle, x, STRINGIFY(x))

/************************************************************************/
/*                          LoadPythonAPI()                             */
/************************************************************************/

#if defined(LOAD_NOCHECK_WITH_NAME) && defined(HAVE_DLFCN_H) && !defined(WIN32)
static LibraryHandle libHandleStatic = nullptr;
#endif

/** Load the subset of the Python C API that we need */
static bool LoadPythonAPI()
{
    CPLMutexHolder oHolder(&ghMutex);

    static bool bInit = false;
    if( bInit )
        return true;

#ifdef LOAD_NOCHECK_WITH_NAME
    // The static here is just to avoid Coverity warning about resource leak.
    LibraryHandle libHandle = nullptr;
    CPL_IGNORE_RET_VAL(libHandle);

    const char* pszPythonSO = CPLGetConfigOption("PYTHONSO", nullptr);
#if defined(HAVE_DLFCN_H) && !defined(WIN32)

    // First try in the current process in case the python symbols would
    // be already loaded
    libHandle = dlopen(nullptr, RTLD_LAZY);
    libHandleStatic = libHandle;
    if( libHandle != nullptr &&
        dlsym(libHandle, "Py_SetProgramName") != nullptr )
    {
        CPLDebug("VRT", "Current process has python symbols loaded");
    }
    else
    {
        libHandle = nullptr;
    }

    // Then try the user provided shared object name
    if( libHandle == nullptr && pszPythonSO != nullptr )
    {
        // coverity[tainted_string]
        libHandle = dlopen(pszPythonSO, RTLD_NOW | RTLD_GLOBAL);
        if( libHandle == nullptr )
        {
            CPLError(CE_Failure, CPLE_AppDefined,
                     "Cannot load %s",
                     pszPythonSO);
            return false;
        }
        if( dlsym(libHandle, "Py_SetProgramName") == nullptr )
        {
            CPLError(CE_Failure, CPLE_AppDefined,
                     "Cannot find Py_SetProgramName symbol in %s",
                     pszPythonSO);
            return false;
        }
    }

    // Then try the PYTHONSO_DEFAULT if defined at compile time
#ifdef PYTHONSO_DEFAULT
    if( libHandle == NULL )
    {
        libHandle = dlopen(PYTHONSO_DEFAULT, RTLD_NOW | RTLD_GLOBAL);
        if( !libHandle )
        {
            CPLDebug("VRT", "%s found", PYTHONSO_DEFAULT);
        }
    }
#endif

#if defined(__MACH__) && defined(__APPLE__)
#define SO_EXT "dylib"
#else
#define IS_SO_EXT
#define SO_EXT "so"
#endif

    const auto tryDlopen = [](CPLString osPythonSO)
    {
        CPLDebug("VRT", "Trying %s", osPythonSO.c_str());
        auto l_libHandle = dlopen(osPythonSO.c_str(), RTLD_NOW | RTLD_GLOBAL);
#ifdef IS_SO_EXT
        if( l_libHandle == nullptr )
        {
            osPythonSO += ".1.0";
            CPLDebug("VRT", "Trying %s", osPythonSO.c_str());
            l_libHandle = dlopen(osPythonSO.c_str(), RTLD_NOW | RTLD_GLOBAL);
        }
#endif
        return l_libHandle;
    };

    // Then try to find the libpython that corresponds to the python binary
    // in the PATH
    if( libHandle == nullptr )
    {
        CPLString osVersion;
        char* pszPath = getenv("PATH");
        if( pszPath != nullptr
#ifdef DEBUG
           // For testing purposes
           && CPLTestBool( CPLGetConfigOption(
                                    "VRT_ENABLE_PYTHON_PATH", "YES") )
#endif
          )
        {
            char** papszTokens = CSLTokenizeString2(pszPath, ":", 0);
            for( int iTry = 0; iTry < 2; ++iTry )
            {
                for( char** papszIter = papszTokens;
                        papszIter != nullptr && *papszIter != nullptr;
                        ++papszIter )
                {
                    struct stat sStat;
                    CPLString osPythonBinary(
                        CPLFormFilename(*papszIter, "python", nullptr));
                    if( iTry == 1 )
                        osPythonBinary += "3";
                    if( lstat(osPythonBinary, &sStat) != 0 )
                        continue;

                    CPLDebug("VRT", "Found %s", osPythonBinary.c_str());

                    if( S_ISLNK(sStat.st_mode)
#ifdef DEBUG
                        // For testing purposes
                        && CPLTestBool( CPLGetConfigOption(
                                    "VRT_ENABLE_PYTHON_SYMLINK", "YES") )
#endif
                        )
                    {
                        // If this is a symlink, hopefully the resolved
                        // name will be like "python2.7"
                        const int nBufSize = 2048;
                        std::vector<char> oFilename(nBufSize);
                        char *szPointerFilename = &oFilename[0];
                        int nBytes = static_cast<int>(
                            readlink( osPythonBinary, szPointerFilename,
                                      nBufSize ) );
                        if (nBytes != -1)
                        {
                            szPointerFilename[std::min(nBytes,
                                                       nBufSize - 1)] = 0;
                            CPLString osFilename(
                                            CPLGetFilename(szPointerFilename));
                            CPLDebug("VRT", "Which is an alias to: %s",
                                     szPointerFilename);
                            if( STARTS_WITH(osFilename, "python") )
                            {
                                osVersion = osFilename.substr(strlen("python"));
                                CPLDebug("VRT",
                                         "Python version from binary name: %s",
                                         osVersion.c_str());
                            }
                        }
                        else
                        {
                            CPLDebug("VRT", "realink(%s) failed",
                                        osPythonBinary.c_str());
                        }
                    }

                    // Otherwise, expensive way: start the binary and ask
                    // it for its version...
                    if( osVersion.empty() )
                    {
                        const char* pszPrintVersion =
                            "import sys; print(str(sys.version_info[0]) +"
                            "'.' + str(sys.version_info[1]))";
                        const char* const apszArgv[] = {
                                osPythonBinary.c_str(), "-c",
                                pszPrintVersion,
                                nullptr };
                        const CPLString osTmpFilename(
                                        "/vsimem/LoadPythonAPI/out.txt");
                        VSILFILE* fout = VSIFOpenL( osTmpFilename, "wb+");
                        if( CPLSpawn( apszArgv, nullptr, fout, FALSE ) == 0 )
                        {
                            char* pszStr = reinterpret_cast<char*>(
                                VSIGetMemFileBuffer( osTmpFilename,
                                                        nullptr, FALSE ));
                            osVersion = pszStr;
                            if( !osVersion.empty() &&
                                osVersion.back() == '\n' )
                            {
                                osVersion.resize(osVersion.size() - 1);
                            }
                            CPLDebug("VRT", "Python version from binary: %s",
                                        osVersion.c_str());
                        }
                        VSIFCloseL(fout);
                        VSIUnlink(osTmpFilename);
                    }
                    break;
                }
                if( !osVersion.empty() )
                    break;
            }
            CSLDestroy(papszTokens);
        }

        if( !osVersion.empty() )
        {
            libHandle = tryDlopen("libpython" + osVersion + "." SO_EXT);
            if( libHandle != nullptr )
            {
                CPLDebug("VRT", "... success");
            }
            else if( osVersion[0] == '3' )
            {
                libHandle = tryDlopen("libpython" + osVersion + "m." SO_EXT);
                if( libHandle != nullptr )
                {
                    CPLDebug("VRT", "... success");
                }
            }
        }
    }

    // Otherwise probe a few known objects.
    // Note: update doc/source/drivers/raster/vrt.rst if change
    if( libHandle == nullptr )
    {
        const char* const apszPythonSO[] = { "libpython2.7." SO_EXT,
                                                "libpython2.6." SO_EXT,
                                                "libpython3.4m." SO_EXT,
                                                "libpython3.5m." SO_EXT,
                                                "libpython3.6m." SO_EXT,
                                                "libpython3.3." SO_EXT,
                                                "libpython3.2." SO_EXT };
        for( size_t i = 0; libHandle == nullptr &&
                            i < CPL_ARRAYSIZE(apszPythonSO); ++i )
        {
            libHandle = tryDlopen(apszPythonSO[i]);
            if( libHandle != nullptr )
                CPLDebug("VRT", "... success");
        }
    }

#elif defined(WIN32)

    // First try in the current process in case the python symbols would
    // be already loaded
    HANDLE hProcess = GetCurrentProcess();
    const size_t nMaxModules = 10000;
    HMODULE* pahModules = static_cast<HMODULE*>(
        CPLMalloc(nMaxModules * sizeof(HMODULE)));
    DWORD nSizeNeeded = 0;

    EnumProcessModules(hProcess, pahModules, nMaxModules * sizeof(HMODULE),
                        &nSizeNeeded);

    const size_t nModules =
        std::min(nMaxModules,
                 static_cast<size_t>(nSizeNeeded) / sizeof(HMODULE));
    for( size_t i = 0; i < nModules; i++ )
    {
        if( GetProcAddress(pahModules[i], "Py_SetProgramName") )
        {
            libHandle = pahModules[i];
            CPLDebug("VRT", "Current process has python symbols loaded");
            break;
        }
    }
    CPLFree(pahModules);

    // Then try the user provided shared object name
    if( libHandle == nullptr && pszPythonSO != nullptr )
    {
        UINT        uOldErrorMode;
        /* Avoid error boxes to pop up (#5211, #5525) */
        uOldErrorMode = SetErrorMode(SEM_NOOPENFILEERRORBOX |
                                     SEM_FAILCRITICALERRORS);

#if (defined(WIN32) && _MSC_VER >= 1310) || __MSVCRT_VERSION__ >= 0x0601
        if( CPLTestBool( CPLGetConfigOption( "GDAL_FILENAME_IS_UTF8", "YES" ) ) )
        {
            wchar_t *pwszFilename =
                CPLRecodeToWChar( pszPythonSO, CPL_ENC_UTF8, CPL_ENC_UCS2 );
            libHandle = LoadLibraryW(pwszFilename);
            CPLFree( pwszFilename );
        }
        else
#endif
        {
            libHandle = LoadLibrary(pszPythonSO);
        }

        SetErrorMode(uOldErrorMode);

        if( libHandle == nullptr )
        {
            CPLError(CE_Failure, CPLE_AppDefined,
                     "Cannot load %s",
                     pszPythonSO);
            return false;
        }
        if( GetProcAddress(libHandle, "Py_SetProgramName") == nullptr )
        {
            CPLError(CE_Failure, CPLE_AppDefined,
                     "Cannot find Py_SetProgramName symbol in %s",
                     pszPythonSO);
            return false;
        }
    }

    // Then try the PYTHONSO_DEFAULT if defined at compile time
#ifdef PYTHONSO_DEFAULT
    if( libHandle == nullptr )
    {
        UINT        uOldErrorMode;
        uOldErrorMode = SetErrorMode(SEM_NOOPENFILEERRORBOX |
                                        SEM_FAILCRITICALERRORS);

        libHandle = LoadLibrary(PYTHONSO_DEFAULT);
        SetErrorMode(uOldErrorMode);
        if( !libHandle )
        {
            CPLDebug("VRT", "%s found", PYTHONSO_DEFAULT);
        }
    }
#endif

    // Then try to find the pythonXY.dll that corresponds to the python binary
    // in the PATH
    if( libHandle == nullptr )
    {
        CPLString osDLLName;
        char* pszPath = getenv("PATH");
        if( pszPath != nullptr
#ifdef DEBUG
           // For testing purposes
           && CPLTestBool( CPLGetConfigOption(
                                    "VRT_ENABLE_PYTHON_PATH", "YES") )
#endif
          )
        {
            char** papszTokens = CSLTokenizeString2(pszPath, ";", 0);
            for( int iTry = 0; iTry < 2; ++iTry )
            {
                for( char** papszIter = papszTokens;
                        papszIter != nullptr && *papszIter != nullptr;
                        ++papszIter )
                {
                    VSIStatBufL sStat;
                    CPLString osPythonBinary(
                            CPLFormFilename(*papszIter, "python.exe", nullptr));
                    if( iTry == 1 )
                        osPythonBinary += "3";
                    if( VSIStatL(osPythonBinary, &sStat) != 0 )
                        continue;

                    CPLDebug("VRT", "Found %s", osPythonBinary.c_str());

                    // In python2.7, the dll is in the same directory as the exe
                    char** papszFiles = VSIReadDir(*papszIter);
                    for( char** papszFileIter = papszFiles;
                                papszFileIter != nullptr && *papszFileIter != nullptr;
                                ++papszFileIter )
                    {
                        if( STARTS_WITH_CI(*papszFileIter, "python") &&
                            EQUAL(CPLGetExtension(*papszFileIter), "dll") )
                        {
                            osDLLName = CPLFormFilename(*papszIter,
                                                        *papszFileIter,
                                                        nullptr);
                            break;
                        }
                    }
                    CSLDestroy(papszFiles);

                    // In python3.2, the dll is in the DLLs subdirectory
                    if( osDLLName.empty() )
                    {
                        CPLString osDLLsDir(
                                CPLFormFilename(*papszIter, "DLLs", nullptr));
                        papszFiles = VSIReadDir( osDLLsDir );
                        for( char** papszFileIter = papszFiles;
                                    papszFileIter != nullptr && *papszFileIter != nullptr;
                                    ++papszFileIter )
                        {
                            if( STARTS_WITH_CI(*papszFileIter, "python") &&
                                EQUAL(CPLGetExtension(*papszFileIter), "dll") )
                            {
                                osDLLName = CPLFormFilename(osDLLsDir,
                                                            *papszFileIter,
                                                            nullptr);
                                break;
                            }
                        }
                        CSLDestroy(papszFiles);
                    }

                    break;
                }
                if( !osDLLName.empty() )
                    break;
            }
            CSLDestroy(papszTokens);
        }

        if( !osDLLName.empty() )
        {
            //CPLDebug("VRT", "Trying %s", osDLLName.c_str());
            UINT        uOldErrorMode;
            uOldErrorMode = SetErrorMode(SEM_NOOPENFILEERRORBOX |
                                            SEM_FAILCRITICALERRORS);
            libHandle = LoadLibrary(osDLLName);
            SetErrorMode(uOldErrorMode);
            if( libHandle != nullptr )
            {
                CPLDebug("VRT", "%s loaded", osDLLName.c_str());
            }
        }
    }

    // Otherwise probe a few known objects
    // Note: update vrt_tutorial.dox if change
    if( libHandle == nullptr )
    {
        const char* const apszPythonSO[] = { "python27.dll",
                                            "python26.dll",
                                            "python34.dll",
                                            "python35.dll",
                                            "python36.dll",
                                            "python33.dll",
                                            "python32.dll" };
        UINT        uOldErrorMode;
        uOldErrorMode = SetErrorMode(SEM_NOOPENFILEERRORBOX |
                                        SEM_FAILCRITICALERRORS);

        for( size_t i = 0; libHandle == nullptr &&
                            i < CPL_ARRAYSIZE(apszPythonSO); ++i )
        {
            CPLDebug("VRT", "Trying %s", apszPythonSO[i]);
            libHandle = LoadLibrary(apszPythonSO[i]);
            if( libHandle != nullptr )
                CPLDebug("VRT", "... success");
        }
        SetErrorMode(uOldErrorMode);
    }
#endif
    if( !libHandle )
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Cannot find python/libpython. You can set the PYTHONSO "
                 "configuration option to point to the a python .so/.dll/.dylib");
        return false;
    }

    LOAD(libHandle, Py_SetProgramName);
    LOAD_NOCHECK(libHandle, PyBuffer_FromReadWriteMemory);
    LOAD_NOCHECK(libHandle, PyBuffer_FillInfo);
    LOAD_NOCHECK(libHandle, PyMemoryView_FromBuffer);
    if( PyBuffer_FromReadWriteMemory == nullptr &&
        (PyBuffer_FillInfo == nullptr || PyMemoryView_FromBuffer == nullptr) )
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Cannot find PyBuffer_FillInfo or "
                 "PyBuffer_FillInfo+PyMemoryView_FromBuffer\n");
        return false;
    }
    LOAD(libHandle, PyTuple_New);
    if( PyBuffer_FromReadWriteMemory )
    {
        // Python 2
        LOAD(libHandle, PyInt_FromLong);
        LOAD(libHandle, PyString_FromStringAndSize);
        LOAD(libHandle, PyString_AsString);
    }
    else
    {
        // Python 3
        LOAD_WITH_NAME(libHandle, PyInt_FromLong, "PyLong_FromLong");
        LOAD_WITH_NAME(libHandle, PyString_FromStringAndSize,
                                        "PyBytes_FromStringAndSize");
        LOAD_WITH_NAME(libHandle, PyString_AsString, "PyBytes_AsString");
    }
    LOAD(libHandle, PyFloat_FromDouble);
    LOAD(libHandle, PyObject_Call);
    LOAD(libHandle, Py_IncRef);
    LOAD(libHandle, Py_DecRef);
    LOAD(libHandle, PyErr_Occurred);
    LOAD(libHandle, PyErr_Print);
    LOAD(libHandle, Py_IsInitialized);
    LOAD(libHandle, Py_InitializeEx);
    LOAD(libHandle, PyEval_InitThreads);
    LOAD(libHandle, PyEval_SaveThread);
    LOAD(libHandle, PyEval_RestoreThread);
    LOAD(libHandle, Py_Finalize);
    LOAD(libHandle, Py_CompileString);
    LOAD(libHandle, PyImport_ExecCodeModule);
    LOAD(libHandle, PyObject_GetAttrString);
    LOAD(libHandle, PyTuple_SetItem);
    LOAD(libHandle, PyObject_Print);
    LOAD(libHandle, PyImport_ImportModule);
    LOAD(libHandle, PyCallable_Check);
    LOAD(libHandle, PyDict_New);
    LOAD(libHandle, PyDict_SetItemString);
    LOAD(libHandle, PyGILState_Ensure);
    LOAD(libHandle, PyGILState_Release);
    LOAD(libHandle, PyErr_Fetch);
    LOAD(libHandle, PyErr_Clear);
    LOAD(libHandle, Py_GetVersion);

    CPLString osPythonVersion(Py_GetVersion());
    osPythonVersion.replaceAll("\r\n", ' ');
    osPythonVersion.replaceAll('\n', ' ');
    CPLDebug("VRT", "Python version used: %s", osPythonVersion.c_str());

#else // LOAD_NOCHECK_WITH_NAME
    CPLError(CE_Failure, CPLE_AppDefined,
             "This platform doesn't support dynamic loading of libraries")
    return false;
#endif // LOAD_NOCHECK_WITH_NAME

    bInit = true;
    return bInit;
}

/************************************************************************/
/*                      GetPyExceptionString()                          */
/************************************************************************/

static CPLString GetPyExceptionString()
{
    PyObject *poPyType = nullptr;
    PyObject *poPyValue = nullptr;
    PyObject *poPyTraceback = nullptr;

    PyErr_Fetch(&poPyType, &poPyValue, &poPyTraceback);
    if( poPyType )
        Py_IncRef(poPyType);
    if( poPyValue )
        Py_IncRef(poPyValue);
    if( poPyTraceback )
        Py_IncRef(poPyTraceback);

    // This is a mess. traceback.format_exception/format_exception_only
    // sometimes throw exceptions themselves !
    CPLString osPythonCode(
        "import traceback\n"
        "\n"
        "def GDALFormatException2(etype, value):\n"
        "    try:\n"
        "       return ''.join(traceback.format_exception_only(etype, value)).encode('UTF-8')\n"
        "    except:\n"
        "       return (str(etype) + ', ' + str(value)).encode('UTF-8')\n"
        "\n"
        "def GDALFormatException3(etype, value, tb):\n"
        //"    print(etype, value, tb)\n"
        "    try:\n"
        "       return ''.join(traceback.format_exception(etype, value, tb)).encode('UTF-8')\n"
        "    except:\n"
        "       return (str(etype) + ', ' + str(value)).encode('UTF-8')\n");

    CPLString osRet("An exception occurred in exception formatting code...");

    static int nCounter = 0;
    CPLString osModuleName( CPLSPrintf("gdal_exception_%d", nCounter));
    PyObject* poCompiledString = Py_CompileString(osPythonCode,
                                                  osModuleName, Py_file_input);
    if( poCompiledString == nullptr || PyErr_Occurred() )
    {
        PyErr_Print();
    }
    else
    {
        PyObject* poModule =
            PyImport_ExecCodeModule(osModuleName, poCompiledString);
        CPLAssert(poModule);

        Py_DecRef(poCompiledString);

        PyObject* poPyGDALFormatException2 = PyObject_GetAttrString(poModule,
                                                "GDALFormatException2" );
        CPLAssert(poPyGDALFormatException2);

        PyObject* poPyGDALFormatException3 = PyObject_GetAttrString(poModule,
                                                "GDALFormatException3" );
        CPLAssert(poPyGDALFormatException3);

        Py_DecRef(poModule);

        PyObject* pyArgs = PyTuple_New( poPyTraceback ? 3 : 2);
        PyTuple_SetItem(pyArgs, 0, poPyType);
        PyTuple_SetItem(pyArgs, 1, poPyValue);
        if( poPyTraceback )
            PyTuple_SetItem(pyArgs, 2, poPyTraceback );
        PyObject* poPyRet = PyObject_Call(
            poPyTraceback ? poPyGDALFormatException3 : poPyGDALFormatException2,
            pyArgs, nullptr );
        Py_DecRef(pyArgs);

        if( PyErr_Occurred() )
        {
            osRet = "An exception occurred in exception formatting code...";
            PyErr_Print();
        }
        else
        {
            osRet = PyString_AsString(poPyRet);
            Py_DecRef(poPyRet);
        }

        Py_DecRef(poPyGDALFormatException2);
        Py_DecRef(poPyGDALFormatException3);
    }

    if( poPyType )
        Py_DecRef(poPyType);
    if( poPyValue )
        Py_DecRef(poPyValue);
    if( poPyTraceback )
        Py_DecRef(poPyTraceback);

    return osRet;
}

/************************************************************************/
/* ==================================================================== */
/*                     VRTDerivedRasterBandPrivateData                  */
/* ==================================================================== */
/************************************************************************/

class VRTDerivedRasterBandPrivateData
{
    public:
        CPLString m_osCode;
        CPLString m_osLanguage;
        int       m_nBufferRadius;
        PyObject* m_poGDALCreateNumpyArray;
        PyObject* m_poUserFunction;
        bool      m_bPythonInitializationDone;
        bool      m_bPythonInitializationSuccess;
        bool      m_bExclusiveLock;
        bool      m_bFirstTime;
        std::vector< std::pair<CPLString,CPLString> > m_oFunctionArgs;

        VRTDerivedRasterBandPrivateData():
            m_osLanguage("C"),
            m_nBufferRadius(0),
            m_poGDALCreateNumpyArray(nullptr),
            m_poUserFunction(nullptr),
            m_bPythonInitializationDone(false),
            m_bPythonInitializationSuccess(false),
            m_bExclusiveLock(false),
            m_bFirstTime(true)
        {
        }

        virtual ~VRTDerivedRasterBandPrivateData()
        {
            if( m_poGDALCreateNumpyArray )
                Py_DecRef(m_poGDALCreateNumpyArray);
            if( m_poUserFunction )
                Py_DecRef(m_poUserFunction);

            CPLMutexHolder oHolder(&ghMutex);
            gnPythonInstanceCounter --;
        }
};

/************************************************************************/
/* ==================================================================== */
/*                            VRT_GIL_Holder                            */
/* ==================================================================== */
/************************************************************************/

class VRT_GIL_Holder
{
        bool             m_bExclusiveLock;
        PyGILState_STATE m_eState;

    public:

        explicit VRT_GIL_Holder(bool bExclusiveLock);
        virtual ~VRT_GIL_Holder();
};

VRT_GIL_Holder::VRT_GIL_Holder(bool bExclusiveLock):
    m_bExclusiveLock(bExclusiveLock)
{
    if( bExclusiveLock )
    {
        if( ghMutex )
            CPLAcquireMutex( ghMutex, 1000.0 );
    }
    m_eState = PyGILState_Ensure();
}

VRT_GIL_Holder::~VRT_GIL_Holder()
{
    PyGILState_Release(m_eState);
    if( m_bExclusiveLock )
    {
        if( ghMutex )
            CPLReleaseMutex( ghMutex );
    }
    else
    {
    }
}

/************************************************************************/
/* ==================================================================== */
/*                          VRTDerivedRasterBand                        */
/* ==================================================================== */
/************************************************************************/

/************************************************************************/
/*                        VRTDerivedRasterBand()                        */
/************************************************************************/

VRTDerivedRasterBand::VRTDerivedRasterBand( GDALDataset *poDSIn, int nBandIn ) :
    VRTSourcedRasterBand( poDSIn, nBandIn ),
    m_poPrivate(nullptr),
    pszFuncName(nullptr),
    eSourceTransferType(GDT_Unknown)
{
    m_poPrivate = new VRTDerivedRasterBandPrivateData;
}

/************************************************************************/
/*                        VRTDerivedRasterBand()                        */
/************************************************************************/

VRTDerivedRasterBand::VRTDerivedRasterBand( GDALDataset *poDSIn, int nBandIn,
                                            GDALDataType eType,
                                            int nXSize, int nYSize ) :
    VRTSourcedRasterBand(poDSIn, nBandIn, eType, nXSize, nYSize),
    m_poPrivate(nullptr),
    pszFuncName(nullptr),
    eSourceTransferType(GDT_Unknown)
{
    m_poPrivate = new VRTDerivedRasterBandPrivateData;
}

/************************************************************************/
/*                       ~VRTDerivedRasterBand()                        */
/************************************************************************/

VRTDerivedRasterBand::~VRTDerivedRasterBand()

{
    CPLFree( pszFuncName );
    delete m_poPrivate;
}

/************************************************************************/
/*                               Cleanup()                              */
/************************************************************************/

void VRTDerivedRasterBand::Cleanup()
{
    if( ghMutex )
        CPLDestroyMutex(ghMutex);
    ghMutex = nullptr;

    if( gnPythonInstanceCounter == 0 && gbHasInitializedPython &&
        CPLTestBool(CPLGetConfigOption("GDAL_VRT_ENABLE_PYTHON_FINALIZE",
                                       "YES")) )
    {
        // We call Py_Finalize at driver destruction, rather at dataset
        // destruction, since numpy crashes when it is reloaded after the next
        // Py_Initialize
        CPLDebug("VRT", "Py_Finalize() = %p", Py_Finalize);
        PyEval_RestoreThread(gphThreadState);
        Py_Finalize();
        gbHasInitializedPython = false;
        gphThreadState = nullptr;
    }
}

/************************************************************************/
/*                           AddPixelFunction()                         */
/************************************************************************/

/*! @endcond */

/**
 * This adds a pixel function to the global list of available pixel
 * functions for derived bands.  Pixel functions must be registered
 * in this way before a derived band tries to access data.
 *
 * Derived bands are stored with only the name of the pixel function
 * that it will apply, and if a pixel function matching the name is not
 * found the IRasterIO() call will do nothing.
 *
 * @param pszFuncName Name used to access pixel function
 * @param pfnNewFunction Pixel function associated with name.  An
 *  existing pixel function registered with the same name will be
 *  replaced with the new one.
 *
 * @return CE_None, invalid (NULL) parameters are currently ignored.
 */
CPLErr CPL_STDCALL
GDALAddDerivedBandPixelFunc( const char *pszFuncName,
                             GDALDerivedPixelFunc pfnNewFunction )
{
    if( pszFuncName == nullptr || pszFuncName[0] == '\0' ||
        pfnNewFunction == nullptr )
    {
      return CE_None;
    }

    osMapPixelFunction[pszFuncName] = pfnNewFunction;

    return CE_None;
}

/*! @cond Doxygen_Suppress */

/**
 * This adds a pixel function to the global list of available pixel
 * functions for derived bands.
 *
 * This is the same as the c function GDALAddDerivedBandPixelFunc()
 *
 * @param pszFuncName Name used to access pixel function
 * @param pfnNewFunction Pixel function associated with name.  An
 *  existing pixel function registered with the same name will be
 *  replaced with the new one.
 *
 * @return CE_None, invalid (NULL) parameters are currently ignored.
 */
CPLErr
VRTDerivedRasterBand::AddPixelFunction(
    const char *pszFuncName, GDALDerivedPixelFunc pfnNewFunction )
{
    return GDALAddDerivedBandPixelFunc(pszFuncName, pfnNewFunction);
}

/************************************************************************/
/*                           GetPixelFunction()                         */
/************************************************************************/

/**
 * Get a pixel function previously registered using the global
 * AddPixelFunction.
 *
 * @param pszFuncName The name associated with the pixel function.
 *
 * @return A derived band pixel function, or NULL if none have been
 * registered for pszFuncName.
 */
GDALDerivedPixelFunc
VRTDerivedRasterBand::GetPixelFunction( const char *pszFuncName )
{
    if( pszFuncName == nullptr || pszFuncName[0] == '\0' )
    {
        return nullptr;
    }

    std::map<CPLString, GDALDerivedPixelFunc>::iterator oIter =
        osMapPixelFunction.find(pszFuncName);

    if( oIter == osMapPixelFunction.end())
        return nullptr;

    return oIter->second;
}

/************************************************************************/
/*                         SetPixelFunctionName()                       */
/************************************************************************/

/**
 * Set the pixel function name to be applied to this derived band.  The
 * name should match a pixel function registered using AddPixelFunction.
 *
 * @param pszFuncNameIn Name of pixel function to be applied to this derived
 * band.
 */
void VRTDerivedRasterBand::SetPixelFunctionName( const char *pszFuncNameIn )
{
    pszFuncName = CPLStrdup( pszFuncNameIn );
}

/************************************************************************/
/*                         SetPixelFunctionLanguage()                   */
/************************************************************************/

/**
 * Set the language of the pixel function.
 *
 * @param pszLanguage Language of the pixel function (only "C" and "Python"
 * are supported currently)
 * @since GDAL 2.3
 */
void VRTDerivedRasterBand::SetPixelFunctionLanguage( const char* pszLanguage )
{
    m_poPrivate->m_osLanguage = pszLanguage;
}

/************************************************************************/
/*                         SetSourceTransferType()                      */
/************************************************************************/

/**
 * Set the transfer type to be used to obtain pixel information from
 * all of the sources.  If unset, the transfer type used will be the
 * same as the derived band data type.  This makes it possible, for
 * example, to pass CFloat32 source pixels to the pixel function, even
 * if the pixel function generates a raster for a derived band that
 * is of type Byte.
 *
 * @param eDataTypeIn Data type to use to obtain pixel information from
 * the sources to be passed to the derived band pixel function.
 */
void VRTDerivedRasterBand::SetSourceTransferType( GDALDataType eDataTypeIn )
{
    eSourceTransferType = eDataTypeIn;
}

/************************************************************************/
/*                           InitializePython()                         */
/************************************************************************/

bool VRTDerivedRasterBand::InitializePython()
{
    if( m_poPrivate->m_bPythonInitializationDone )
        return m_poPrivate->m_bPythonInitializationSuccess;

    m_poPrivate->m_bPythonInitializationDone = true;
    m_poPrivate->m_bPythonInitializationSuccess = false;

    const CPLString osPythonFullname( pszFuncName ? pszFuncName : "" );
    const size_t nIdxDot = osPythonFullname.rfind(".");
    CPLString osPythonModule;
    CPLString osPythonFunction;
    if( nIdxDot != std::string::npos )
    {
        osPythonModule = osPythonFullname.substr(0, nIdxDot);
        osPythonFunction = osPythonFullname.substr(nIdxDot+1);
    }
    else
    {
        osPythonFunction = osPythonFullname;
    }

#ifndef GDAL_VRT_DISABLE_PYTHON
    const char* pszPythonEnabled =
                            CPLGetConfigOption("GDAL_VRT_ENABLE_PYTHON", nullptr);
#else
    const char* pszPythonEnabled = "NO";
#endif
    const CPLString osPythonEnabled(pszPythonEnabled ? pszPythonEnabled :
                                            GDAL_VRT_ENABLE_PYTHON_DEFAULT);

    if( EQUAL(osPythonEnabled, "TRUSTED_MODULES") )
    {
        bool bIsTrustedModule = false;
        const CPLString osVRTTrustedModules(
                    CPLGetConfigOption( "GDAL_VRT_PYTHON_TRUSTED_MODULES", "") );
        if( !osPythonModule.empty() )
        {
            char** papszTrustedModules = CSLTokenizeString2(
                                                osVRTTrustedModules, ",", 0 );
            for( char** papszIter = papszTrustedModules;
                !bIsTrustedModule && papszIter && *papszIter;
                ++papszIter )
            {
                const char* pszIterModule = *papszIter;
                size_t nIterModuleLen = strlen(pszIterModule);
                if( nIterModuleLen > 2 &&
                    strncmp(pszIterModule + nIterModuleLen - 2, ".*", 2) == 0 )
                {
                    bIsTrustedModule =
                        (strncmp( osPythonModule, pszIterModule,
                                                  nIterModuleLen - 2 ) == 0) &&
                        (osPythonModule.size() == nIterModuleLen - 2 ||
                         (osPythonModule.size() >= nIterModuleLen &&
                          osPythonModule[nIterModuleLen-1] == '.') );
                }
                else if( nIterModuleLen >= 1 &&
                        pszIterModule[nIterModuleLen-1] == '*' )
                {
                    bIsTrustedModule = (strncmp( osPythonModule, pszIterModule,
                                                nIterModuleLen - 1 ) == 0);
                }
                else
                {
                    bIsTrustedModule =
                                (strcmp(osPythonModule, pszIterModule) == 0);
                }
            }
            CSLDestroy(papszTrustedModules);
        }

        if( !bIsTrustedModule )
        {
            if( osPythonModule.empty() )
            {
                CPLError(CE_Failure, CPLE_AppDefined,
                         "Python code needs to be executed, but it uses online code "
                         "in the VRT whereas the current policy is to trust only "
                         "code from external trusted modules (defined in the "
                         "GDAL_VRT_PYTHON_TRUSTED_MODULES configuration option). "
                         "If you trust the code in %s, you can set the "
                         "GDAL_VRT_ENABLE_PYTHON configuration option to YES.",
                         GetDataset() ? GetDataset()->GetDescription() :
                                    "(unknown VRT)");
            }
            else if( osVRTTrustedModules.empty() )
            {
                CPLError(CE_Failure, CPLE_AppDefined,
                         "Python code needs to be executed, but it uses code "
                         "from module '%s', whereas the current policy is to "
                         "trust only code from modules defined in the "
                         "GDAL_VRT_PYTHON_TRUSTED_MODULES configuration option, "
                         "which is currently unset. "
                         "If you trust the code in '%s', you can add module '%s' "
                         "to GDAL_VRT_PYTHON_TRUSTED_MODULES (or set the "
                         "GDAL_VRT_ENABLE_PYTHON configuration option to YES).",
                         osPythonModule.c_str(),
                         GetDataset() ? GetDataset()->GetDescription() :
                                    "(unknown VRT)",
                         osPythonModule.c_str());
            }
            else
            {
                CPLError(CE_Failure, CPLE_AppDefined,
                         "Python code needs to be executed, but it uses code "
                         "from module '%s', whereas the current policy is to "
                         "trust only code from modules '%s' (defined in the "
                         "GDAL_VRT_PYTHON_TRUSTED_MODULES configuration option). "
                         "If you trust the code in '%s', you can add module '%s' "
                         "to GDAL_VRT_PYTHON_TRUSTED_MODULES (or set the "
                         "GDAL_VRT_ENABLE_PYTHON configuration option to YES).",
                         osPythonModule.c_str(),
                         osVRTTrustedModules.c_str(),
                         GetDataset() ? GetDataset()->GetDescription() :
                                    "(unknown VRT)",
                         osPythonModule.c_str());
            }
            return false;
        }
    }

#ifdef disabled_because_this_is_probably_broken_by_design
    // See https://lwn.net/Articles/574215/
    // and http://nedbatchelder.com/blog/201206/eval_really_is_dangerous.html
    else if( EQUAL(osPythonEnabled, "IF_SAFE") )
    {
        bool bSafe = true;
        // If the function comes from another module, then we don't know
        if( !osPythonModule.empty() )
        {
            CPLDebug("VRT", "Python function is from another module");
            bSafe = false;
        }

        CPLString osCode(m_poPrivate->m_osCode);

        // Reject all imports except a few trusted modules
        const char* const apszTrustedImports[] = {
                "import math",
                "from math import",
                "import numpy", // caution: numpy has lots of I/O functions !
                "from numpy import",
                // TODO: not sure if importing arbitrary stuff from numba is OK
                // so let's just restrict to jit.
                "from numba import jit",

                // Not imports but still whitelisted, whereas other __ is banned
                "__init__",
                "__call__",
        };
        for( size_t i = 0; i < CPL_ARRAYSIZE(apszTrustedImports); ++i )
        {
            osCode.replaceAll(CPLString(apszTrustedImports[i]), "");
        }

        // Some dangerous built-in functions or numpy functions
        const char* const apszUntrusted[] = { "import", // and __import__
                                              "eval",
                                              "compile",
                                              "open",
                                              "load", // reload, numpy.load
                                              "file", // and exec_file, numpy.fromfile, numpy.tofile
                                              "input", // and raw_input
                                              "save", // numpy.save
                                              "memmap", // numpy.memmap
                                              "DataSource", // numpy.DataSource
                                              "genfromtxt", // numpy.genfromtxt
                                              "getattr",
                                              "ctypeslib", // numpy.ctypeslib
                                              "testing", // numpy.testing
                                              "dump", // numpy.ndarray.dump
                                              "fromregex", // numpy.fromregex
                                              "__"
                                             };
        for( size_t i = 0; i < CPL_ARRAYSIZE(apszUntrusted); ++i )
        {
            if( osCode.find(apszUntrusted[i]) != std::string::npos )
            {
                CPLDebug("VRT", "Found '%s' word in Python code",
                         apszUntrusted[i]);
                bSafe = false;
            }
        }

        if( !bSafe )
        {
            CPLError(CE_Failure, CPLE_AppDefined,
                     "Python code needs to be executed, but we cannot verify "
                     "if it is safe, so this is disabled by default. "
                     "If you trust the code in %s, you can set the "
                     "GDAL_VRT_ENABLE_PYTHON configuration option to YES.",
                     GetDataset() ? GetDataset()->GetDescription() :
                                    "(unknown VRT)");
            return false;
        }
    }
#endif //disabled_because_this_is_probably_broken_by_design

    else if( !EQUAL(osPythonEnabled, "YES") &&
             !EQUAL(osPythonEnabled, "ON") &&
             !EQUAL(osPythonEnabled, "TRUE") )
    {
        if( pszPythonEnabled == nullptr )
        {
            // Note: this is dead code with our current default policy
            // GDAL_VRT_ENABLE_PYTHON == "TRUSTED_MODULES"
            CPLError(CE_Failure, CPLE_AppDefined,
                 "Python code needs to be executed, but this is "
                 "disabled by default. If you trust the code in %s, "
                 "you can set the GDAL_VRT_ENABLE_PYTHON configuration "
                 "option to YES.",
                GetDataset() ? GetDataset()->GetDescription() :
                                                    "(unknown VRT)");
        }
        else
        {
            CPLError(CE_Failure, CPLE_AppDefined,
                    "Python code in %s needs to be executed, but this has been "
                    "explicitly disabled.",
                     GetDataset() ? GetDataset()->GetDescription() :
                                                            "(unknown VRT)");
        }
        return false;
    }

    if( !LoadPythonAPI() )
        return false;

    // Make sure the python interpreter is initialized
    {
        CPLMutexHolder oHolder(&ghMutex);
        int bIsInitialized = Py_IsInitialized();
        if( !bIsInitialized)
        {
            gbHasInitializedPython = true;
            Py_InitializeEx(0);
            CPLDebug("VRT", "Py_Initialize()");
            PyEval_InitThreads();
            gphThreadState = PyEval_SaveThread();
        }
        gnPythonInstanceCounter ++;
    }

    // Whether we should just use our own global mutex, in addition to Python
    // GIL locking.
    m_poPrivate->m_bExclusiveLock =
        CPLTestBool(CPLGetConfigOption("GDAL_VRT_PYTHON_EXCLUSIVE_LOCK", "NO"));

    // numba jit'ification doesn't seem to be thread-safe, so force use of
    // lock now and at first execution of function. Later executions seem to
    // be thread-safe. This problem doesn't seem to appear for code in
    // regular files
    const bool bUseExclusiveLock = m_poPrivate->m_bExclusiveLock ||
                    m_poPrivate->m_osCode.find("@jit") != std::string::npos;
    VRT_GIL_Holder oHolder(bUseExclusiveLock);

    // As we don't want to depend on numpy C API/ABI, we use a trick to build
    // a numpy array object. We define a Python function to which we pass a
    // Python buffer object.

    // We need to build a unique module name, otherwise this will crash in
    // multithreaded use cases.
    CPLString osModuleName( CPLSPrintf("gdal_vrt_module_%p", this) );
    PyObject* poCompiledString = Py_CompileString(
        ("import numpy\n"
        "def GDALCreateNumpyArray(buffer, dtype, height, width):\n"
        "    return numpy.frombuffer(buffer, str(dtype.decode('ascii')))."
                                                "reshape([height, width])\n"
        "\n" + m_poPrivate->m_osCode).c_str(),
        osModuleName, Py_file_input);
    if( poCompiledString == nullptr || PyErr_Occurred() )
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Couldn't compile code:\n%s",
                 GetPyExceptionString().c_str());
        return false;
    }
    PyObject* poModule =
        PyImport_ExecCodeModule(osModuleName, poCompiledString);
    Py_DecRef(poCompiledString);

    if( poModule == nullptr || PyErr_Occurred() )
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "%s", GetPyExceptionString().c_str());
        return false;
    }

    // Fetch user computation function
    if( !osPythonModule.empty() )
    {
        PyObject* poUserModule = PyImport_ImportModule(osPythonModule);
        if (poUserModule == nullptr || PyErr_Occurred())
        {
            CPLError(CE_Failure, CPLE_AppDefined,
                 "%s", GetPyExceptionString().c_str());
            Py_DecRef(poModule);
            return false;
        }
        m_poPrivate->m_poUserFunction = PyObject_GetAttrString(poUserModule,
                                                            osPythonFunction );
        Py_DecRef(poUserModule);
    }
    else
    {
        m_poPrivate->m_poUserFunction = PyObject_GetAttrString(poModule,
                                            osPythonFunction );
    }
    if (m_poPrivate->m_poUserFunction == nullptr || PyErr_Occurred())
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "%s", GetPyExceptionString().c_str());
        Py_DecRef(poModule);
        return false;
    }
    if( !PyCallable_Check(m_poPrivate->m_poUserFunction) )
    {
        CPLError(CE_Failure, CPLE_AppDefined, "Object '%s' is not callable",
                 osPythonFunction.c_str());
        Py_DecRef(poModule);
        return false;
    }

    // Fetch our GDALCreateNumpyArray python function
    m_poPrivate->m_poGDALCreateNumpyArray =
        PyObject_GetAttrString(poModule, "GDALCreateNumpyArray" );
    if (m_poPrivate->m_poGDALCreateNumpyArray == nullptr || PyErr_Occurred())
    {
        // Shouldn't happen normally...
        CPLError(CE_Failure, CPLE_AppDefined,
                 "%s", GetPyExceptionString().c_str());
        Py_DecRef(poModule);
        return false;
    }
    Py_DecRef(poModule);

    m_poPrivate->m_bPythonInitializationSuccess = true;
    return true;
}

/************************************************************************/
/*                             IRasterIO()                              */
/************************************************************************/

/**
 * Read/write a region of image data for this band.
 *
 * Each of the sources for this derived band will be read and passed to
 * the derived band pixel function.  The pixel function is responsible
 * for applying whatever algorithm is necessary to generate this band's
 * pixels from the sources.
 *
 * The sources will be read using the transfer type specified for sources
 * using SetSourceTransferType().  If no transfer type has been set for
 * this derived band, the band's data type will be used as the transfer type.
 *
 * @see gdalrasterband
 *
 * @param eRWFlag Either GF_Read to read a region of data, or GT_Write to
 * write a region of data.
 *
 * @param nXOff The pixel offset to the top left corner of the region
 * of the band to be accessed.  This would be zero to start from the left side.
 *
 * @param nYOff The line offset to the top left corner of the region
 * of the band to be accessed.  This would be zero to start from the top.
 *
 * @param nXSize The width of the region of the band to be accessed in pixels.
 *
 * @param nYSize The height of the region of the band to be accessed in lines.
 *
 * @param pData The buffer into which the data should be read, or from which
 * it should be written.  This buffer must contain at least nBufXSize *
 * nBufYSize words of type eBufType.  It is organized in left to right,
 * top to bottom pixel order.  Spacing is controlled by the nPixelSpace,
 * and nLineSpace parameters.
 *
 * @param nBufXSize The width of the buffer image into which the desired
 * region is to be read, or from which it is to be written.
 *
 * @param nBufYSize The height of the buffer image into which the desired
 * region is to be read, or from which it is to be written.
 *
 * @param eBufType The type of the pixel values in the pData data buffer.  The
 * pixel values will automatically be translated to/from the GDALRasterBand
 * data type as needed.
 *
 * @param nPixelSpace The byte offset from the start of one pixel value in
 * pData to the start of the next pixel value within a scanline.  If defaulted
 * (0) the size of the datatype eBufType is used.
 *
 * @param nLineSpace The byte offset from the start of one scanline in
 * pData to the start of the next.  If defaulted the size of the datatype
 * eBufType * nBufXSize is used.
 *
 * @return CE_Failure if the access fails, otherwise CE_None.
 */
CPLErr VRTDerivedRasterBand::IRasterIO( GDALRWFlag eRWFlag,
                                        int nXOff, int nYOff, int nXSize,
                                        int nYSize, void * pData, int nBufXSize,
                                        int nBufYSize, GDALDataType eBufType,
                                        GSpacing nPixelSpace,
                                        GSpacing nLineSpace,
                                        GDALRasterIOExtraArg* psExtraArg )
{
    if( eRWFlag == GF_Write )
    {
        CPLError( CE_Failure, CPLE_AppDefined,
                  "Writing through VRTSourcedRasterBand is not supported." );
        return CE_Failure;
    }

    const int nBufTypeSize = GDALGetDataTypeSizeBytes(eBufType);
    GDALDataType eSrcType = eSourceTransferType;
    if( eSrcType == GDT_Unknown || eSrcType >= GDT_TypeCount ) {
        eSrcType = eBufType;
    }
    const int nSrcTypeSize = GDALGetDataTypeSizeBytes(eSrcType);

/* -------------------------------------------------------------------- */
/*      Initialize the buffer to some background value. Use the         */
/*      nodata value if available.                                      */
/* -------------------------------------------------------------------- */
    if( bSkipBufferInitialization )
    {
        // Do nothing
    }
    else if( nPixelSpace == nBufTypeSize &&
        (!m_bNoDataValueSet || m_dfNoDataValue == 0) ) {
        memset( pData, 0,
                static_cast<size_t>(nBufXSize * nBufYSize * nPixelSpace) );
    }
    else if( m_bNoDataValueSet )
    {
        double dfWriteValue = m_dfNoDataValue;

        for( int iLine = 0; iLine < nBufYSize; iLine++ )
        {
            GDALCopyWords(
                &dfWriteValue, GDT_Float64, 0,
                reinterpret_cast<GByte *>( pData ) + nLineSpace * iLine,
                eBufType, static_cast<int>(nPixelSpace), nBufXSize );
        }
    }

/* -------------------------------------------------------------------- */
/*      Do we have overviews that would be appropriate to satisfy       */
/*      this request?                                                   */
/* -------------------------------------------------------------------- */
    if( (nBufXSize < nXSize || nBufYSize < nYSize)
        && GetOverviewCount() > 0 )
    {
        if( OverviewRasterIO(
               eRWFlag, nXOff, nYOff, nXSize, nYSize,
               pData, nBufXSize, nBufYSize,
               eBufType, nPixelSpace, nLineSpace, psExtraArg ) == CE_None )
            return CE_None;
    }

    /* ---- Get pixel function for band ---- */
    GDALDerivedPixelFunc pfnPixelFunc = nullptr;

    if( EQUAL(m_poPrivate->m_osLanguage, "C") )
    {
        pfnPixelFunc = VRTDerivedRasterBand::GetPixelFunction(pszFuncName);
        if( pfnPixelFunc == nullptr )
        {
            CPLError( CE_Failure, CPLE_IllegalArg,
                    "VRTDerivedRasterBand::IRasterIO:"
                    "Derived band pixel function '%s' not registered.",
                    this->pszFuncName) ;
            return CE_Failure;
        }
    }

    /* TODO: It would be nice to use a MallocBlock function for each
       individual buffer that would recycle blocks of memory from a
       cache by reassigning blocks that are nearly the same size.
       A corresponding FreeBlock might only truly free if the total size
       of freed blocks gets to be too great of a percentage of the size
       of the allocated blocks. */

    // Get buffers for each source.
    const int nBufferRadius = m_poPrivate->m_nBufferRadius;
    if( nBufferRadius > (INT_MAX - nBufXSize) / 2 ||
        nBufferRadius > (INT_MAX - nBufYSize) / 2 )
    {
        return CE_Failure;
    }
    const int nExtBufXSize = nBufXSize + 2 * nBufferRadius;
    const int nExtBufYSize = nBufYSize + 2 * nBufferRadius;
    void **pBuffers
        = reinterpret_cast<void **>( CPLMalloc(sizeof(void *) * nSources) );
    for( int iSource = 0; iSource < nSources; iSource++ ) {
        pBuffers[iSource] =
            VSI_MALLOC3_VERBOSE(nSrcTypeSize, nExtBufXSize, nExtBufYSize);
        if( pBuffers[iSource] == nullptr )
        {
            for (int i = 0; i < iSource; i++) {
                VSIFree(pBuffers[i]);
            }
            CPLFree(pBuffers);
            return CE_Failure;
        }

        /* ------------------------------------------------------------ */
        /* #4045: Initialize the newly allocated buffers before handing */
        /* them off to the sources. These buffers are packed, so we     */
        /* don't need any special line-by-line handling when a nonzero  */
        /* nodata value is set.                                         */
        /* ------------------------------------------------------------ */
        if( !m_bNoDataValueSet || m_dfNoDataValue == 0 )
        {
            memset( pBuffers[iSource], 0, static_cast<size_t>(nSrcTypeSize) *
                    nExtBufXSize * nExtBufYSize );
        }
        else
        {
            GDALCopyWords( &m_dfNoDataValue, GDT_Float64, 0,
                           reinterpret_cast<GByte *>( pBuffers[iSource] ),
                           eSrcType, nSrcTypeSize,
                           nExtBufXSize * nExtBufYSize );
        }
    }

    GDALRasterIOExtraArg sExtraArg;
    GDALCopyRasterIOExtraArg(&sExtraArg, psExtraArg);

    int nXShiftInBuffer = 0;
    int nYShiftInBuffer = 0;
    int nExtBufXSizeReq = nExtBufXSize;
    int nExtBufYSizeReq = nExtBufYSize;

    int nXOffExt = nXOff;
    int nYOffExt = nYOff;
    int nXSizeExt = nXSize;
    int nYSizeExt = nYSize;

    if( nBufferRadius )
    {
        double dfXRatio = static_cast<double>(nXSize) / nBufXSize;
        double dfYRatio = static_cast<double>(nYSize) / nBufYSize;

        if( !sExtraArg.bFloatingPointWindowValidity )
        {
            sExtraArg.dfXOff = nXOff;
            sExtraArg.dfYOff = nYOff;
            sExtraArg.dfXSize = nXSize;
            sExtraArg.dfYSize = nYSize;
        }

        sExtraArg.dfXOff -= dfXRatio * nBufferRadius;
        sExtraArg.dfYOff -= dfYRatio * nBufferRadius;
        sExtraArg.dfXSize += 2 * dfXRatio * nBufferRadius;
        sExtraArg.dfYSize += 2 * dfYRatio * nBufferRadius;
        if( sExtraArg.dfXOff < 0 )
        {
            nXShiftInBuffer = -static_cast<int>(sExtraArg.dfXOff / dfXRatio);
            nExtBufXSizeReq -= nXShiftInBuffer;
            sExtraArg.dfXSize += sExtraArg.dfXOff;
            sExtraArg.dfXOff = 0;
        }
        if( sExtraArg.dfYOff < 0 )
        {
            nYShiftInBuffer = -static_cast<int>(sExtraArg.dfYOff / dfYRatio);
            nExtBufYSizeReq -= nYShiftInBuffer;
            sExtraArg.dfYSize += sExtraArg.dfYOff;
            sExtraArg.dfYOff = 0;
        }
        if( sExtraArg.dfXOff + sExtraArg.dfXSize > nRasterXSize )
        {
            nExtBufXSizeReq -= static_cast<int>((sExtraArg.dfXOff +
                        sExtraArg.dfXSize - nRasterXSize) / dfXRatio);
            sExtraArg.dfXSize = nRasterXSize - sExtraArg.dfXOff;
        }
        if( sExtraArg.dfYOff + sExtraArg.dfYSize > nRasterYSize )
        {
            nExtBufYSizeReq -= static_cast<int>((sExtraArg.dfYOff +
                        sExtraArg.dfYSize - nRasterYSize) / dfYRatio);
            sExtraArg.dfYSize = nRasterYSize - sExtraArg.dfYOff;
        }

        nXOffExt = static_cast<int>(sExtraArg.dfXOff);
        nYOffExt = static_cast<int>(sExtraArg.dfYOff);
        nXSizeExt = std::min(static_cast<int>(sExtraArg.dfXSize + 0.5),
                             nRasterXSize - nXOffExt);
        nYSizeExt = std::min(static_cast<int>(sExtraArg.dfYSize + 0.5),
                             nRasterYSize - nYOffExt);
    }

    // Load values for sources into packed buffers.
    CPLErr eErr = CE_None;
    for( int iSource = 0; iSource < nSources && eErr == CE_None; iSource++ ) {
        GByte* pabyBuffer = reinterpret_cast<GByte*>(pBuffers[iSource]);
        eErr = reinterpret_cast<VRTSource *>( papoSources[iSource] )->RasterIO(
            eSrcType,
            nXOffExt, nYOffExt, nXSizeExt, nYSizeExt,
            pabyBuffer + (nYShiftInBuffer * nExtBufXSize +
                                            nXShiftInBuffer) * nSrcTypeSize,
            nExtBufXSizeReq, nExtBufYSizeReq,
            eSrcType,
            nSrcTypeSize,
            nSrcTypeSize * nExtBufXSize,
            &sExtraArg );

        // Extend first lines
        for( int iY = 0; iY < nYShiftInBuffer; iY++ )
        {
            memcpy( pabyBuffer + iY * nExtBufXSize * nSrcTypeSize,
                    pabyBuffer + nYShiftInBuffer * nExtBufXSize * nSrcTypeSize,
                    nExtBufXSize * nSrcTypeSize );
        }
        // Extend last lines
        for( int iY = nYShiftInBuffer + nExtBufYSizeReq; iY < nExtBufYSize; iY++ )
        {
            memcpy( pabyBuffer + iY * nExtBufXSize * nSrcTypeSize,
                    pabyBuffer + (nYShiftInBuffer + nExtBufYSizeReq - 1) *
                                                    nExtBufXSize * nSrcTypeSize,
                    nExtBufXSize * nSrcTypeSize );
        }
        // Extend first cols
        if( nXShiftInBuffer )
        {
            for( int iY = 0; iY < nExtBufYSize; iY ++ )
            {
                for( int iX = 0; iX < nXShiftInBuffer; iX++ )
                {
                    memcpy( pabyBuffer + (iY * nExtBufXSize + iX) * nSrcTypeSize,
                            pabyBuffer + (iY * nExtBufXSize +
                                                nXShiftInBuffer) * nSrcTypeSize,
                            nSrcTypeSize );
                }
            }
        }
        // Extent last cols
        if( nXShiftInBuffer + nExtBufXSizeReq < nExtBufXSize )
        {
            for( int iY = 0; iY < nExtBufYSize; iY ++ )
            {
                for( int iX = nXShiftInBuffer + nExtBufXSizeReq;
                         iX < nExtBufXSize; iX++ )
                {
                    memcpy( pabyBuffer + (iY * nExtBufXSize + iX) * nSrcTypeSize,
                            pabyBuffer + (iY * nExtBufXSize + nXShiftInBuffer +
                                            nExtBufXSizeReq - 1) * nSrcTypeSize,
                            nSrcTypeSize );
                }
            }
        }
    }

    // Apply pixel function.
    if( eErr == CE_None && EQUAL(m_poPrivate->m_osLanguage, "Python") )
    {
        eErr = CE_Failure;

        // numpy doesn't have native cint16/cint32
        if( eSrcType == GDT_CInt16 || eSrcType == GDT_CInt32 )
        {
            CPLError(CE_Failure, CPLE_AppDefined,
                     "CInt16/CInt32 data type not supported for SourceTransferType");
            goto end;
        }
        if( eDataType == GDT_CInt16 || eDataType == GDT_CInt32 )
        {
            CPLError(CE_Failure, CPLE_AppDefined,
                     "CInt16/CInt32 data type not supported for data type");
            goto end;
        }

        if( !InitializePython() )
            goto end;

        GByte* pabyTmpBuffer = nullptr;
        // Do we need a temporary buffer or can we use directly the output
        // buffer ?
        if( nBufferRadius != 0 ||
            eDataType != eBufType ||
            nPixelSpace != nBufTypeSize ||
            nLineSpace != static_cast<GSpacing>(nBufTypeSize) * nBufXSize )
        {
            pabyTmpBuffer = reinterpret_cast<GByte*>(VSI_CALLOC_VERBOSE(
                            static_cast<size_t>(nExtBufXSize) * nExtBufYSize,
                            GDALGetDataTypeSizeBytes(eDataType)));
            if( !pabyTmpBuffer )
                goto end;
        }

        {
        const bool bUseExclusiveLock = m_poPrivate->m_bExclusiveLock ||
                    ( m_poPrivate->m_bFirstTime &&
                    m_poPrivate->m_osCode.find("@jit") != std::string::npos);
        m_poPrivate->m_bFirstTime = false;
        VRT_GIL_Holder oHolder(bUseExclusiveLock);

        // Prepare target numpy array
        PyObject* poPyDstArray = GDALCreateNumpyArray(
                                    m_poPrivate->m_poGDALCreateNumpyArray,
                                    pabyTmpBuffer ? pabyTmpBuffer : pData,
                                    eDataType,
                                    nExtBufYSize,
                                    nExtBufXSize);
        if( !poPyDstArray )
        {
            VSIFree(pabyTmpBuffer);
            goto end;
        }

        // Wrap source buffers as input numpy arrays
        PyObject* pyArgInputArray = PyTuple_New(nSources);
        for( int i = 0; i < nSources; i++ )
        {
            GByte* pabyBuffer = reinterpret_cast<GByte*>(pBuffers[i]);
            PyObject* poPySrcArray = GDALCreateNumpyArray(
                        m_poPrivate->m_poGDALCreateNumpyArray,
                        pabyBuffer,
                        eSrcType,
                        nExtBufYSize,
                        nExtBufXSize);
            CPLAssert(poPySrcArray);
            PyTuple_SetItem(pyArgInputArray, i, poPySrcArray);
        }

        // Create arguments
        PyObject* pyArgs = PyTuple_New(10);
        PyTuple_SetItem(pyArgs, 0, pyArgInputArray);
        PyTuple_SetItem(pyArgs, 1, poPyDstArray);
        PyTuple_SetItem(pyArgs, 2, PyInt_FromLong(nXOff));
        PyTuple_SetItem(pyArgs, 3, PyInt_FromLong(nYOff));
        PyTuple_SetItem(pyArgs, 4, PyInt_FromLong(nXSize));
        PyTuple_SetItem(pyArgs, 5, PyInt_FromLong(nYSize));
        PyTuple_SetItem(pyArgs, 6, PyInt_FromLong(nRasterXSize));
        PyTuple_SetItem(pyArgs, 7, PyInt_FromLong(nRasterYSize));
        PyTuple_SetItem(pyArgs, 8, PyInt_FromLong(nBufferRadius));

        double adfGeoTransform[6];
        adfGeoTransform[0] = 0;
        adfGeoTransform[1] = 1;
        adfGeoTransform[2] = 0;
        adfGeoTransform[3] = 0;
        adfGeoTransform[4] = 0;
        adfGeoTransform[5] = 1;
        if( GetDataset() )
            GetDataset()->GetGeoTransform(adfGeoTransform);
        PyObject* pyGT = PyTuple_New(6);
        for(int i = 0; i < 6; i++ )
            PyTuple_SetItem(pyGT, i, PyFloat_FromDouble(adfGeoTransform[i]));
        PyTuple_SetItem(pyArgs, 9, pyGT);

        // Prepare kwargs
        PyObject* pyKwargs = PyDict_New();
        for( size_t i = 0; i < m_poPrivate->m_oFunctionArgs.size(); ++i )
        {
            const char* pszKey =
                m_poPrivate->m_oFunctionArgs[i].first.c_str();
            const char* pszValue =
                m_poPrivate->m_oFunctionArgs[i].second.c_str();
            PyDict_SetItemString(pyKwargs, pszKey,
                PyString_FromStringAndSize(pszValue, strlen(pszValue)));
        }

        // Call user function
        PyObject* pRetValue = PyObject_Call(
                                        m_poPrivate->m_poUserFunction,
                                        pyArgs, pyKwargs);

        Py_DecRef(pyArgs);
        Py_DecRef(pyKwargs);

        if (PyErr_Occurred())
        {
            CPLError(CE_Failure, CPLE_AppDefined,
                     "%s", GetPyExceptionString().c_str());
        }
        else
        {
            eErr = CE_None;
        }
        if( pRetValue )
            Py_DecRef(pRetValue);
        } // End of GIL section

        if( pabyTmpBuffer )
        {
            // Copy numpy destination array to user buffer
            for( int iY = 0; iY < nBufYSize; iY++ )
            {
                size_t nSrcOffset = (static_cast<size_t>(iY + nBufferRadius) *
                    nExtBufXSize + nBufferRadius) *
                    GDALGetDataTypeSizeBytes(eDataType);
                GDALCopyWords(pabyTmpBuffer + nSrcOffset,
                              eDataType,
                              GDALGetDataTypeSizeBytes(eDataType),
                              reinterpret_cast<GByte*>(pData) + iY * nLineSpace,
                              eBufType,
                              static_cast<int>(nPixelSpace),
                              nBufXSize);
            }

            VSIFree(pabyTmpBuffer);
        }
    }
    else if( eErr == CE_None && pfnPixelFunc != nullptr ) {
        eErr = pfnPixelFunc( reinterpret_cast<void **>( pBuffers ), nSources,
                             pData, nBufXSize, nBufYSize,
                             eSrcType, eBufType, static_cast<int>(nPixelSpace),
                             static_cast<int>(nLineSpace) );
    }
end:
    // Release buffers.
    for ( int iSource = 0; iSource < nSources; iSource++ ) {
        VSIFree(pBuffers[iSource]);
    }
    CPLFree(pBuffers);

    return eErr;
}

/************************************************************************/
/*                         IGetDataCoverageStatus()                     */
/************************************************************************/

int  VRTDerivedRasterBand::IGetDataCoverageStatus( int /* nXOff */,
                                                   int /* nYOff */,
                                                   int /* nXSize */,
                                                   int /* nYSize */,
                                                   int /* nMaskFlagStop */,
                                                   double* pdfDataPct)
{
    if( pdfDataPct != nullptr )
        *pdfDataPct = -1.0;
    return GDAL_DATA_COVERAGE_STATUS_UNIMPLEMENTED | GDAL_DATA_COVERAGE_STATUS_DATA;
}

/************************************************************************/
/*                              XMLInit()                               */
/************************************************************************/

CPLErr VRTDerivedRasterBand::XMLInit( CPLXMLNode *psTree,
                                      const char *pszVRTPath,
                                      void* pUniqueHandle )

{
    const CPLErr eErr = VRTSourcedRasterBand::XMLInit( psTree, pszVRTPath,
                                                       pUniqueHandle );
    if( eErr != CE_None )
        return eErr;

    // Read derived pixel function type.
    SetPixelFunctionName( CPLGetXMLValue( psTree, "PixelFunctionType", nullptr ) );
    if( pszFuncName == nullptr || EQUAL(pszFuncName, "") )
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "PixelFunctionType missing");
        return CE_Failure;
    }

    m_poPrivate->m_osLanguage = CPLGetXMLValue( psTree,
                                                "PixelFunctionLanguage", "C" );
    if( !EQUAL(m_poPrivate->m_osLanguage, "C") &&
        !EQUAL(m_poPrivate->m_osLanguage, "Python") )
    {
        CPLError(CE_Failure, CPLE_NotSupported,
                 "Unsupported PixelFunctionLanguage");
        return CE_Failure;
    }

    m_poPrivate->m_osCode =
                        CPLGetXMLValue( psTree, "PixelFunctionCode", "" );
    if( !m_poPrivate->m_osCode.empty() &&
        !EQUAL(m_poPrivate->m_osLanguage, "Python") )
    {
        CPLError(CE_Failure, CPLE_NotSupported,
                 "PixelFunctionCode can only be used with Python");
        return CE_Failure;
    }

    m_poPrivate->m_nBufferRadius =
                        atoi(CPLGetXMLValue( psTree, "BufferRadius", "0" ));
    if( m_poPrivate->m_nBufferRadius < 0 ||
        m_poPrivate->m_nBufferRadius > 1024 )
    {
        CPLError(CE_Failure, CPLE_AppDefined, "Invalid value for BufferRadius");
        return CE_Failure;
    }
    if( m_poPrivate->m_nBufferRadius != 0 &&
        !EQUAL(m_poPrivate->m_osLanguage, "Python") )
    {
        CPLError(CE_Failure, CPLE_NotSupported,
                 "BufferRadius can only be used with Python");
        return CE_Failure;
    }

    CPLXMLNode* psArgs = CPLGetXMLNode( psTree, "PixelFunctionArguments" );
    if( psArgs != nullptr )
    {
        if( !EQUAL(m_poPrivate->m_osLanguage, "Python") )
        {
            CPLError(CE_Failure, CPLE_NotSupported,
                     "PixelFunctionArguments can only be used with Python");
            return CE_Failure;
        }
        for( CPLXMLNode* psIter = psArgs->psChild;
                         psIter != nullptr;
                         psIter = psIter->psNext )
        {
            if( psIter->eType == CXT_Attribute )
            {
                m_poPrivate->m_oFunctionArgs.push_back(
                    std::pair<CPLString,CPLString>(psIter->pszValue,
                                                   psIter->psChild->pszValue));
            }
        }
    }

    // Read optional source transfer data type.
    const char *pszTypeName = CPLGetXMLValue(psTree, "SourceTransferType", nullptr);
    if( pszTypeName != nullptr )
    {
        eSourceTransferType = GDALGetDataTypeByName( pszTypeName );
    }

    return CE_None;
}

/************************************************************************/
/*                           SerializeToXML()                           */
/************************************************************************/

CPLXMLNode *VRTDerivedRasterBand::SerializeToXML( const char *pszVRTPath )
{
    CPLXMLNode *psTree = VRTSourcedRasterBand::SerializeToXML( pszVRTPath );

/* -------------------------------------------------------------------- */
/*      Set subclass.                                                   */
/* -------------------------------------------------------------------- */
    CPLCreateXMLNode(
        CPLCreateXMLNode( psTree, CXT_Attribute, "subClass" ),
        CXT_Text, "VRTDerivedRasterBand" );

    /* ---- Encode DerivedBand-specific fields ---- */
    if( !EQUAL( m_poPrivate->m_osLanguage, "C" ) )
    {
        CPLSetXMLValue( psTree, "PixelFunctionLanguage",
                        m_poPrivate->m_osLanguage );
    }
    if( pszFuncName != nullptr && strlen(pszFuncName) > 0 )
        CPLSetXMLValue( psTree, "PixelFunctionType", pszFuncName );
    if( !m_poPrivate->m_oFunctionArgs.empty() )
    {
        CPLXMLNode* psArgs = CPLCreateXMLNode( psTree, CXT_Element,
                                               "PixelFunctionArguments" );
        for( size_t i = 0; i < m_poPrivate->m_oFunctionArgs.size(); ++i )
        {
            const char* pszKey =
                m_poPrivate->m_oFunctionArgs[i].first.c_str();
            const char* pszValue =
                m_poPrivate->m_oFunctionArgs[i].second.c_str();
            CPLCreateXMLNode(
                CPLCreateXMLNode( psArgs, CXT_Attribute, pszKey ),
                                  CXT_Text, pszValue );
        }
    }
    if( !m_poPrivate->m_osCode.empty() )
    {
        if( m_poPrivate->m_osCode.find("<![CDATA[") == std::string::npos )
        {
            CPLCreateXMLNode(
                CPLCreateXMLNode( psTree,
                                  CXT_Element, "PixelFunctionCode" ),
                 CXT_Literal,
                 ("<![CDATA[" + m_poPrivate->m_osCode + "]]>").c_str() );
        }
        else
        {
            CPLSetXMLValue( psTree, "PixelFunctionCode",
                            m_poPrivate->m_osCode );
        }
    }
    if( m_poPrivate->m_nBufferRadius != 0 )
        CPLSetXMLValue( psTree, "BufferRadius",
                        CPLSPrintf("%d",m_poPrivate->m_nBufferRadius) );
    if( this->eSourceTransferType != GDT_Unknown)
        CPLSetXMLValue( psTree, "SourceTransferType",
                        GDALGetDataTypeName( eSourceTransferType ) );

    return psTree;
}

/************************************************************************/
/*                             GetMinimum()                             */
/************************************************************************/

double VRTDerivedRasterBand::GetMinimum( int *pbSuccess )
{
    return GDALRasterBand::GetMinimum(pbSuccess);
}

/************************************************************************/
/*                             GetMaximum()                             */
/************************************************************************/

double VRTDerivedRasterBand::GetMaximum( int *pbSuccess )
{
    return GDALRasterBand::GetMaximum(pbSuccess);
}

/************************************************************************/
/*                       ComputeRasterMinMax()                          */
/************************************************************************/

CPLErr VRTDerivedRasterBand::ComputeRasterMinMax( int bApproxOK, double* adfMinMax )
{
    return GDALRasterBand::ComputeRasterMinMax( bApproxOK, adfMinMax );
}

/************************************************************************/
/*                         ComputeStatistics()                          */
/************************************************************************/

CPLErr
VRTDerivedRasterBand::ComputeStatistics( int bApproxOK,
                                         double *pdfMin, double *pdfMax,
                                         double *pdfMean, double *pdfStdDev,
                                         GDALProgressFunc pfnProgress,
                                         void *pProgressData )

{
    return GDALRasterBand::ComputeStatistics(  bApproxOK,
                                            pdfMin, pdfMax,
                                            pdfMean, pdfStdDev,
                                            pfnProgress, pProgressData );
}

/************************************************************************/
/*                            GetHistogram()                            */
/************************************************************************/

CPLErr VRTDerivedRasterBand::GetHistogram( double dfMin, double dfMax,
                                           int nBuckets, GUIntBig *panHistogram,
                                           int bIncludeOutOfRange, int bApproxOK,
                                           GDALProgressFunc pfnProgress,
                                           void *pProgressData )

{
    return VRTRasterBand::GetHistogram( dfMin, dfMax,
                                            nBuckets, panHistogram,
                                            bIncludeOutOfRange, bApproxOK,
                                            pfnProgress, pProgressData );
}

/*! @endcond */