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