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