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