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 184 : VRTDerivedRasterBandPrivateData() = default;
154 :
155 368 : virtual ~VRTDerivedRasterBandPrivateData()
156 184 : {
157 184 : if (m_poGDALCreateNumpyArray)
158 43 : Py_DecRef(m_poGDALCreateNumpyArray);
159 184 : if (m_poUserFunction)
160 44 : Py_DecRef(m_poUserFunction);
161 368 : }
162 : };
163 :
164 : /************************************************************************/
165 : /* ==================================================================== */
166 : /* VRTDerivedRasterBand */
167 : /* ==================================================================== */
168 : /************************************************************************/
169 :
170 : /************************************************************************/
171 : /* VRTDerivedRasterBand() */
172 : /************************************************************************/
173 :
174 163 : VRTDerivedRasterBand::VRTDerivedRasterBand(GDALDataset *poDSIn, int nBandIn)
175 : : VRTSourcedRasterBand(poDSIn, nBandIn), m_poPrivate(nullptr),
176 163 : pszFuncName(nullptr), eSourceTransferType(GDT_Unknown)
177 : {
178 163 : m_poPrivate = new VRTDerivedRasterBandPrivateData;
179 163 : }
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 368 : VRTDerivedRasterBand::~VRTDerivedRasterBand()
200 :
201 : {
202 184 : CPLFree(pszFuncName);
203 184 : delete m_poPrivate;
204 368 : }
205 :
206 : /************************************************************************/
207 : /* Cleanup() */
208 : /************************************************************************/
209 :
210 933 : void VRTDerivedRasterBand::Cleanup()
211 : {
212 933 : }
213 :
214 : /************************************************************************/
215 : /* GetGlobalMapPixelFunction() */
216 : /************************************************************************/
217 :
218 : static std::map<std::string,
219 : std::pair<VRTDerivedRasterBand::PixelFunc, std::string>> &
220 36286 : GetGlobalMapPixelFunction()
221 : {
222 : static std::map<std::string,
223 : std::pair<VRTDerivedRasterBand::PixelFunc, std::string>>
224 36286 : gosMapPixelFunction;
225 36286 : 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 19382 : CPLErr CPL_STDCALL GDALAddDerivedBandPixelFunc(
251 : const char *pszName, GDALDerivedPixelFunc pfnNewFunction)
252 : {
253 19382 : if (pszName == nullptr || pszName[0] == '\0' || pfnNewFunction == nullptr)
254 : {
255 0 : return CE_None;
256 : }
257 :
258 38764 : GetGlobalMapPixelFunction()[pszName] = {
259 58 : [pfnNewFunction](void **papoSources, int nSources, void *pData,
260 : int nBufXSize, int nBufYSize, GDALDataType eSrcType,
261 : GDALDataType eBufType, int nPixelSpace, int nLineSpace,
262 58 : CSLConstList papszFunctionArgs)
263 : {
264 : (void)papszFunctionArgs;
265 58 : return pfnNewFunction(papoSources, nSources, pData, nBufXSize,
266 : nBufYSize, eSrcType, eBufType, nPixelSpace,
267 58 : nLineSpace);
268 : },
269 58146 : ""};
270 :
271 19382 : 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 16798 : CPLErr CPL_STDCALL GDALAddDerivedBandPixelFuncWithArgs(
293 : const char *pszName, GDALDerivedPixelFuncWithArgs pfnNewFunction,
294 : const char *pszMetadata)
295 : {
296 16798 : if (!pszName || pszName[0] == '\0' || !pfnNewFunction)
297 : {
298 0 : return CE_None;
299 : }
300 :
301 33596 : GetGlobalMapPixelFunction()[pszName] = {pfnNewFunction,
302 50394 : pszMetadata ? pszMetadata : ""};
303 :
304 16798 : 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 106 : VRTDerivedRasterBand::GetPixelFunction(const char *pszFuncNameIn)
352 : {
353 106 : if (pszFuncNameIn == nullptr || pszFuncNameIn[0] == '\0')
354 : {
355 0 : return nullptr;
356 : }
357 :
358 106 : const auto &oMapPixelFunction = GetGlobalMapPixelFunction();
359 106 : const auto oIter = oMapPixelFunction.find(pszFuncNameIn);
360 :
361 106 : if (oIter == oMapPixelFunction.end())
362 1 : return nullptr;
363 :
364 105 : 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 185 : void VRTDerivedRasterBand::SetPixelFunctionName(const char *pszFuncNameIn)
379 : {
380 185 : CPLFree(pszFuncName);
381 185 : pszFuncName = CPLStrdup(pszFuncNameIn);
382 185 : }
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 46 : CPLErr VRTDerivedRasterBand::GetPixelFunctionArguments(
759 : const CPLString &osMetadata,
760 : std::vector<std::pair<CPLString, CPLString>> &oAdditionalArgs)
761 : {
762 :
763 92 : auto poArgs = CPLXMLTreeCloser(CPLParseXMLString(osMetadata));
764 92 : if (poArgs != nullptr && poArgs->eType == CXT_Element &&
765 46 : !strcmp(poArgs->pszValue, "PixelFunctionArgumentsList"))
766 : {
767 117 : for (CPLXMLNode *psIter = poArgs->psChild; psIter != nullptr;
768 71 : psIter = psIter->psNext)
769 : {
770 72 : if (psIter->eType == CXT_Element &&
771 72 : !strcmp(psIter->pszValue, "Argument"))
772 : {
773 72 : CPLString osName, osType, osValue;
774 72 : auto pszName = CPLGetXMLValue(psIter, "name", nullptr);
775 72 : if (pszName != nullptr)
776 61 : osName = pszName;
777 72 : auto pszType = CPLGetXMLValue(psIter, "type", nullptr);
778 72 : if (pszType != nullptr)
779 72 : osType = pszType;
780 72 : auto pszValue = CPLGetXMLValue(psIter, "value", nullptr);
781 72 : if (pszValue != nullptr)
782 16 : osValue = pszValue;
783 72 : if (osType == "constant" && osValue != "" && osName != "")
784 1 : oAdditionalArgs.push_back(
785 2 : std::pair<CPLString, CPLString>(osName, osValue));
786 72 : 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 45 : 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 493 : 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 493 : 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 492 : const int nBufTypeSize = GDALGetDataTypeSizeBytes(eBufType);
897 492 : GDALDataType eSrcType = eSourceTransferType;
898 492 : if (eSrcType == GDT_Unknown || eSrcType >= GDT_TypeCount)
899 : {
900 409 : eSrcType = eBufType;
901 : }
902 492 : const int nSrcTypeSize = GDALGetDataTypeSizeBytes(eSrcType);
903 :
904 : /* -------------------------------------------------------------------- */
905 : /* Initialize the buffer to some background value. Use the */
906 : /* nodata value if available. */
907 : /* -------------------------------------------------------------------- */
908 492 : if (SkipBufferInitialization())
909 : {
910 : // Do nothing
911 : }
912 437 : else if (nPixelSpace == nBufTypeSize &&
913 437 : (!m_bNoDataValueSet || m_dfNoDataValue == 0))
914 : {
915 436 : memset(pData, 0,
916 436 : 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 492 : 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 492 : const std::pair<PixelFunc, std::string> *poPixelFunc = nullptr;
944 984 : std::vector<std::pair<CPLString, CPLString>> oAdditionalArgs;
945 :
946 492 : if (EQUAL(m_poPrivate->m_osLanguage, "C"))
947 : {
948 106 : poPixelFunc = VRTDerivedRasterBand::GetPixelFunction(pszFuncName);
949 106 : 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 105 : if (poPixelFunc->second != "")
959 : {
960 92 : if (GetPixelFunctionArguments(poPixelFunc->second,
961 46 : 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 490 : const int nBufferRadius = m_poPrivate->m_nBufferRadius;
977 490 : if (nBufferRadius > (INT_MAX - nBufXSize) / 2 ||
978 490 : nBufferRadius > (INT_MAX - nBufYSize) / 2)
979 : {
980 0 : return CE_Failure;
981 : }
982 490 : const int nExtBufXSize = nBufXSize + 2 * nBufferRadius;
983 490 : const int nExtBufYSize = nBufYSize + 2 * nBufferRadius;
984 490 : int nBufferCount = 0;
985 : void **pBuffers =
986 490 : static_cast<void **>(CPLMalloc(sizeof(void *) * nSources));
987 980 : std::vector<int> anMapBufferIdxToSourceIdx(nSources);
988 667 : for (int iSource = 0; iSource < nSources; iSource++)
989 : {
990 185 : 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 173 : anMapBufferIdxToSourceIdx[nBufferCount] = iSource;
1021 346 : pBuffers[nBufferCount] =
1022 173 : VSI_MALLOC3_VERBOSE(nSrcTypeSize, nExtBufXSize, nExtBufYSize);
1023 173 : 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 173 : if (!m_bNoDataValueSet || m_dfNoDataValue == 0)
1040 : {
1041 170 : memset(pBuffers[nBufferCount], 0,
1042 170 : static_cast<size_t>(nSrcTypeSize) * nExtBufXSize *
1043 170 : 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 173 : ++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 490 : if (nBufferCount == 0 && m_poPrivate->m_bSkipNonContributingSources)
1059 : {
1060 1 : CPLFree(pBuffers);
1061 1 : return CE_None;
1062 : }
1063 :
1064 : GDALRasterIOExtraArg sExtraArg;
1065 489 : GDALCopyRasterIOExtraArg(&sExtraArg, psExtraArg);
1066 :
1067 489 : int nXShiftInBuffer = 0;
1068 489 : int nYShiftInBuffer = 0;
1069 489 : int nExtBufXSizeReq = nExtBufXSize;
1070 489 : int nExtBufYSizeReq = nExtBufYSize;
1071 :
1072 489 : int nXOffExt = nXOff;
1073 489 : int nYOffExt = nYOff;
1074 489 : int nXSizeExt = nXSize;
1075 489 : int nYSizeExt = nYSize;
1076 :
1077 489 : 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 489 : CPLErr eErr = CE_None;
1133 489 : VRTSource::WorkingState oWorkingState;
1134 662 : for (int iBuffer = 0; iBuffer < nBufferCount && eErr == CE_None; iBuffer++)
1135 : {
1136 173 : const int iSource = anMapBufferIdxToSourceIdx[iBuffer];
1137 173 : GByte *pabyBuffer = static_cast<GByte *>(pBuffers[iBuffer]);
1138 173 : eErr = static_cast<VRTSource *>(papoSources[iSource])
1139 346 : ->RasterIO(
1140 : eSrcType, nXOffExt, nYOffExt, nXSizeExt, nYSizeExt,
1141 173 : pabyBuffer +
1142 173 : (nYShiftInBuffer * nExtBufXSize + nXShiftInBuffer) *
1143 : nSrcTypeSize,
1144 : nExtBufXSizeReq, nExtBufYSizeReq, eSrcType, nSrcTypeSize,
1145 173 : static_cast<GSpacing>(nSrcTypeSize) * nExtBufXSize,
1146 173 : &sExtraArg, oWorkingState);
1147 :
1148 : // Extend first lines
1149 182 : 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 182 : 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 173 : 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 173 : 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 489 : 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 41177 : for (int iY = 0; iY < nBufYSize; iY++)
1338 : {
1339 : size_t nSrcOffset =
1340 40837 : (static_cast<size_t>(iY + nBufferRadius) * nExtBufXSize +
1341 40837 : nBufferRadius) *
1342 40837 : GDALGetDataTypeSizeBytes(eDataType);
1343 40830 : GDALCopyWords(pabyTmpBuffer + nSrcOffset, eDataType,
1344 : GDALGetDataTypeSizeBytes(eDataType),
1345 40824 : static_cast<GByte *>(pData) + iY * nLineSpace,
1346 : eBufType, static_cast<int>(nPixelSpace),
1347 : nBufXSize);
1348 : }
1349 :
1350 340 : VSIFree(pabyTmpBuffer);
1351 : }
1352 : }
1353 104 : else if (eErr == CE_None && poPixelFunc != nullptr)
1354 : {
1355 104 : char **papszArgs = nullptr;
1356 :
1357 104 : oAdditionalArgs.insert(oAdditionalArgs.end(),
1358 104 : m_poPrivate->m_oFunctionArgs.begin(),
1359 208 : m_poPrivate->m_oFunctionArgs.end());
1360 152 : for (const auto &oArg : oAdditionalArgs)
1361 : {
1362 48 : const char *pszKey = oArg.first.c_str();
1363 48 : const char *pszValue = oArg.second.c_str();
1364 48 : papszArgs = CSLSetNameValue(papszArgs, pszKey, pszValue);
1365 : }
1366 :
1367 104 : eErr = (poPixelFunc->first)(
1368 : static_cast<void **>(pBuffers), nBufferCount, pData, nBufXSize,
1369 : nBufYSize, eSrcType, eBufType, static_cast<int>(nPixelSpace),
1370 : static_cast<int>(nLineSpace), papszArgs);
1371 :
1372 104 : CSLDestroy(papszArgs);
1373 : }
1374 0 : end:
1375 : // Release buffers.
1376 662 : for (int iBuffer = 0; iBuffer < nBufferCount; iBuffer++)
1377 : {
1378 173 : VSIFree(pBuffers[iBuffer]);
1379 : }
1380 489 : CPLFree(pBuffers);
1381 :
1382 489 : return eErr;
1383 : }
1384 :
1385 : /************************************************************************/
1386 : /* IGetDataCoverageStatus() */
1387 : /************************************************************************/
1388 :
1389 1 : int VRTDerivedRasterBand::IGetDataCoverageStatus(
1390 : int /* nXOff */, int /* nYOff */, int /* nXSize */, int /* nYSize */,
1391 : int /* nMaskFlagStop */, double *pdfDataPct)
1392 : {
1393 1 : if (pdfDataPct != nullptr)
1394 0 : *pdfDataPct = -1.0;
1395 : return GDAL_DATA_COVERAGE_STATUS_UNIMPLEMENTED |
1396 1 : GDAL_DATA_COVERAGE_STATUS_DATA;
1397 : }
1398 :
1399 : /************************************************************************/
1400 : /* XMLInit() */
1401 : /************************************************************************/
1402 :
1403 163 : CPLErr VRTDerivedRasterBand::XMLInit(const CPLXMLNode *psTree,
1404 : const char *pszVRTPath,
1405 : VRTMapSharedResources &oMapSharedSources)
1406 :
1407 : {
1408 : const CPLErr eErr =
1409 163 : VRTSourcedRasterBand::XMLInit(psTree, pszVRTPath, oMapSharedSources);
1410 163 : if (eErr != CE_None)
1411 0 : return eErr;
1412 :
1413 : // Read derived pixel function type.
1414 163 : SetPixelFunctionName(CPLGetXMLValue(psTree, "PixelFunctionType", nullptr));
1415 163 : if (pszFuncName == nullptr || EQUAL(pszFuncName, ""))
1416 : {
1417 1 : CPLError(CE_Failure, CPLE_AppDefined, "PixelFunctionType missing");
1418 1 : return CE_Failure;
1419 : }
1420 :
1421 162 : m_poPrivate->m_osLanguage =
1422 162 : CPLGetXMLValue(psTree, "PixelFunctionLanguage", "C");
1423 229 : if (!EQUAL(m_poPrivate->m_osLanguage, "C") &&
1424 67 : !EQUAL(m_poPrivate->m_osLanguage, "Python"))
1425 : {
1426 1 : CPLError(CE_Failure, CPLE_NotSupported,
1427 : "Unsupported PixelFunctionLanguage");
1428 1 : return CE_Failure;
1429 : }
1430 :
1431 161 : m_poPrivate->m_osCode = CPLGetXMLValue(psTree, "PixelFunctionCode", "");
1432 198 : if (!m_poPrivate->m_osCode.empty() &&
1433 37 : !EQUAL(m_poPrivate->m_osLanguage, "Python"))
1434 : {
1435 1 : CPLError(CE_Failure, CPLE_NotSupported,
1436 : "PixelFunctionCode can only be used with Python");
1437 1 : return CE_Failure;
1438 : }
1439 :
1440 160 : m_poPrivate->m_nBufferRadius =
1441 160 : atoi(CPLGetXMLValue(psTree, "BufferRadius", "0"));
1442 160 : if (m_poPrivate->m_nBufferRadius < 0 || m_poPrivate->m_nBufferRadius > 1024)
1443 : {
1444 1 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid value for BufferRadius");
1445 1 : return CE_Failure;
1446 : }
1447 170 : if (m_poPrivate->m_nBufferRadius != 0 &&
1448 11 : !EQUAL(m_poPrivate->m_osLanguage, "Python"))
1449 : {
1450 1 : CPLError(CE_Failure, CPLE_NotSupported,
1451 : "BufferRadius can only be used with Python");
1452 1 : return CE_Failure;
1453 : }
1454 :
1455 : const CPLXMLNode *const psArgs =
1456 158 : CPLGetXMLNode(psTree, "PixelFunctionArguments");
1457 158 : if (psArgs != nullptr)
1458 : {
1459 90 : for (const CPLXMLNode *psIter = psArgs->psChild; psIter;
1460 58 : psIter = psIter->psNext)
1461 : {
1462 58 : if (psIter->eType == CXT_Attribute)
1463 : {
1464 58 : m_poPrivate->m_oFunctionArgs.push_back(
1465 58 : std::pair<CPLString, CPLString>(psIter->pszValue,
1466 58 : psIter->psChild->pszValue));
1467 : }
1468 : }
1469 : }
1470 :
1471 : // Read optional source transfer data type.
1472 : const char *pszTypeName =
1473 158 : CPLGetXMLValue(psTree, "SourceTransferType", nullptr);
1474 158 : if (pszTypeName != nullptr)
1475 : {
1476 74 : eSourceTransferType = GDALGetDataTypeByName(pszTypeName);
1477 : }
1478 :
1479 : // Whether to skip non contributing sources
1480 : const char *pszSkipNonContributiongSources =
1481 158 : CPLGetXMLValue(psTree, "SkipNonContributingSources", nullptr);
1482 158 : if (pszSkipNonContributiongSources)
1483 : {
1484 2 : m_poPrivate->m_bSkipNonContributingSourcesSpecified = true;
1485 4 : m_poPrivate->m_bSkipNonContributingSources =
1486 2 : CPLTestBool(pszSkipNonContributiongSources);
1487 : }
1488 :
1489 158 : return CE_None;
1490 : }
1491 :
1492 : /************************************************************************/
1493 : /* SerializeToXML() */
1494 : /************************************************************************/
1495 :
1496 7 : CPLXMLNode *VRTDerivedRasterBand::SerializeToXML(const char *pszVRTPath,
1497 : bool &bHasWarnedAboutRAMUsage,
1498 : size_t &nAccRAMUsage)
1499 : {
1500 7 : CPLXMLNode *psTree = VRTSourcedRasterBand::SerializeToXML(
1501 : pszVRTPath, bHasWarnedAboutRAMUsage, nAccRAMUsage);
1502 :
1503 : /* -------------------------------------------------------------------- */
1504 : /* Set subclass. */
1505 : /* -------------------------------------------------------------------- */
1506 7 : CPLCreateXMLNode(CPLCreateXMLNode(psTree, CXT_Attribute, "subClass"),
1507 : CXT_Text, "VRTDerivedRasterBand");
1508 :
1509 : /* ---- Encode DerivedBand-specific fields ---- */
1510 7 : if (!EQUAL(m_poPrivate->m_osLanguage, "C"))
1511 : {
1512 5 : CPLSetXMLValue(psTree, "PixelFunctionLanguage",
1513 5 : m_poPrivate->m_osLanguage);
1514 : }
1515 7 : if (pszFuncName != nullptr && strlen(pszFuncName) > 0)
1516 6 : CPLSetXMLValue(psTree, "PixelFunctionType", pszFuncName);
1517 7 : if (!m_poPrivate->m_oFunctionArgs.empty())
1518 : {
1519 : CPLXMLNode *psArgs =
1520 1 : CPLCreateXMLNode(psTree, CXT_Element, "PixelFunctionArguments");
1521 3 : for (size_t i = 0; i < m_poPrivate->m_oFunctionArgs.size(); ++i)
1522 : {
1523 2 : const char *pszKey = m_poPrivate->m_oFunctionArgs[i].first.c_str();
1524 : const char *pszValue =
1525 2 : m_poPrivate->m_oFunctionArgs[i].second.c_str();
1526 2 : CPLCreateXMLNode(CPLCreateXMLNode(psArgs, CXT_Attribute, pszKey),
1527 : CXT_Text, pszValue);
1528 : }
1529 : }
1530 7 : if (!m_poPrivate->m_osCode.empty())
1531 : {
1532 4 : if (m_poPrivate->m_osCode.find("<![CDATA[") == std::string::npos)
1533 : {
1534 4 : CPLCreateXMLNode(
1535 : CPLCreateXMLNode(psTree, CXT_Element, "PixelFunctionCode"),
1536 : CXT_Literal,
1537 8 : ("<![CDATA[" + m_poPrivate->m_osCode + "]]>").c_str());
1538 : }
1539 : else
1540 : {
1541 0 : CPLSetXMLValue(psTree, "PixelFunctionCode", m_poPrivate->m_osCode);
1542 : }
1543 : }
1544 7 : if (m_poPrivate->m_nBufferRadius != 0)
1545 1 : CPLSetXMLValue(psTree, "BufferRadius",
1546 1 : CPLSPrintf("%d", m_poPrivate->m_nBufferRadius));
1547 7 : if (this->eSourceTransferType != GDT_Unknown)
1548 2 : CPLSetXMLValue(psTree, "SourceTransferType",
1549 : GDALGetDataTypeName(eSourceTransferType));
1550 :
1551 7 : if (m_poPrivate->m_bSkipNonContributingSourcesSpecified)
1552 : {
1553 1 : CPLSetXMLValue(psTree, "SkipNonContributingSources",
1554 1 : m_poPrivate->m_bSkipNonContributingSources ? "true"
1555 : : "false");
1556 : }
1557 :
1558 7 : return psTree;
1559 : }
1560 :
1561 : /************************************************************************/
1562 : /* GetMinimum() */
1563 : /************************************************************************/
1564 :
1565 5 : double VRTDerivedRasterBand::GetMinimum(int *pbSuccess)
1566 : {
1567 5 : return GDALRasterBand::GetMinimum(pbSuccess);
1568 : }
1569 :
1570 : /************************************************************************/
1571 : /* GetMaximum() */
1572 : /************************************************************************/
1573 :
1574 5 : double VRTDerivedRasterBand::GetMaximum(int *pbSuccess)
1575 : {
1576 5 : return GDALRasterBand::GetMaximum(pbSuccess);
1577 : }
1578 :
1579 : /************************************************************************/
1580 : /* ComputeRasterMinMax() */
1581 : /************************************************************************/
1582 :
1583 1 : CPLErr VRTDerivedRasterBand::ComputeRasterMinMax(int bApproxOK,
1584 : double *adfMinMax)
1585 : {
1586 1 : return GDALRasterBand::ComputeRasterMinMax(bApproxOK, adfMinMax);
1587 : }
1588 :
1589 : /************************************************************************/
1590 : /* ComputeStatistics() */
1591 : /************************************************************************/
1592 :
1593 1 : CPLErr VRTDerivedRasterBand::ComputeStatistics(int bApproxOK, double *pdfMin,
1594 : double *pdfMax, double *pdfMean,
1595 : double *pdfStdDev,
1596 : GDALProgressFunc pfnProgress,
1597 : void *pProgressData)
1598 :
1599 : {
1600 1 : return GDALRasterBand::ComputeStatistics(bApproxOK, pdfMin, pdfMax, pdfMean,
1601 : pdfStdDev, pfnProgress,
1602 1 : pProgressData);
1603 : }
1604 :
1605 : /************************************************************************/
1606 : /* GetHistogram() */
1607 : /************************************************************************/
1608 :
1609 1 : CPLErr VRTDerivedRasterBand::GetHistogram(double dfMin, double dfMax,
1610 : int nBuckets, GUIntBig *panHistogram,
1611 : int bIncludeOutOfRange, int bApproxOK,
1612 : GDALProgressFunc pfnProgress,
1613 : void *pProgressData)
1614 :
1615 : {
1616 1 : return VRTRasterBand::GetHistogram(dfMin, dfMax, nBuckets, panHistogram,
1617 : bIncludeOutOfRange, bApproxOK,
1618 1 : pfnProgress, pProgressData);
1619 : }
1620 :
1621 : /*! @endcond */
|
