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