Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Virtual GDAL Datasets
4 : * Purpose: Implementation of VRTSourcedRasterBand
5 : * Author: Frank Warmerdam <warmerdam@pobox.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2001, Frank Warmerdam <warmerdam@pobox.com>
9 : * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "cpl_port.h"
15 : #include "gdal_vrt.h"
16 : #include "vrtdataset.h"
17 :
18 : #include <algorithm>
19 : #include <cassert>
20 : #include <cmath>
21 : #include <cstddef>
22 : #include <cstdio>
23 : #include <cstdlib>
24 : #include <cstring>
25 : #include <limits>
26 : #include <mutex>
27 : #include <set>
28 : #include <string>
29 :
30 : #include "cpl_conv.h"
31 : #include "cpl_error.h"
32 : #include "cpl_hash_set.h"
33 : #include "cpl_minixml.h"
34 : #include "cpl_progress.h"
35 : #include "cpl_quad_tree.h"
36 : #include "cpl_string.h"
37 : #include "cpl_vsi.h"
38 : #include "gdal.h"
39 : #include "gdal_priv.h"
40 : #include "gdal_thread_pool.h"
41 : #include "ogr_geometry.h"
42 :
43 : /*! @cond Doxygen_Suppress */
44 :
45 : /************************************************************************/
46 : /* ==================================================================== */
47 : /* VRTSourcedRasterBand */
48 : /* ==================================================================== */
49 : /************************************************************************/
50 :
51 : /************************************************************************/
52 : /* VRTSourcedRasterBand() */
53 : /************************************************************************/
54 :
55 1247 : VRTSourcedRasterBand::VRTSourcedRasterBand(GDALDataset *poDSIn, int nBandIn)
56 : {
57 1247 : VRTRasterBand::Initialize(poDSIn->GetRasterXSize(),
58 : poDSIn->GetRasterYSize());
59 :
60 1247 : poDS = poDSIn;
61 1247 : nBand = nBandIn;
62 1247 : }
63 :
64 : /************************************************************************/
65 : /* VRTSourcedRasterBand() */
66 : /************************************************************************/
67 :
68 0 : VRTSourcedRasterBand::VRTSourcedRasterBand(GDALDataType eType, int nXSize,
69 0 : int nYSize)
70 : {
71 0 : VRTRasterBand::Initialize(nXSize, nYSize);
72 :
73 0 : eDataType = eType;
74 0 : }
75 :
76 : /************************************************************************/
77 : /* VRTSourcedRasterBand() */
78 : /************************************************************************/
79 :
80 317 : VRTSourcedRasterBand::VRTSourcedRasterBand(GDALDataset *poDSIn, int nBandIn,
81 : GDALDataType eType, int nXSize,
82 317 : int nYSize)
83 317 : : VRTSourcedRasterBand(poDSIn, nBandIn, eType, nXSize, nYSize, 0, 0)
84 : {
85 317 : }
86 :
87 : /************************************************************************/
88 : /* VRTSourcedRasterBand() */
89 : /************************************************************************/
90 :
91 6402 : VRTSourcedRasterBand::VRTSourcedRasterBand(GDALDataset *poDSIn, int nBandIn,
92 : GDALDataType eType, int nXSize,
93 : int nYSize, int nBlockXSizeIn,
94 6402 : int nBlockYSizeIn)
95 : {
96 6402 : VRTRasterBand::Initialize(nXSize, nYSize);
97 :
98 6402 : poDS = poDSIn;
99 6402 : nBand = nBandIn;
100 :
101 6402 : eDataType = eType;
102 6402 : if (nBlockXSizeIn > 0)
103 6065 : nBlockXSize = nBlockXSizeIn;
104 6402 : if (nBlockYSizeIn > 0)
105 6068 : nBlockYSize = nBlockYSizeIn;
106 6402 : }
107 :
108 : /************************************************************************/
109 : /* ~VRTSourcedRasterBand() */
110 : /************************************************************************/
111 :
112 15092 : VRTSourcedRasterBand::~VRTSourcedRasterBand()
113 :
114 : {
115 7649 : VRTSourcedRasterBand::CloseDependentDatasets();
116 7649 : CSLDestroy(m_papszSourceList);
117 15092 : }
118 :
119 : /************************************************************************/
120 : /* CanIRasterIOBeForwardedToEachSource() */
121 : /************************************************************************/
122 :
123 48389 : bool VRTSourcedRasterBand::CanIRasterIOBeForwardedToEachSource(
124 : GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
125 : int nBufXSize, int nBufYSize, GDALRasterIOExtraArg *psExtraArg) const
126 : {
127 32403 : const auto IsNonNearestInvolved = [this, psExtraArg]
128 : {
129 5175 : if (psExtraArg->eResampleAlg != GRIORA_NearestNeighbour)
130 : {
131 15 : return true;
132 : }
133 14478 : for (int i = 0; i < nSources; i++)
134 : {
135 9326 : if (papoSources[i]->GetType() == VRTComplexSource::GetTypeStatic())
136 : {
137 3424 : auto *const poSource =
138 3424 : static_cast<VRTComplexSource *>(papoSources[i]);
139 3424 : if (!poSource->GetResampling().empty())
140 8 : return true;
141 : }
142 : }
143 5152 : return false;
144 48389 : };
145 :
146 : // If resampling with non-nearest neighbour, we need to be careful
147 : // if the VRT band exposes a nodata value, but the sources do not have it.
148 : // To also avoid edge effects on sources when downsampling, use the
149 : // base implementation of IRasterIO() (that is acquiring sources at their
150 : // nominal resolution, and then downsampling), but only if none of the
151 : // contributing sources have overviews.
152 48389 : if (eRWFlag == GF_Read && (nXSize != nBufXSize || nYSize != nBufYSize) &&
153 96778 : nSources != 0 && IsNonNearestInvolved())
154 : {
155 23 : bool bSourceHasOverviews = false;
156 23 : const bool bIsDownsampling = (nBufXSize < nXSize && nBufYSize < nYSize);
157 23 : int nContributingSources = 0;
158 23 : bool bSourceFullySatisfiesRequest = true;
159 37 : for (int i = 0; i < nSources; i++)
160 : {
161 26 : if (!papoSources[i]->IsSimpleSource())
162 : {
163 0 : return false;
164 : }
165 : else
166 : {
167 26 : VRTSimpleSource *const poSource =
168 26 : static_cast<VRTSimpleSource *>(papoSources[i]);
169 :
170 26 : if (poSource->GetType() == VRTComplexSource::GetTypeStatic())
171 : {
172 16 : auto *const poComplexSource =
173 : static_cast<VRTComplexSource *>(poSource);
174 16 : if (!poComplexSource->GetResampling().empty())
175 : {
176 : const int lMaskFlags =
177 : const_cast<VRTSourcedRasterBand *>(this)
178 8 : ->GetMaskFlags();
179 6 : if ((lMaskFlags != GMF_ALL_VALID &&
180 16 : lMaskFlags != GMF_NODATA) ||
181 2 : IsMaskBand())
182 : {
183 : // Unfortunately this will prevent using overviews
184 : // of the sources, but it is unpractical to use
185 : // them without serious implementation complications
186 12 : return false;
187 : }
188 : }
189 : }
190 :
191 18 : double dfXOff = nXOff;
192 18 : double dfYOff = nYOff;
193 18 : double dfXSize = nXSize;
194 18 : double dfYSize = nYSize;
195 18 : if (psExtraArg->bFloatingPointWindowValidity)
196 : {
197 1 : dfXOff = psExtraArg->dfXOff;
198 1 : dfYOff = psExtraArg->dfYOff;
199 1 : dfXSize = psExtraArg->dfXSize;
200 1 : dfYSize = psExtraArg->dfYSize;
201 : }
202 :
203 : // The window we will actually request from the source raster
204 : // band.
205 18 : double dfReqXOff = 0.0;
206 18 : double dfReqYOff = 0.0;
207 18 : double dfReqXSize = 0.0;
208 18 : double dfReqYSize = 0.0;
209 18 : int nReqXOff = 0;
210 18 : int nReqYOff = 0;
211 18 : int nReqXSize = 0;
212 18 : int nReqYSize = 0;
213 :
214 : // The window we will actual set _within_ the pData buffer.
215 18 : int nOutXOff = 0;
216 18 : int nOutYOff = 0;
217 18 : int nOutXSize = 0;
218 18 : int nOutYSize = 0;
219 :
220 18 : bool bError = false;
221 18 : if (!poSource->GetSrcDstWindow(
222 : dfXOff, dfYOff, dfXSize, dfYSize, nBufXSize, nBufYSize,
223 : &dfReqXOff, &dfReqYOff, &dfReqXSize, &dfReqYSize,
224 : &nReqXOff, &nReqYOff, &nReqXSize, &nReqYSize, &nOutXOff,
225 : &nOutYOff, &nOutXSize, &nOutYSize, bError))
226 : {
227 0 : continue;
228 : }
229 18 : auto poBand = poSource->GetRasterBand();
230 18 : if (poBand == nullptr)
231 : {
232 0 : return false;
233 : }
234 18 : ++nContributingSources;
235 18 : if (!(nOutXOff == 0 && nOutYOff == 0 &&
236 17 : nOutXSize == nBufXSize && nOutYSize == nBufYSize))
237 4 : bSourceFullySatisfiesRequest = false;
238 18 : if (m_bNoDataValueSet)
239 : {
240 7 : int bSrcHasNoData = FALSE;
241 : const double dfSrcNoData =
242 7 : poBand->GetNoDataValue(&bSrcHasNoData);
243 7 : if (!bSrcHasNoData || dfSrcNoData != m_dfNoDataValue)
244 : {
245 4 : return false;
246 : }
247 : }
248 14 : if (bIsDownsampling)
249 : {
250 11 : if (poBand->GetOverviewCount() != 0)
251 : {
252 4 : bSourceHasOverviews = true;
253 : }
254 : }
255 : }
256 : }
257 11 : if (bIsDownsampling && !bSourceHasOverviews &&
258 3 : (nContributingSources > 1 || !bSourceFullySatisfiesRequest))
259 : {
260 2 : return false;
261 : }
262 : }
263 48375 : return true;
264 : }
265 :
266 : /************************************************************************/
267 : /* CanMultiThreadRasterIO() */
268 : /************************************************************************/
269 :
270 418 : bool VRTSourcedRasterBand::CanMultiThreadRasterIO(
271 : double dfXOff, double dfYOff, double dfXSize, double dfYSize,
272 : int &nContributingSources) const
273 : {
274 418 : int iLastSource = 0;
275 : CPLRectObj sSourceBounds;
276 418 : CPLQuadTree *hQuadTree = nullptr;
277 418 : bool bRet = true;
278 418 : std::set<std::string> oSetDSName;
279 :
280 418 : nContributingSources = 0;
281 1197 : for (int iSource = 0; iSource < nSources; iSource++)
282 : {
283 781 : const auto poSource = papoSources[iSource];
284 781 : if (!poSource->IsSimpleSource())
285 : {
286 0 : bRet = false;
287 0 : break;
288 : }
289 781 : const auto poSimpleSource = cpl::down_cast<VRTSimpleSource *>(poSource);
290 781 : if (poSimpleSource->DstWindowIntersects(dfXOff, dfYOff, dfXSize,
291 : dfYSize))
292 : {
293 : // Only build hQuadTree if there are 2 or more sources
294 516 : if (nContributingSources == 1)
295 : {
296 : std::string &oFirstSrcDSName =
297 26 : cpl::down_cast<VRTSimpleSource *>(papoSources[iLastSource])
298 26 : ->m_osSrcDSName;
299 26 : oSetDSName.insert(oFirstSrcDSName);
300 :
301 : CPLRectObj sGlobalBounds;
302 26 : sGlobalBounds.minx = dfXOff;
303 26 : sGlobalBounds.miny = dfYOff;
304 26 : sGlobalBounds.maxx = dfXOff + dfXSize;
305 26 : sGlobalBounds.maxy = dfYOff + dfYSize;
306 26 : hQuadTree = CPLQuadTreeCreate(&sGlobalBounds, nullptr);
307 :
308 26 : CPLQuadTreeInsertWithBounds(
309 : hQuadTree,
310 : reinterpret_cast<void *>(
311 26 : static_cast<uintptr_t>(iLastSource)),
312 : &sSourceBounds);
313 : }
314 :
315 : // Check there are not several sources with the same name, to avoid
316 : // the same GDALDataset* to be used from multiple threads. We may
317 : // be a bit too pessimistic, for example if working with unnamed
318 : // Memory datasets, but that would involve comparing
319 : // poSource->GetRasterBandNoOpen()->GetDataset()
320 516 : if (oSetDSName.find(poSimpleSource->m_osSrcDSName) !=
321 1032 : oSetDSName.end())
322 : {
323 0 : bRet = false;
324 2 : break;
325 : }
326 516 : oSetDSName.insert(poSimpleSource->m_osSrcDSName);
327 :
328 : double dfSourceXOff;
329 : double dfSourceYOff;
330 : double dfSourceXSize;
331 : double dfSourceYSize;
332 516 : poSimpleSource->GetDstWindow(dfSourceXOff, dfSourceYOff,
333 : dfSourceXSize, dfSourceYSize);
334 516 : constexpr double EPSILON = 1e-1;
335 516 : sSourceBounds.minx = dfSourceXOff + EPSILON;
336 516 : sSourceBounds.miny = dfSourceYOff + EPSILON;
337 516 : sSourceBounds.maxx = dfSourceXOff + dfSourceXSize - EPSILON;
338 516 : sSourceBounds.maxy = dfSourceYOff + dfSourceYSize - EPSILON;
339 516 : iLastSource = iSource;
340 :
341 516 : if (hQuadTree)
342 : {
343 : // Check that the new source doesn't overlap an existing one.
344 278 : if (CPLQuadTreeHasMatch(hQuadTree, &sSourceBounds))
345 : {
346 2 : bRet = false;
347 2 : break;
348 : }
349 :
350 276 : CPLQuadTreeInsertWithBounds(
351 : hQuadTree,
352 276 : reinterpret_cast<void *>(static_cast<uintptr_t>(iSource)),
353 : &sSourceBounds);
354 : }
355 :
356 514 : ++nContributingSources;
357 : }
358 : }
359 :
360 418 : if (hQuadTree)
361 26 : CPLQuadTreeDestroy(hQuadTree);
362 :
363 836 : return bRet;
364 : }
365 :
366 : /************************************************************************/
367 : /* IRasterIO() */
368 : /************************************************************************/
369 :
370 48391 : CPLErr VRTSourcedRasterBand::IRasterIO(
371 : GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
372 : void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
373 : GSpacing nPixelSpace, GSpacing nLineSpace, GDALRasterIOExtraArg *psExtraArg)
374 :
375 : {
376 48391 : if (eRWFlag == GF_Write)
377 : {
378 0 : CPLError(CE_Failure, CPLE_AppDefined,
379 : "Writing through VRTSourcedRasterBand is not supported.");
380 0 : return CE_Failure;
381 : }
382 :
383 96782 : const std::string osFctId("VRTSourcedRasterBand::IRasterIO");
384 96782 : GDALAntiRecursionGuard oGuard(osFctId);
385 48391 : if (oGuard.GetCallDepth() >= 32)
386 : {
387 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
388 0 : return CE_Failure;
389 : }
390 :
391 145173 : GDALAntiRecursionGuard oGuard2(oGuard, poDS->GetDescription());
392 : // Allow 2 recursion depths on the same dataset for non-nearest resampling
393 48391 : if (oGuard2.GetCallDepth() > 2)
394 : {
395 1 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
396 1 : return CE_Failure;
397 : }
398 :
399 : /* ==================================================================== */
400 : /* Do we have overviews that would be appropriate to satisfy */
401 : /* this request? */
402 : /* ==================================================================== */
403 48390 : auto l_poDS = dynamic_cast<VRTDataset *>(poDS);
404 48390 : if (l_poDS &&
405 96698 : l_poDS->m_apoOverviews.empty() && // do not use virtual overviews
406 96780 : (nBufXSize < nXSize || nBufYSize < nYSize) && GetOverviewCount() > 0)
407 : {
408 2 : if (OverviewRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize, pData,
409 : nBufXSize, nBufYSize, eBufType, nPixelSpace,
410 2 : nLineSpace, psExtraArg) == CE_None)
411 1 : return CE_None;
412 : }
413 :
414 : // If resampling with non-nearest neighbour, we need to be careful
415 : // if the VRT band exposes a nodata value, but the sources do not have it.
416 : // To also avoid edge effects on sources when downsampling, use the
417 : // base implementation of IRasterIO() (that is acquiring sources at their
418 : // nominal resolution, and then downsampling), but only if none of the
419 : // contributing sources have overviews.
420 48389 : if (l_poDS && !CanIRasterIOBeForwardedToEachSource(
421 : eRWFlag, nXOff, nYOff, nXSize, nYSize, nBufXSize,
422 : nBufYSize, psExtraArg))
423 : {
424 14 : const bool bBackupEnabledOverviews = l_poDS->AreOverviewsEnabled();
425 14 : if (!l_poDS->m_apoOverviews.empty() && l_poDS->AreOverviewsEnabled())
426 : {
427 : // Disable use of implicit overviews to avoid infinite
428 : // recursion
429 1 : l_poDS->SetEnableOverviews(false);
430 : }
431 :
432 14 : const auto eResampleAlgBackup = psExtraArg->eResampleAlg;
433 14 : if (psExtraArg->eResampleAlg == GRIORA_NearestNeighbour)
434 : {
435 16 : std::string osResampling;
436 24 : for (int i = 0; i < nSources; i++)
437 : {
438 16 : if (papoSources[i]->GetType() ==
439 16 : VRTComplexSource::GetTypeStatic())
440 : {
441 16 : auto *const poComplexSource =
442 16 : static_cast<VRTComplexSource *>(papoSources[i]);
443 16 : if (!poComplexSource->GetResampling().empty())
444 : {
445 16 : if (i == 0)
446 8 : osResampling = poComplexSource->GetResampling();
447 8 : else if (osResampling !=
448 8 : poComplexSource->GetResampling())
449 : {
450 0 : osResampling.clear();
451 0 : break;
452 : }
453 : }
454 : }
455 : }
456 8 : if (!osResampling.empty())
457 8 : psExtraArg->eResampleAlg =
458 8 : GDALRasterIOGetResampleAlg(osResampling.c_str());
459 : }
460 :
461 14 : const auto eErr = GDALRasterBand::IRasterIO(
462 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
463 : eBufType, nPixelSpace, nLineSpace, psExtraArg);
464 :
465 14 : psExtraArg->eResampleAlg = eResampleAlgBackup;
466 14 : l_poDS->SetEnableOverviews(bBackupEnabledOverviews);
467 14 : return eErr;
468 : }
469 :
470 : /* -------------------------------------------------------------------- */
471 : /* Initialize the buffer to some background value. Use the */
472 : /* nodata value if available. */
473 : /* -------------------------------------------------------------------- */
474 48375 : if (SkipBufferInitialization())
475 : {
476 : // Do nothing
477 : }
478 17741 : else if (nPixelSpace == GDALGetDataTypeSizeBytes(eBufType) &&
479 17530 : !(m_bNoDataValueSet && m_dfNoDataValue != 0.0) &&
480 52626 : !(m_bNoDataSetAsInt64 && m_nNoDataValueInt64 != 0) &&
481 17355 : !(m_bNoDataSetAsUInt64 && m_nNoDataValueUInt64 != 0))
482 : {
483 17354 : if (nLineSpace == nBufXSize * nPixelSpace)
484 : {
485 17276 : memset(pData, 0, static_cast<size_t>(nBufYSize * nLineSpace));
486 : }
487 : else
488 : {
489 10648 : for (int iLine = 0; iLine < nBufYSize; iLine++)
490 : {
491 10570 : memset(static_cast<GByte *>(pData) +
492 10570 : static_cast<GIntBig>(iLine) * nLineSpace,
493 10570 : 0, static_cast<size_t>(nBufXSize * nPixelSpace));
494 : }
495 : }
496 : }
497 387 : else if (m_bNoDataSetAsInt64)
498 : {
499 7 : for (int iLine = 0; iLine < nBufYSize; iLine++)
500 : {
501 6 : GDALCopyWords(&m_nNoDataValueInt64, GDT_Int64, 0,
502 : static_cast<GByte *>(pData) +
503 6 : static_cast<GIntBig>(nLineSpace) * iLine,
504 : eBufType, static_cast<int>(nPixelSpace), nBufXSize);
505 : }
506 : }
507 386 : else if (m_bNoDataSetAsUInt64)
508 : {
509 7 : for (int iLine = 0; iLine < nBufYSize; iLine++)
510 : {
511 6 : GDALCopyWords(&m_nNoDataValueUInt64, GDT_UInt64, 0,
512 : static_cast<GByte *>(pData) +
513 6 : static_cast<GIntBig>(nLineSpace) * iLine,
514 : eBufType, static_cast<int>(nPixelSpace), nBufXSize);
515 : }
516 : }
517 : else
518 : {
519 385 : double dfWriteValue = 0.0;
520 385 : if (m_bNoDataValueSet)
521 186 : dfWriteValue = m_dfNoDataValue;
522 :
523 47692 : for (int iLine = 0; iLine < nBufYSize; iLine++)
524 : {
525 47307 : GDALCopyWords(&dfWriteValue, GDT_Float64, 0,
526 : static_cast<GByte *>(pData) +
527 47307 : static_cast<GIntBig>(nLineSpace) * iLine,
528 : eBufType, static_cast<int>(nPixelSpace), nBufXSize);
529 : }
530 : }
531 :
532 : /* -------------------------------------------------------------------- */
533 : /* Overlay each source in turn over top this. */
534 : /* -------------------------------------------------------------------- */
535 48375 : CPLErr eErr = CE_None;
536 :
537 48375 : double dfXOff = nXOff;
538 48375 : double dfYOff = nYOff;
539 48375 : double dfXSize = nXSize;
540 48375 : double dfYSize = nYSize;
541 48375 : if (psExtraArg->bFloatingPointWindowValidity)
542 : {
543 248 : dfXOff = psExtraArg->dfXOff;
544 248 : dfYOff = psExtraArg->dfYOff;
545 248 : dfXSize = psExtraArg->dfXSize;
546 248 : dfYSize = psExtraArg->dfYSize;
547 : }
548 :
549 48375 : if (l_poDS)
550 48375 : l_poDS->m_bMultiThreadedRasterIOLastUsed = false;
551 :
552 48375 : int nContributingSources = 0;
553 48375 : int nMaxThreads = 0;
554 48375 : constexpr int MINIMUM_PIXEL_COUNT_FOR_THREADED_IO = 1000 * 1000;
555 96750 : if (l_poDS &&
556 48375 : (static_cast<int64_t>(nBufXSize) * nBufYSize >=
557 48066 : MINIMUM_PIXEL_COUNT_FOR_THREADED_IO ||
558 48066 : static_cast<int64_t>(nXSize) * nYSize >=
559 316 : MINIMUM_PIXEL_COUNT_FOR_THREADED_IO) &&
560 316 : CanMultiThreadRasterIO(dfXOff, dfYOff, dfXSize, dfYSize,
561 315 : nContributingSources) &&
562 96770 : nContributingSources > 1 &&
563 20 : (nMaxThreads = VRTDataset::GetNumThreads(l_poDS)) > 1)
564 : {
565 18 : l_poDS->m_bMultiThreadedRasterIOLastUsed = true;
566 18 : l_poDS->m_oMapSharedSources.InitMutex();
567 :
568 18 : std::atomic<bool> bSuccess = true;
569 18 : CPLWorkerThreadPool *psThreadPool = GDALGetGlobalThreadPool(
570 18 : std::min(nContributingSources, nMaxThreads));
571 : const int nThreads =
572 18 : std::min(nContributingSources, psThreadPool->GetThreadCount());
573 18 : CPLDebugOnly("VRT",
574 : "IRasterIO(): use optimized "
575 : "multi-threaded code path for mosaic. "
576 : "Using %d threads",
577 : nThreads);
578 :
579 : {
580 36 : std::lock_guard oLock(l_poDS->m_oQueueWorkingStates.oMutex);
581 18 : if (l_poDS->m_oQueueWorkingStates.oStates.size() <
582 18 : static_cast<size_t>(nThreads))
583 : {
584 6 : l_poDS->m_oQueueWorkingStates.oStates.resize(nThreads);
585 : }
586 117 : for (int i = 0; i < nThreads; ++i)
587 : {
588 99 : if (!l_poDS->m_oQueueWorkingStates.oStates[i])
589 75 : l_poDS->m_oQueueWorkingStates.oStates[i] =
590 150 : std::make_unique<VRTSource::WorkingState>();
591 : }
592 : }
593 :
594 18 : auto oQueue = psThreadPool->CreateJobQueue();
595 18 : std::atomic<int> nCompletedJobs = 0;
596 180 : for (int iSource = 0; iSource < nSources; iSource++)
597 : {
598 162 : auto poSource = papoSources[iSource];
599 162 : if (!poSource->IsSimpleSource())
600 0 : continue;
601 162 : auto poSimpleSource = cpl::down_cast<VRTSimpleSource *>(poSource);
602 162 : if (poSimpleSource->DstWindowIntersects(dfXOff, dfYOff, dfXSize,
603 : dfYSize))
604 : {
605 99 : auto psJob = new RasterIOJob();
606 99 : psJob->pbSuccess = &bSuccess;
607 99 : psJob->pnCompletedJobs = &nCompletedJobs;
608 99 : psJob->poQueueWorkingStates = &(l_poDS->m_oQueueWorkingStates);
609 99 : psJob->eVRTBandDataType = eDataType;
610 99 : psJob->nXOff = nXOff;
611 99 : psJob->nYOff = nYOff;
612 99 : psJob->nXSize = nXSize;
613 99 : psJob->nYSize = nYSize;
614 99 : psJob->pData = pData;
615 99 : psJob->nBufXSize = nBufXSize;
616 99 : psJob->nBufYSize = nBufYSize;
617 99 : psJob->eBufType = eBufType;
618 99 : psJob->nPixelSpace = nPixelSpace;
619 99 : psJob->nLineSpace = nLineSpace;
620 99 : psJob->psExtraArg = psExtraArg;
621 99 : psJob->poSource = poSimpleSource;
622 :
623 99 : if (!oQueue->SubmitJob(RasterIOJob::Func, psJob))
624 : {
625 0 : delete psJob;
626 0 : bSuccess = false;
627 0 : break;
628 : }
629 : }
630 : }
631 :
632 80 : while (oQueue->WaitEvent())
633 : {
634 : // Quite rough progress callback. We could do better by counting
635 : // the number of contributing pixels.
636 62 : if (psExtraArg->pfnProgress)
637 : {
638 102 : psExtraArg->pfnProgress(double(nCompletedJobs.load()) /
639 : nContributingSources,
640 : "", psExtraArg->pProgressData);
641 : }
642 : }
643 :
644 18 : eErr = bSuccess ? CE_None : CE_Failure;
645 : }
646 : else
647 : {
648 48357 : GDALProgressFunc const pfnProgressGlobal = psExtraArg->pfnProgress;
649 48357 : void *const pProgressDataGlobal = psExtraArg->pProgressData;
650 :
651 48357 : VRTSource::WorkingState oWorkingState;
652 114033 : for (int iSource = 0; eErr == CE_None && iSource < nSources; iSource++)
653 : {
654 65676 : psExtraArg->pfnProgress = GDALScaledProgress;
655 131352 : psExtraArg->pProgressData = GDALCreateScaledProgress(
656 65676 : 1.0 * iSource / nSources, 1.0 * (iSource + 1) / nSources,
657 : pfnProgressGlobal, pProgressDataGlobal);
658 65676 : if (psExtraArg->pProgressData == nullptr)
659 62600 : psExtraArg->pfnProgress = nullptr;
660 :
661 65676 : eErr = papoSources[iSource]->RasterIO(
662 : eDataType, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize,
663 : nBufYSize, eBufType, nPixelSpace, nLineSpace, psExtraArg,
664 65676 : l_poDS ? l_poDS->m_oWorkingState : oWorkingState);
665 :
666 65676 : GDALDestroyScaledProgress(psExtraArg->pProgressData);
667 : }
668 :
669 48357 : psExtraArg->pfnProgress = pfnProgressGlobal;
670 48357 : psExtraArg->pProgressData = pProgressDataGlobal;
671 : }
672 :
673 48375 : if (eErr == CE_None && psExtraArg->pfnProgress)
674 : {
675 2943 : psExtraArg->pfnProgress(1.0, "", psExtraArg->pProgressData);
676 : }
677 :
678 48375 : return eErr;
679 : }
680 :
681 : /************************************************************************/
682 : /* VRTSourcedRasterBand::RasterIOJob::Func() */
683 : /************************************************************************/
684 :
685 99 : void VRTSourcedRasterBand::RasterIOJob::Func(void *pData)
686 : {
687 : auto psJob =
688 198 : std::unique_ptr<RasterIOJob>(static_cast<RasterIOJob *>(pData));
689 99 : if (*psJob->pbSuccess)
690 : {
691 99 : GDALRasterIOExtraArg sArg = *(psJob->psExtraArg);
692 99 : sArg.pfnProgress = nullptr;
693 99 : sArg.pProgressData = nullptr;
694 :
695 99 : std::unique_ptr<VRTSource::WorkingState> poWorkingState;
696 : {
697 198 : std::lock_guard oLock(psJob->poQueueWorkingStates->oMutex);
698 : poWorkingState =
699 99 : std::move(psJob->poQueueWorkingStates->oStates.back());
700 99 : psJob->poQueueWorkingStates->oStates.pop_back();
701 99 : CPLAssert(poWorkingState.get());
702 : }
703 :
704 198 : if (psJob->poSource->RasterIO(
705 99 : psJob->eVRTBandDataType, psJob->nXOff, psJob->nYOff,
706 99 : psJob->nXSize, psJob->nYSize, psJob->pData, psJob->nBufXSize,
707 99 : psJob->nBufYSize, psJob->eBufType, psJob->nPixelSpace,
708 198 : psJob->nLineSpace, &sArg, *(poWorkingState.get())) != CE_None)
709 : {
710 0 : *psJob->pbSuccess = false;
711 : }
712 :
713 : {
714 198 : std::lock_guard oLock(psJob->poQueueWorkingStates->oMutex);
715 198 : psJob->poQueueWorkingStates->oStates.push_back(
716 99 : std::move(poWorkingState));
717 : }
718 : }
719 :
720 99 : ++(*psJob->pnCompletedJobs);
721 99 : }
722 :
723 : /************************************************************************/
724 : /* IGetDataCoverageStatus() */
725 : /************************************************************************/
726 :
727 1336 : int VRTSourcedRasterBand::IGetDataCoverageStatus(int nXOff, int nYOff,
728 : int nXSize, int nYSize,
729 : int nMaskFlagStop,
730 : double *pdfDataPct)
731 : {
732 1336 : if (pdfDataPct)
733 8 : *pdfDataPct = -1.0;
734 :
735 : // Particular case for a single simple source covering the whole dataset
736 2657 : if (nSources == 1 && papoSources[0]->IsSimpleSource() &&
737 1321 : papoSources[0]->GetType() == VRTSimpleSource::GetTypeStatic())
738 : {
739 1294 : VRTSimpleSource *poSource =
740 1294 : static_cast<VRTSimpleSource *>(papoSources[0]);
741 :
742 1294 : GDALRasterBand *poBand = poSource->GetRasterBand();
743 1294 : if (!poBand)
744 0 : poBand = poSource->GetMaskBandMainBand();
745 1294 : if (!poBand)
746 : {
747 : return GDAL_DATA_COVERAGE_STATUS_UNIMPLEMENTED |
748 301 : GDAL_DATA_COVERAGE_STATUS_DATA;
749 : }
750 :
751 : /* Check that it uses the full source dataset */
752 1294 : double dfReqXOff = 0.0;
753 1294 : double dfReqYOff = 0.0;
754 1294 : double dfReqXSize = 0.0;
755 1294 : double dfReqYSize = 0.0;
756 1294 : int nReqXOff = 0;
757 1294 : int nReqYOff = 0;
758 1294 : int nReqXSize = 0;
759 1294 : int nReqYSize = 0;
760 1294 : int nOutXOff = 0;
761 1294 : int nOutYOff = 0;
762 1294 : int nOutXSize = 0;
763 1294 : int nOutYSize = 0;
764 1294 : bool bError = false;
765 1294 : if (poSource->GetSrcDstWindow(
766 1294 : 0, 0, GetXSize(), GetYSize(), GetXSize(), GetYSize(),
767 : &dfReqXOff, &dfReqYOff, &dfReqXSize, &dfReqYSize, &nReqXOff,
768 : &nReqYOff, &nReqXSize, &nReqYSize, &nOutXOff, &nOutYOff,
769 1294 : &nOutXSize, &nOutYSize, bError) &&
770 1294 : nReqXOff == 0 && nReqYOff == 0 && nReqXSize == GetXSize() &&
771 371 : nReqXSize == poBand->GetXSize() && nReqYSize == GetYSize() &&
772 302 : nReqYSize == poBand->GetYSize() && nOutXOff == 0 && nOutYOff == 0 &&
773 2588 : nOutXSize == GetXSize() && nOutYSize == GetYSize())
774 : {
775 301 : return poBand->GetDataCoverageStatus(nXOff, nYOff, nXSize, nYSize,
776 301 : nMaskFlagStop, pdfDataPct);
777 : }
778 : }
779 :
780 : #ifndef HAVE_GEOS
781 : return GDAL_DATA_COVERAGE_STATUS_UNIMPLEMENTED |
782 : GDAL_DATA_COVERAGE_STATUS_DATA;
783 : #else
784 1035 : int nStatus = 0;
785 :
786 2070 : auto poPolyNonCoveredBySources = std::make_unique<OGRPolygon>();
787 : {
788 2070 : auto poLR = std::make_unique<OGRLinearRing>();
789 1035 : poLR->addPoint(nXOff, nYOff);
790 1035 : poLR->addPoint(nXOff, nYOff + nYSize);
791 1035 : poLR->addPoint(nXOff + nXSize, nYOff + nYSize);
792 1035 : poLR->addPoint(nXOff + nXSize, nYOff);
793 1035 : poLR->addPoint(nXOff, nYOff);
794 1035 : poPolyNonCoveredBySources->addRingDirectly(poLR.release());
795 : }
796 :
797 1045 : for (int iSource = 0; iSource < nSources; iSource++)
798 : {
799 1040 : if (!papoSources[iSource]->IsSimpleSource())
800 : {
801 : return GDAL_DATA_COVERAGE_STATUS_UNIMPLEMENTED |
802 4 : GDAL_DATA_COVERAGE_STATUS_DATA;
803 : }
804 1036 : VRTSimpleSource *poSS =
805 1036 : static_cast<VRTSimpleSource *>(papoSources[iSource]);
806 : // Check if the AOI is fully inside the source
807 1036 : double dfDstXOff = std::max(0.0, poSS->m_dfDstXOff);
808 1036 : double dfDstYOff = std::max(0.0, poSS->m_dfDstYOff);
809 1036 : double dfDstXSize = poSS->m_dfDstXSize;
810 1036 : double dfDstYSize = poSS->m_dfDstYSize;
811 1036 : auto l_poBand = poSS->GetRasterBand();
812 1036 : if (!l_poBand)
813 0 : continue;
814 1036 : if (dfDstXSize == -1)
815 8 : dfDstXSize = l_poBand->GetXSize() - dfDstXOff;
816 1036 : if (dfDstYSize == -1)
817 8 : dfDstYSize = l_poBand->GetYSize() - dfDstYOff;
818 :
819 1036 : if (nXOff >= dfDstXOff && nYOff >= dfDstYOff &&
820 1027 : nXOff + nXSize <= dfDstXOff + dfDstXSize &&
821 1011 : nYOff + nYSize <= dfDstYOff + dfDstYSize)
822 : {
823 1010 : if (pdfDataPct)
824 1 : *pdfDataPct = 100.0;
825 1010 : return GDAL_DATA_COVERAGE_STATUS_DATA;
826 : }
827 : // Check intersection of bounding boxes.
828 26 : if (dfDstXOff + dfDstXSize > nXOff && dfDstYOff + dfDstYSize > nYOff &&
829 20 : dfDstXOff < nXOff + nXSize && dfDstYOff < nYOff + nYSize)
830 : {
831 18 : nStatus |= GDAL_DATA_COVERAGE_STATUS_DATA;
832 18 : if (poPolyNonCoveredBySources)
833 : {
834 18 : OGRPolygon oPolySource;
835 18 : auto poLR = std::make_unique<OGRLinearRing>();
836 18 : poLR->addPoint(dfDstXOff, dfDstYOff);
837 18 : poLR->addPoint(dfDstXOff, dfDstYOff + dfDstYSize);
838 18 : poLR->addPoint(dfDstXOff + dfDstXSize, dfDstYOff + dfDstYSize);
839 18 : poLR->addPoint(dfDstXOff + dfDstXSize, dfDstYOff);
840 18 : poLR->addPoint(dfDstXOff, dfDstYOff);
841 18 : oPolySource.addRingDirectly(poLR.release());
842 : auto poRes = std::unique_ptr<OGRGeometry>(
843 18 : poPolyNonCoveredBySources->Difference(&oPolySource));
844 18 : if (poRes && poRes->IsEmpty())
845 : {
846 1 : if (pdfDataPct)
847 1 : *pdfDataPct = 100.0;
848 1 : return GDAL_DATA_COVERAGE_STATUS_DATA;
849 : }
850 17 : else if (poRes && poRes->getGeometryType() == wkbPolygon)
851 : {
852 17 : poPolyNonCoveredBySources.reset(
853 : poRes.release()->toPolygon());
854 : }
855 : else
856 : {
857 0 : poPolyNonCoveredBySources.reset();
858 : }
859 : }
860 : }
861 25 : if (nMaskFlagStop != 0 && (nStatus & nMaskFlagStop) != 0)
862 : {
863 15 : return nStatus;
864 : }
865 : }
866 5 : if (poPolyNonCoveredBySources)
867 : {
868 5 : if (!poPolyNonCoveredBySources->IsEmpty())
869 5 : nStatus |= GDAL_DATA_COVERAGE_STATUS_EMPTY;
870 5 : if (pdfDataPct)
871 2 : *pdfDataPct = 100.0 * (1.0 - poPolyNonCoveredBySources->get_Area() /
872 2 : nXSize / nYSize);
873 : }
874 5 : return nStatus;
875 : #endif // HAVE_GEOS
876 : }
877 :
878 : /************************************************************************/
879 : /* IReadBlock() */
880 : /************************************************************************/
881 :
882 668 : CPLErr VRTSourcedRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
883 : void *pImage)
884 :
885 : {
886 668 : const int nPixelSize = GDALGetDataTypeSizeBytes(eDataType);
887 :
888 668 : int nReadXSize = 0;
889 668 : if ((nBlockXOff + 1) * nBlockXSize > GetXSize())
890 0 : nReadXSize = GetXSize() - nBlockXOff * nBlockXSize;
891 : else
892 668 : nReadXSize = nBlockXSize;
893 :
894 668 : int nReadYSize = 0;
895 668 : if ((nBlockYOff + 1) * nBlockYSize > GetYSize())
896 0 : nReadYSize = GetYSize() - nBlockYOff * nBlockYSize;
897 : else
898 668 : nReadYSize = nBlockYSize;
899 :
900 : GDALRasterIOExtraArg sExtraArg;
901 668 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
902 :
903 668 : return IRasterIO(
904 668 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize, nReadXSize,
905 : nReadYSize, pImage, nReadXSize, nReadYSize, eDataType, nPixelSize,
906 1336 : static_cast<GSpacing>(nPixelSize) * nBlockXSize, &sExtraArg);
907 : }
908 :
909 : /************************************************************************/
910 : /* CPLGettimeofday() */
911 : /************************************************************************/
912 :
913 : #if defined(_WIN32) && !defined(__CYGWIN__)
914 : #include <sys/timeb.h>
915 :
916 : namespace
917 : {
918 : struct CPLTimeVal
919 : {
920 : time_t tv_sec; /* seconds */
921 : long tv_usec; /* and microseconds */
922 : };
923 : } // namespace
924 :
925 : static int CPLGettimeofday(struct CPLTimeVal *tp, void * /* timezonep*/)
926 : {
927 : struct _timeb theTime;
928 :
929 : _ftime(&theTime);
930 : tp->tv_sec = static_cast<time_t>(theTime.time);
931 : tp->tv_usec = theTime.millitm * 1000;
932 : return 0;
933 : }
934 : #else
935 : #include <sys/time.h> /* for gettimeofday() */
936 : #define CPLTimeVal timeval
937 : #define CPLGettimeofday(t, u) gettimeofday(t, u)
938 : #endif
939 :
940 : /************************************************************************/
941 : /* CanUseSourcesMinMaxImplementations() */
942 : /************************************************************************/
943 :
944 107 : bool VRTSourcedRasterBand::CanUseSourcesMinMaxImplementations()
945 : {
946 : const char *pszUseSources =
947 107 : CPLGetConfigOption("VRT_MIN_MAX_FROM_SOURCES", nullptr);
948 107 : if (pszUseSources)
949 0 : return CPLTestBool(pszUseSources);
950 :
951 : // Use heuristics to determine if we are going to use the source
952 : // GetMinimum() or GetMaximum() implementation: all the sources must be
953 : // "simple" sources with a dataset description that match a "regular" file
954 : // on the filesystem, whose open time and GetMinimum()/GetMaximum()
955 : // implementations we hope to be fast enough.
956 : // In case of doubt return FALSE.
957 : struct CPLTimeVal tvStart;
958 107 : memset(&tvStart, 0, sizeof(CPLTimeVal));
959 107 : if (nSources > 1)
960 10 : CPLGettimeofday(&tvStart, nullptr);
961 235 : for (int iSource = 0; iSource < nSources; iSource++)
962 : {
963 131 : if (!(papoSources[iSource]->IsSimpleSource()))
964 2 : return false;
965 129 : VRTSimpleSource *const poSimpleSource =
966 129 : static_cast<VRTSimpleSource *>(papoSources[iSource]);
967 129 : const char *pszFilename = poSimpleSource->m_osSrcDSName.c_str();
968 : // /vsimem/ should be fast.
969 129 : if (STARTS_WITH(pszFilename, "/vsimem/"))
970 10 : continue;
971 : // but not other /vsi filesystems
972 119 : if (STARTS_WITH(pszFilename, "/vsi"))
973 0 : return false;
974 119 : char ch = '\0';
975 : // We will assume that filenames that are only with ascii characters
976 : // are real filenames and so we will not try to 'stat' them.
977 1815 : for (int i = 0; (ch = pszFilename[i]) != '\0'; i++)
978 : {
979 1718 : if (!((ch >= 'a' && ch <= 'z') || (ch >= 'A' && ch <= 'Z') ||
980 409 : (ch >= '0' && ch <= '9') || ch == ':' || ch == '/' ||
981 143 : ch == '\\' || ch == ' ' || ch == '.' || ch == '_'))
982 5 : break;
983 : }
984 119 : if (ch != '\0')
985 : {
986 : // Otherwise do a real filesystem check.
987 : VSIStatBuf sStat;
988 5 : if (VSIStat(pszFilename, &sStat) != 0)
989 1 : return false;
990 4 : if (nSources > 1)
991 : {
992 : struct CPLTimeVal tvCur;
993 4 : CPLGettimeofday(&tvCur, nullptr);
994 4 : if (tvCur.tv_sec - tvStart.tv_sec +
995 4 : (tvCur.tv_usec - tvStart.tv_usec) * 1e-6 >
996 : 1)
997 0 : return false;
998 : }
999 : }
1000 : }
1001 104 : return true;
1002 : }
1003 :
1004 : /************************************************************************/
1005 : /* GetMinimum() */
1006 : /************************************************************************/
1007 :
1008 59 : double VRTSourcedRasterBand::GetMinimum(int *pbSuccess)
1009 : {
1010 59 : const char *const pszValue = GetMetadataItem("STATISTICS_MINIMUM");
1011 59 : if (pszValue != nullptr)
1012 : {
1013 5 : if (pbSuccess != nullptr)
1014 5 : *pbSuccess = TRUE;
1015 :
1016 5 : return CPLAtofM(pszValue);
1017 : }
1018 :
1019 54 : if (!CanUseSourcesMinMaxImplementations())
1020 2 : return GDALRasterBand::GetMinimum(pbSuccess);
1021 :
1022 104 : const std::string osFctId("VRTSourcedRasterBand::GetMinimum");
1023 104 : GDALAntiRecursionGuard oGuard(osFctId);
1024 52 : if (oGuard.GetCallDepth() >= 32)
1025 : {
1026 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
1027 0 : if (pbSuccess != nullptr)
1028 0 : *pbSuccess = FALSE;
1029 0 : return 0;
1030 : }
1031 :
1032 156 : GDALAntiRecursionGuard oGuard2(oGuard, poDS->GetDescription());
1033 52 : if (oGuard2.GetCallDepth() >= 2)
1034 : {
1035 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
1036 0 : if (pbSuccess != nullptr)
1037 0 : *pbSuccess = FALSE;
1038 0 : return 0;
1039 : }
1040 :
1041 : struct CPLTimeVal tvStart;
1042 52 : memset(&tvStart, 0, sizeof(CPLTimeVal));
1043 52 : if (nSources > 1)
1044 5 : CPLGettimeofday(&tvStart, nullptr);
1045 52 : double dfMin = 0;
1046 57 : for (int iSource = 0; iSource < nSources; iSource++)
1047 : {
1048 53 : int bSuccess = FALSE;
1049 : double dfSourceMin =
1050 53 : papoSources[iSource]->GetMinimum(GetXSize(), GetYSize(), &bSuccess);
1051 53 : if (!bSuccess)
1052 : {
1053 48 : dfMin = GDALRasterBand::GetMinimum(pbSuccess);
1054 48 : return dfMin;
1055 : }
1056 :
1057 5 : if (iSource == 0 || dfSourceMin < dfMin)
1058 : {
1059 4 : dfMin = dfSourceMin;
1060 4 : if (dfMin == 0 && eDataType == GDT_Byte)
1061 0 : break;
1062 : }
1063 5 : if (nSources > 1)
1064 : {
1065 : struct CPLTimeVal tvCur;
1066 2 : CPLGettimeofday(&tvCur, nullptr);
1067 2 : if (tvCur.tv_sec - tvStart.tv_sec +
1068 2 : (tvCur.tv_usec - tvStart.tv_usec) * 1e-6 >
1069 : 1)
1070 : {
1071 0 : return GDALRasterBand::GetMinimum(pbSuccess);
1072 : }
1073 : }
1074 : }
1075 :
1076 4 : if (pbSuccess != nullptr)
1077 4 : *pbSuccess = TRUE;
1078 :
1079 4 : return dfMin;
1080 : }
1081 :
1082 : /************************************************************************/
1083 : /* GetMaximum() */
1084 : /************************************************************************/
1085 :
1086 57 : double VRTSourcedRasterBand::GetMaximum(int *pbSuccess)
1087 : {
1088 57 : const char *const pszValue = GetMetadataItem("STATISTICS_MAXIMUM");
1089 57 : if (pszValue != nullptr)
1090 : {
1091 4 : if (pbSuccess != nullptr)
1092 4 : *pbSuccess = TRUE;
1093 :
1094 4 : return CPLAtofM(pszValue);
1095 : }
1096 :
1097 53 : if (!CanUseSourcesMinMaxImplementations())
1098 1 : return GDALRasterBand::GetMaximum(pbSuccess);
1099 :
1100 104 : const std::string osFctId("VRTSourcedRasterBand::GetMaximum");
1101 104 : GDALAntiRecursionGuard oGuard(osFctId);
1102 52 : if (oGuard.GetCallDepth() >= 32)
1103 : {
1104 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
1105 0 : if (pbSuccess != nullptr)
1106 0 : *pbSuccess = FALSE;
1107 0 : return 0;
1108 : }
1109 :
1110 156 : GDALAntiRecursionGuard oGuard2(oGuard, poDS->GetDescription());
1111 52 : if (oGuard2.GetCallDepth() >= 2)
1112 : {
1113 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
1114 0 : if (pbSuccess != nullptr)
1115 0 : *pbSuccess = FALSE;
1116 0 : return 0;
1117 : }
1118 :
1119 : struct CPLTimeVal tvStart;
1120 52 : memset(&tvStart, 0, sizeof(CPLTimeVal));
1121 52 : if (nSources > 1)
1122 5 : CPLGettimeofday(&tvStart, nullptr);
1123 52 : double dfMax = 0;
1124 54 : for (int iSource = 0; iSource < nSources; iSource++)
1125 : {
1126 52 : int bSuccess = FALSE;
1127 : const double dfSourceMax =
1128 52 : papoSources[iSource]->GetMaximum(GetXSize(), GetYSize(), &bSuccess);
1129 52 : if (!bSuccess)
1130 : {
1131 48 : dfMax = GDALRasterBand::GetMaximum(pbSuccess);
1132 48 : return dfMax;
1133 : }
1134 :
1135 4 : if (iSource == 0 || dfSourceMax > dfMax)
1136 : {
1137 4 : dfMax = dfSourceMax;
1138 4 : if (dfMax == 255.0 && eDataType == GDT_Byte)
1139 2 : break;
1140 : }
1141 2 : if (nSources > 1)
1142 : {
1143 : struct CPLTimeVal tvCur;
1144 0 : CPLGettimeofday(&tvCur, nullptr);
1145 0 : if (tvCur.tv_sec - tvStart.tv_sec +
1146 0 : (tvCur.tv_usec - tvStart.tv_usec) * 1e-6 >
1147 : 1)
1148 : {
1149 0 : return GDALRasterBand::GetMaximum(pbSuccess);
1150 : }
1151 : }
1152 : }
1153 :
1154 4 : if (pbSuccess != nullptr)
1155 4 : *pbSuccess = TRUE;
1156 :
1157 4 : return dfMax;
1158 : }
1159 :
1160 : /************************************************************************/
1161 : /* IsMosaicOfNonOverlappingSimpleSourcesOfFullRasterNoResAndTypeChange() */
1162 : /************************************************************************/
1163 :
1164 : /* Returns true if the VRT raster band consists of non-overlapping simple
1165 : * sources or complex sources that don't change values, and use the full extent
1166 : * of the source band.
1167 : */
1168 49 : bool VRTSourcedRasterBand::
1169 : IsMosaicOfNonOverlappingSimpleSourcesOfFullRasterNoResAndTypeChange(
1170 : bool bAllowMaxValAdjustment) const
1171 : {
1172 49 : bool bRet = true;
1173 : CPLRectObj sGlobalBounds;
1174 49 : sGlobalBounds.minx = 0;
1175 49 : sGlobalBounds.miny = 0;
1176 49 : sGlobalBounds.maxx = nRasterXSize;
1177 49 : sGlobalBounds.maxy = nRasterYSize;
1178 49 : CPLQuadTree *hQuadTree = CPLQuadTreeCreate(&sGlobalBounds, nullptr);
1179 114 : for (int i = 0; i < nSources; ++i)
1180 : {
1181 74 : if (!papoSources[i]->IsSimpleSource())
1182 : {
1183 1 : bRet = false;
1184 9 : break;
1185 : }
1186 :
1187 73 : auto poSimpleSource = cpl::down_cast<VRTSimpleSource *>(papoSources[i]);
1188 73 : const char *pszType = poSimpleSource->GetType();
1189 73 : if (pszType == VRTSimpleSource::GetTypeStatic())
1190 : {
1191 : // ok
1192 : }
1193 5 : else if (pszType == VRTComplexSource::GetTypeStatic())
1194 : {
1195 : auto poComplexSource =
1196 5 : cpl::down_cast<VRTComplexSource *>(papoSources[i]);
1197 5 : if (!poComplexSource->AreValuesUnchanged())
1198 : {
1199 2 : bRet = false;
1200 2 : break;
1201 : }
1202 : }
1203 : else
1204 : {
1205 0 : bRet = false;
1206 0 : break;
1207 : }
1208 :
1209 71 : if (!bAllowMaxValAdjustment && poSimpleSource->NeedMaxValAdjustment())
1210 : {
1211 2 : bRet = false;
1212 2 : break;
1213 : }
1214 :
1215 69 : auto poSimpleSourceBand = poSimpleSource->GetRasterBand();
1216 138 : if (poSimpleSourceBand == nullptr ||
1217 69 : poSimpleSourceBand->GetRasterDataType() != eDataType)
1218 : {
1219 0 : bRet = false;
1220 0 : break;
1221 : }
1222 :
1223 69 : double dfReqXOff = 0.0;
1224 69 : double dfReqYOff = 0.0;
1225 69 : double dfReqXSize = 0.0;
1226 69 : double dfReqYSize = 0.0;
1227 69 : int nReqXOff = 0;
1228 69 : int nReqYOff = 0;
1229 69 : int nReqXSize = 0;
1230 69 : int nReqYSize = 0;
1231 69 : int nOutXOff = 0;
1232 69 : int nOutYOff = 0;
1233 69 : int nOutXSize = 0;
1234 69 : int nOutYSize = 0;
1235 :
1236 69 : bool bError = false;
1237 207 : if (!poSimpleSource->GetSrcDstWindow(
1238 69 : 0, 0, nRasterXSize, nRasterYSize, nRasterXSize, nRasterYSize,
1239 : &dfReqXOff, &dfReqYOff, &dfReqXSize, &dfReqYSize, &nReqXOff,
1240 : &nReqYOff, &nReqXSize, &nReqYSize, &nOutXOff, &nOutYOff,
1241 69 : &nOutXSize, &nOutYSize, bError) ||
1242 69 : nReqXOff != 0 || nReqYOff != 0 ||
1243 68 : nReqXSize != poSimpleSourceBand->GetXSize() ||
1244 68 : nReqYSize != poSimpleSourceBand->GetYSize() ||
1245 138 : nOutXSize != nReqXSize || nOutYSize != nReqYSize)
1246 : {
1247 1 : bRet = false;
1248 1 : break;
1249 : }
1250 :
1251 : CPLRectObj sBounds;
1252 68 : constexpr double EPSILON = 1e-1;
1253 68 : sBounds.minx = nOutXOff + EPSILON;
1254 68 : sBounds.miny = nOutYOff + EPSILON;
1255 68 : sBounds.maxx = nOutXOff + nOutXSize - EPSILON;
1256 68 : sBounds.maxy = nOutYOff + nOutYSize - EPSILON;
1257 :
1258 : // Check that the new source doesn't overlap an existing one.
1259 68 : if (CPLQuadTreeHasMatch(hQuadTree, &sBounds))
1260 : {
1261 3 : bRet = false;
1262 3 : break;
1263 : }
1264 :
1265 65 : CPLQuadTreeInsertWithBounds(
1266 65 : hQuadTree, reinterpret_cast<void *>(static_cast<uintptr_t>(i)),
1267 : &sBounds);
1268 : }
1269 49 : CPLQuadTreeDestroy(hQuadTree);
1270 :
1271 49 : return bRet;
1272 : }
1273 :
1274 : /************************************************************************/
1275 : /* ComputeRasterMinMax() */
1276 : /************************************************************************/
1277 :
1278 22 : CPLErr VRTSourcedRasterBand::ComputeRasterMinMax(int bApproxOK,
1279 : double *adfMinMax)
1280 : {
1281 : /* -------------------------------------------------------------------- */
1282 : /* Does the driver already know the min/max? */
1283 : /* -------------------------------------------------------------------- */
1284 22 : if (bApproxOK)
1285 : {
1286 3 : int bSuccessMin = FALSE;
1287 3 : int bSuccessMax = FALSE;
1288 :
1289 3 : const double dfMin = GetMinimum(&bSuccessMin);
1290 3 : const double dfMax = GetMaximum(&bSuccessMax);
1291 :
1292 3 : if (bSuccessMin && bSuccessMax)
1293 : {
1294 1 : adfMinMax[0] = dfMin;
1295 1 : adfMinMax[1] = dfMax;
1296 1 : return CE_None;
1297 : }
1298 : }
1299 :
1300 42 : const std::string osFctId("VRTSourcedRasterBand::ComputeRasterMinMax");
1301 42 : GDALAntiRecursionGuard oGuard(osFctId);
1302 21 : if (oGuard.GetCallDepth() >= 32)
1303 : {
1304 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
1305 0 : return CE_Failure;
1306 : }
1307 :
1308 63 : GDALAntiRecursionGuard oGuard2(oGuard, poDS->GetDescription());
1309 21 : if (oGuard2.GetCallDepth() >= 2)
1310 : {
1311 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
1312 0 : return CE_Failure;
1313 : }
1314 :
1315 : /* -------------------------------------------------------------------- */
1316 : /* If we have overview bands, use them for min/max. */
1317 : /* -------------------------------------------------------------------- */
1318 21 : if (bApproxOK && GetOverviewCount() > 0 && !HasArbitraryOverviews())
1319 : {
1320 : GDALRasterBand *const poBand =
1321 0 : GetRasterSampleOverview(GDALSTAT_APPROX_NUMSAMPLES);
1322 :
1323 0 : if (poBand != nullptr && poBand != this)
1324 : {
1325 0 : auto l_poDS = dynamic_cast<VRTDataset *>(poDS);
1326 0 : if (l_poDS && !l_poDS->m_apoOverviews.empty() &&
1327 0 : dynamic_cast<VRTSourcedRasterBand *>(poBand) != nullptr)
1328 : {
1329 0 : auto apoTmpOverviews = std::move(l_poDS->m_apoOverviews);
1330 0 : l_poDS->m_apoOverviews.clear();
1331 0 : auto eErr = poBand->GDALRasterBand::ComputeRasterMinMax(
1332 : TRUE, adfMinMax);
1333 0 : l_poDS->m_apoOverviews = std::move(apoTmpOverviews);
1334 0 : return eErr;
1335 : }
1336 : else
1337 : {
1338 0 : return poBand->ComputeRasterMinMax(TRUE, adfMinMax);
1339 : }
1340 : }
1341 : }
1342 :
1343 21 : if (IsMosaicOfNonOverlappingSimpleSourcesOfFullRasterNoResAndTypeChange(
1344 : /*bAllowMaxValAdjustment = */ true))
1345 : {
1346 18 : CPLDebugOnly(
1347 : "VRT", "ComputeRasterMinMax(): use optimized code path for mosaic");
1348 :
1349 18 : uint64_t nCoveredArea = 0;
1350 :
1351 : // If source bands have nodata value, we can't use source band's
1352 : // ComputeRasterMinMax() as we don't know if there are pixels actually
1353 : // at the nodata value, so use ComputeStatistics() instead that takes
1354 : // into account that aspect.
1355 18 : bool bUseComputeStatistics = false;
1356 37 : for (int i = 0; i < nSources; ++i)
1357 : {
1358 : auto poSimpleSource =
1359 23 : cpl::down_cast<VRTSimpleSource *>(papoSources[i]);
1360 23 : auto poSimpleSourceBand = poSimpleSource->GetRasterBand();
1361 23 : int bHasNoData = FALSE;
1362 23 : CPL_IGNORE_RET_VAL(poSimpleSourceBand->GetNoDataValue(&bHasNoData));
1363 23 : if (bHasNoData)
1364 : {
1365 4 : bUseComputeStatistics = true;
1366 4 : break;
1367 : }
1368 19 : nCoveredArea +=
1369 19 : static_cast<uint64_t>(poSimpleSourceBand->GetXSize()) *
1370 19 : poSimpleSourceBand->GetYSize();
1371 : }
1372 :
1373 18 : if (bUseComputeStatistics)
1374 : {
1375 : CPLErr eErr;
1376 4 : std::string osLastErrorMsg;
1377 : {
1378 8 : CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);
1379 4 : CPLErrorReset();
1380 : eErr =
1381 8 : ComputeStatistics(bApproxOK, &adfMinMax[0], &adfMinMax[1],
1382 4 : nullptr, nullptr, nullptr, nullptr);
1383 4 : if (eErr == CE_Failure)
1384 : {
1385 1 : osLastErrorMsg = CPLGetLastErrorMsg();
1386 : }
1387 : }
1388 4 : if (eErr == CE_Failure)
1389 : {
1390 1 : if (strstr(osLastErrorMsg.c_str(), "no valid pixels found") !=
1391 : nullptr)
1392 : {
1393 1 : ReportError(CE_Failure, CPLE_AppDefined,
1394 : "Failed to compute min/max, no valid pixels "
1395 : "found in sampling.");
1396 : }
1397 : else
1398 : {
1399 0 : ReportError(CE_Failure, CPLE_AppDefined, "%s",
1400 : osLastErrorMsg.c_str());
1401 : }
1402 : }
1403 4 : return eErr;
1404 : }
1405 :
1406 14 : bool bSignedByte = false;
1407 14 : if (eDataType == GDT_Byte)
1408 : {
1409 8 : EnablePixelTypeSignedByteWarning(false);
1410 : const char *pszPixelType =
1411 8 : GetMetadataItem("PIXELTYPE", "IMAGE_STRUCTURE");
1412 8 : EnablePixelTypeSignedByteWarning(true);
1413 8 : bSignedByte =
1414 8 : pszPixelType != nullptr && EQUAL(pszPixelType, "SIGNEDBYTE");
1415 : }
1416 :
1417 14 : double dfGlobalMin = std::numeric_limits<double>::max();
1418 14 : double dfGlobalMax = -std::numeric_limits<double>::max();
1419 :
1420 : // If the mosaic doesn't cover the whole VRT raster, take into account
1421 : // VRT nodata value.
1422 14 : if (nCoveredArea < static_cast<uint64_t>(nRasterXSize) * nRasterYSize)
1423 : {
1424 2 : if (m_bNoDataValueSet && m_bHideNoDataValue)
1425 : {
1426 1 : if (IsNoDataValueInDataTypeRange())
1427 : {
1428 1 : dfGlobalMin = std::min(dfGlobalMin, m_dfNoDataValue);
1429 1 : dfGlobalMax = std::max(dfGlobalMax, m_dfNoDataValue);
1430 : }
1431 : }
1432 1 : else if (!m_bNoDataValueSet)
1433 : {
1434 0 : dfGlobalMin = std::min(dfGlobalMin, 0.0);
1435 0 : dfGlobalMax = std::max(dfGlobalMax, 0.0);
1436 : }
1437 : }
1438 :
1439 31 : for (int i = 0; i < nSources; ++i)
1440 : {
1441 : auto poSimpleSource =
1442 18 : cpl::down_cast<VRTSimpleSource *>(papoSources[i]);
1443 18 : double adfMinMaxSource[2] = {0};
1444 :
1445 18 : auto poSimpleSourceBand = poSimpleSource->GetRasterBand();
1446 36 : CPLErr eErr = poSimpleSourceBand->ComputeRasterMinMax(
1447 18 : bApproxOK, adfMinMaxSource);
1448 18 : if (eErr == CE_Failure)
1449 : {
1450 1 : return CE_Failure;
1451 : }
1452 : else
1453 : {
1454 18 : if (poSimpleSource->NeedMaxValAdjustment())
1455 : {
1456 2 : const double dfMaxValue =
1457 2 : static_cast<double>(poSimpleSource->m_nMaxValue);
1458 2 : adfMinMaxSource[0] =
1459 2 : std::min(adfMinMaxSource[0], dfMaxValue);
1460 2 : adfMinMaxSource[1] =
1461 2 : std::min(adfMinMaxSource[1], dfMaxValue);
1462 : }
1463 :
1464 18 : if (m_bNoDataValueSet && !m_bHideNoDataValue &&
1465 3 : m_dfNoDataValue >= adfMinMaxSource[0] &&
1466 1 : m_dfNoDataValue <= adfMinMaxSource[1])
1467 : {
1468 1 : return GDALRasterBand::ComputeRasterMinMax(bApproxOK,
1469 1 : adfMinMax);
1470 : }
1471 :
1472 17 : dfGlobalMin = std::min(dfGlobalMin, adfMinMaxSource[0]);
1473 17 : dfGlobalMax = std::max(dfGlobalMax, adfMinMaxSource[1]);
1474 : }
1475 :
1476 : // Early exit if we know we reached theoretical bounds
1477 17 : if (eDataType == GDT_Byte && !bSignedByte && dfGlobalMin == 0.0 &&
1478 0 : dfGlobalMax == 255.0)
1479 : {
1480 0 : break;
1481 : }
1482 : }
1483 :
1484 13 : if (dfGlobalMin > dfGlobalMax)
1485 : {
1486 0 : adfMinMax[0] = 0.0;
1487 0 : adfMinMax[1] = 0.0;
1488 0 : ReportError(CE_Failure, CPLE_AppDefined,
1489 : "Failed to compute min/max, no valid pixels found in "
1490 : "sampling.");
1491 0 : return CE_Failure;
1492 : }
1493 :
1494 13 : adfMinMax[0] = dfGlobalMin;
1495 13 : adfMinMax[1] = dfGlobalMax;
1496 13 : return CE_None;
1497 : }
1498 : else
1499 : {
1500 3 : return GDALRasterBand::ComputeRasterMinMax(bApproxOK, adfMinMax);
1501 : }
1502 : }
1503 :
1504 : // #define naive_update_not_used
1505 :
1506 : namespace
1507 : {
1508 : struct Context
1509 : {
1510 : CPL_DISALLOW_COPY_ASSIGN(Context)
1511 : Context() = default;
1512 :
1513 : // Protected by mutex
1514 : std::mutex oMutex{};
1515 : uint64_t nTotalIteratedPixels = 0;
1516 : uint64_t nLastReportedPixels = 0;
1517 : bool bFailure = false;
1518 : bool bFallbackToBase = false;
1519 : // End of protected by mutex
1520 :
1521 : bool bApproxOK = false;
1522 : GDALProgressFunc pfnProgress = nullptr;
1523 : void *pProgressData = nullptr;
1524 :
1525 : // VRTSourcedRasterBand parameters
1526 : double dfNoDataValue = 0;
1527 : bool bNoDataValueSet = false;
1528 : bool bHideNoDataValue = false;
1529 :
1530 : double dfGlobalMin = std::numeric_limits<double>::max();
1531 : double dfGlobalMax = -std::numeric_limits<double>::max();
1532 : #ifdef naive_update_not_used
1533 : // This native method uses the fact that stddev = sqrt(sum_of_squares/N -
1534 : // mean^2) and that thus sum_of_squares = N * (stddev^2 + mean^2)
1535 : double dfGlobalSum = 0;
1536 : double dfGlobalSumSquare = 0;
1537 : #else
1538 : // This method uses
1539 : // https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Parallel_algorithm
1540 : // which is more numerically robust
1541 : double dfGlobalMean = 0;
1542 : double dfGlobalM2 = 0;
1543 : #endif
1544 : uint64_t nGlobalValidPixels = 0;
1545 : uint64_t nTotalPixelsOfSources = 0;
1546 : };
1547 : } // namespace
1548 :
1549 3 : static void UpdateStatsWithConstantValue(Context &sContext, double dfVal,
1550 : uint64_t nPixelCount)
1551 : {
1552 3 : sContext.dfGlobalMin = std::min(sContext.dfGlobalMin, dfVal);
1553 3 : sContext.dfGlobalMax = std::max(sContext.dfGlobalMax, dfVal);
1554 : #ifdef naive_update_not_used
1555 : sContext.dfGlobalSum += dfVal * nPixelCount;
1556 : sContext.dfGlobalSumSquare += dfVal * dfVal * nPixelCount;
1557 : #else
1558 3 : const auto nNewGlobalValidPixels =
1559 3 : sContext.nGlobalValidPixels + nPixelCount;
1560 3 : const double dfDelta = dfVal - sContext.dfGlobalMean;
1561 3 : sContext.dfGlobalMean += nPixelCount * dfDelta / nNewGlobalValidPixels;
1562 3 : sContext.dfGlobalM2 += dfDelta * dfDelta * nPixelCount *
1563 3 : sContext.nGlobalValidPixels / nNewGlobalValidPixels;
1564 : #endif
1565 3 : sContext.nGlobalValidPixels += nPixelCount;
1566 3 : }
1567 :
1568 : /************************************************************************/
1569 : /* ComputeStatistics() */
1570 : /************************************************************************/
1571 :
1572 30 : CPLErr VRTSourcedRasterBand::ComputeStatistics(int bApproxOK, double *pdfMin,
1573 : double *pdfMax, double *pdfMean,
1574 : double *pdfStdDev,
1575 : GDALProgressFunc pfnProgress,
1576 : void *pProgressData)
1577 :
1578 : {
1579 60 : const std::string osFctId("VRTSourcedRasterBand::ComputeStatistics");
1580 60 : GDALAntiRecursionGuard oGuard(osFctId);
1581 30 : if (oGuard.GetCallDepth() >= 32)
1582 : {
1583 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
1584 0 : return CE_Failure;
1585 : }
1586 :
1587 90 : GDALAntiRecursionGuard oGuard2(oGuard, poDS->GetDescription());
1588 30 : if (oGuard2.GetCallDepth() >= 2)
1589 : {
1590 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
1591 0 : return CE_Failure;
1592 : }
1593 :
1594 : /* -------------------------------------------------------------------- */
1595 : /* If we have overview bands, use them for statistics. */
1596 : /* -------------------------------------------------------------------- */
1597 30 : if (bApproxOK && GetOverviewCount() > 0 && !HasArbitraryOverviews())
1598 : {
1599 : GDALRasterBand *const poBand =
1600 2 : GetRasterSampleOverview(GDALSTAT_APPROX_NUMSAMPLES);
1601 :
1602 2 : if (poBand != nullptr && poBand != this)
1603 : {
1604 2 : auto l_poDS = dynamic_cast<VRTDataset *>(poDS);
1605 : CPLErr eErr;
1606 3 : if (l_poDS && !l_poDS->m_apoOverviews.empty() &&
1607 1 : dynamic_cast<VRTSourcedRasterBand *>(poBand) != nullptr)
1608 : {
1609 2 : auto apoTmpOverviews = std::move(l_poDS->m_apoOverviews);
1610 1 : l_poDS->m_apoOverviews.clear();
1611 1 : eErr = poBand->GDALRasterBand::ComputeStatistics(
1612 : TRUE, pdfMin, pdfMax, pdfMean, pdfStdDev, pfnProgress,
1613 : pProgressData);
1614 1 : l_poDS->m_apoOverviews = std::move(apoTmpOverviews);
1615 : }
1616 : else
1617 : {
1618 1 : eErr = poBand->ComputeStatistics(TRUE, pdfMin, pdfMax, pdfMean,
1619 : pdfStdDev, pfnProgress,
1620 1 : pProgressData);
1621 : }
1622 2 : if (eErr == CE_None && pdfMin && pdfMax && pdfMean && pdfStdDev)
1623 : {
1624 2 : SetMetadataItem("STATISTICS_APPROXIMATE", "YES");
1625 2 : SetMetadataItem(
1626 : "STATISTICS_VALID_PERCENT",
1627 2 : poBand->GetMetadataItem("STATISTICS_VALID_PERCENT"));
1628 2 : SetStatistics(*pdfMin, *pdfMax, *pdfMean, *pdfStdDev);
1629 : }
1630 2 : return eErr;
1631 : }
1632 : }
1633 :
1634 28 : if (IsMosaicOfNonOverlappingSimpleSourcesOfFullRasterNoResAndTypeChange(
1635 : /*bAllowMaxValAdjustment = */ false))
1636 : {
1637 22 : Context sContext;
1638 22 : sContext.bApproxOK = bApproxOK;
1639 22 : sContext.dfNoDataValue = m_dfNoDataValue;
1640 22 : sContext.bNoDataValueSet = m_bNoDataValueSet;
1641 22 : sContext.bHideNoDataValue = m_bHideNoDataValue;
1642 22 : sContext.pfnProgress = pfnProgress;
1643 22 : sContext.pProgressData = pProgressData;
1644 :
1645 : struct Job
1646 : {
1647 : Context *psContext = nullptr;
1648 : GDALRasterBand *poRasterBand = nullptr;
1649 : uint64_t nPixelCount = 0;
1650 : uint64_t nLastAddedPixels = 0;
1651 : uint64_t nValidPixels = 0;
1652 : double dfMin = 0;
1653 : double dfMax = 0;
1654 : double dfMean = 0;
1655 : double dfStdDev = 0;
1656 :
1657 58 : static int CPL_STDCALL ProgressFunc(double dfComplete,
1658 : const char *pszMessage,
1659 : void *pProgressArg)
1660 : {
1661 58 : auto psJob = static_cast<Job *>(pProgressArg);
1662 58 : auto psContext = psJob->psContext;
1663 58 : const uint64_t nNewAddedPixels =
1664 : dfComplete == 1.0
1665 58 : ? psJob->nPixelCount
1666 : : static_cast<uint64_t>(
1667 54 : dfComplete * psJob->nPixelCount + 0.5);
1668 : const auto nUpdateThreshold =
1669 116 : std::min(psContext->nTotalPixelsOfSources / 1000,
1670 58 : static_cast<uint64_t>(1000 * 1000));
1671 116 : std::lock_guard<std::mutex> oLock(psContext->oMutex);
1672 58 : psContext->nTotalIteratedPixels +=
1673 58 : (nNewAddedPixels - psJob->nLastAddedPixels);
1674 58 : psJob->nLastAddedPixels = nNewAddedPixels;
1675 58 : if (psContext->nTotalIteratedPixels ==
1676 58 : psContext->nTotalPixelsOfSources)
1677 : {
1678 2 : psContext->nLastReportedPixels =
1679 2 : psContext->nTotalIteratedPixels;
1680 2 : return psContext->pfnProgress(1.0, pszMessage,
1681 2 : psContext->pProgressData);
1682 : }
1683 56 : else if (psContext->nTotalIteratedPixels -
1684 56 : psContext->nLastReportedPixels >
1685 : nUpdateThreshold)
1686 : {
1687 48 : psContext->nLastReportedPixels =
1688 48 : psContext->nTotalIteratedPixels;
1689 96 : return psContext->pfnProgress(
1690 48 : static_cast<double>(psContext->nTotalIteratedPixels) /
1691 48 : psContext->nTotalPixelsOfSources,
1692 48 : pszMessage, psContext->pProgressData);
1693 : }
1694 8 : return 1;
1695 : }
1696 :
1697 32 : static void UpdateStats(const Job *psJob)
1698 : {
1699 32 : const auto nValidPixels = psJob->nValidPixels;
1700 32 : auto psContext = psJob->psContext;
1701 32 : if (nValidPixels > 0)
1702 : {
1703 24 : psContext->dfGlobalMin =
1704 24 : std::min(psContext->dfGlobalMin, psJob->dfMin);
1705 24 : psContext->dfGlobalMax =
1706 24 : std::max(psContext->dfGlobalMax, psJob->dfMax);
1707 : #ifdef naive_update_not_used
1708 : psContext->dfGlobalSum += nValidPixels * psJob->dfMean;
1709 : psContext->dfGlobalSumSquare +=
1710 : nValidPixels * (psJob->dfStdDev * psJob->dfStdDev +
1711 : psJob->dfMean * psJob->dfMean);
1712 : psContext->nGlobalValidPixels += nValidPixels;
1713 : #else
1714 24 : const auto nNewGlobalValidPixels =
1715 24 : psContext->nGlobalValidPixels + nValidPixels;
1716 24 : const double dfDelta =
1717 24 : psJob->dfMean - psContext->dfGlobalMean;
1718 24 : psContext->dfGlobalMean +=
1719 24 : nValidPixels * dfDelta / nNewGlobalValidPixels;
1720 24 : psContext->dfGlobalM2 +=
1721 24 : nValidPixels * psJob->dfStdDev * psJob->dfStdDev +
1722 24 : dfDelta * dfDelta * nValidPixels *
1723 24 : psContext->nGlobalValidPixels /
1724 : nNewGlobalValidPixels;
1725 24 : psContext->nGlobalValidPixels = nNewGlobalValidPixels;
1726 : #endif
1727 : }
1728 32 : int bHasNoData = FALSE;
1729 : const double dfNoDataValue =
1730 32 : psJob->poRasterBand->GetNoDataValue(&bHasNoData);
1731 9 : if (nValidPixels < psJob->nPixelCount && bHasNoData &&
1732 50 : !std::isnan(dfNoDataValue) &&
1733 9 : (!psContext->bNoDataValueSet ||
1734 7 : dfNoDataValue != psContext->dfNoDataValue))
1735 : {
1736 : const auto eBandDT =
1737 2 : psJob->poRasterBand->GetRasterDataType();
1738 : // Check that the band nodata value is in the range of the
1739 : // original raster type
1740 : GByte abyTempBuffer[2 * sizeof(double)];
1741 2 : CPLAssert(GDALGetDataTypeSizeBytes(eBandDT) <=
1742 : static_cast<int>(sizeof(abyTempBuffer)));
1743 2 : GDALCopyWords(&dfNoDataValue, GDT_Float64, 0,
1744 : &abyTempBuffer[0], eBandDT, 0, 1);
1745 2 : double dfNoDataValueAfter = dfNoDataValue;
1746 2 : GDALCopyWords(&abyTempBuffer[0], eBandDT, 0,
1747 : &dfNoDataValueAfter, GDT_Float64, 0, 1);
1748 4 : if (!std::isfinite(dfNoDataValue) ||
1749 2 : std::fabs(dfNoDataValueAfter - dfNoDataValue) < 1.0)
1750 : {
1751 2 : UpdateStatsWithConstantValue(
1752 : *psContext, dfNoDataValueAfter,
1753 2 : psJob->nPixelCount - nValidPixels);
1754 : }
1755 : }
1756 32 : }
1757 : };
1758 :
1759 35 : const auto JobRunner = [](void *pData)
1760 : {
1761 35 : auto psJob = static_cast<Job *>(pData);
1762 35 : auto psContext = psJob->psContext;
1763 : {
1764 35 : std::lock_guard<std::mutex> oLock(psContext->oMutex);
1765 35 : if (psContext->bFallbackToBase || psContext->bFailure)
1766 1 : return;
1767 : }
1768 :
1769 34 : auto poSimpleSourceBand = psJob->poRasterBand;
1770 34 : psJob->nPixelCount =
1771 34 : static_cast<uint64_t>(poSimpleSourceBand->GetXSize()) *
1772 34 : poSimpleSourceBand->GetYSize();
1773 :
1774 34 : CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);
1775 34 : CPLErr eErr = poSimpleSourceBand->ComputeStatistics(
1776 34 : psContext->bApproxOK, &psJob->dfMin, &psJob->dfMax,
1777 : &psJob->dfMean, &psJob->dfStdDev,
1778 40 : psContext->pfnProgress == nullptr ||
1779 6 : psContext->pfnProgress == GDALDummyProgress
1780 : ? GDALDummyProgress
1781 : : Job::ProgressFunc,
1782 34 : psJob);
1783 : const char *pszValidPercent =
1784 34 : poSimpleSourceBand->GetMetadataItem("STATISTICS_VALID_PERCENT");
1785 34 : psJob->nValidPixels =
1786 : pszValidPercent
1787 34 : ? static_cast<uint64_t>(CPLAtof(pszValidPercent) *
1788 34 : psJob->nPixelCount / 100.0)
1789 : : psJob->nPixelCount;
1790 34 : if (eErr == CE_Failure)
1791 : {
1792 16 : if (pszValidPercent != nullptr &&
1793 8 : CPLAtof(pszValidPercent) == 0.0)
1794 : {
1795 : // ok: no valid sample
1796 : }
1797 : else
1798 : {
1799 0 : std::lock_guard<std::mutex> oLock(psContext->oMutex);
1800 0 : psContext->bFailure = true;
1801 : }
1802 : }
1803 : else
1804 : {
1805 26 : int bHasNoData = FALSE;
1806 26 : CPL_IGNORE_RET_VAL(
1807 26 : psJob->poRasterBand->GetNoDataValue(&bHasNoData));
1808 26 : if (!bHasNoData && psContext->bNoDataValueSet &&
1809 8 : !psContext->bHideNoDataValue &&
1810 6 : psContext->dfNoDataValue >= psJob->dfMin &&
1811 2 : psContext->dfNoDataValue <= psJob->dfMax)
1812 : {
1813 2 : std::lock_guard<std::mutex> oLock(psContext->oMutex);
1814 2 : psJob->psContext->bFallbackToBase = true;
1815 2 : return;
1816 : }
1817 : }
1818 : };
1819 :
1820 22 : CPLWorkerThreadPool *poThreadPool = nullptr;
1821 : int nThreads =
1822 22 : nSources > 1
1823 22 : ? VRTDataset::GetNumThreads(dynamic_cast<VRTDataset *>(poDS))
1824 22 : : 0;
1825 22 : if (nThreads > 1024)
1826 0 : nThreads = 1024; // to please Coverity
1827 22 : if (nThreads > 1)
1828 : {
1829 : // Check that all sources refer to different datasets
1830 : // before allowing multithreaded access
1831 : // If the datasets belong to the MEM driver, check GDALDataset*
1832 : // pointer values. Otherwise use dataset name.
1833 24 : std::set<std::string> oSetDatasetNames;
1834 24 : std::set<GDALDataset *> oSetDatasetPointers;
1835 36 : for (int i = 0; i < nSources; ++i)
1836 : {
1837 : auto poSimpleSource =
1838 24 : cpl::down_cast<VRTSimpleSource *>(papoSources[i]);
1839 24 : assert(poSimpleSource);
1840 24 : auto poSimpleSourceBand = poSimpleSource->GetRasterBand();
1841 24 : assert(poSimpleSourceBand);
1842 24 : auto poSourceDataset = poSimpleSourceBand->GetDataset();
1843 24 : if (poSourceDataset == nullptr)
1844 : {
1845 0 : nThreads = 0;
1846 0 : break;
1847 : }
1848 24 : auto poDriver = poSourceDataset->GetDriver();
1849 24 : if (poDriver && EQUAL(poDriver->GetDescription(), "MEM"))
1850 : {
1851 24 : if (oSetDatasetPointers.find(poSourceDataset) !=
1852 48 : oSetDatasetPointers.end())
1853 : {
1854 0 : nThreads = 0;
1855 0 : break;
1856 : }
1857 24 : oSetDatasetPointers.insert(poSourceDataset);
1858 : }
1859 : else
1860 : {
1861 0 : if (oSetDatasetNames.find(
1862 0 : poSourceDataset->GetDescription()) !=
1863 0 : oSetDatasetNames.end())
1864 : {
1865 0 : nThreads = 0;
1866 0 : break;
1867 : }
1868 0 : oSetDatasetNames.insert(poSourceDataset->GetDescription());
1869 : }
1870 : }
1871 12 : if (nThreads > 1)
1872 : {
1873 12 : poThreadPool = GDALGetGlobalThreadPool(nThreads);
1874 : }
1875 : }
1876 :
1877 : // Compute total number of pixels of sources
1878 57 : for (int i = 0; i < nSources; ++i)
1879 : {
1880 35 : auto poSimpleSource =
1881 35 : static_cast<VRTSimpleSource *>(papoSources[i]);
1882 35 : assert(poSimpleSource);
1883 35 : auto poSimpleSourceBand = poSimpleSource->GetRasterBand();
1884 35 : assert(poSimpleSourceBand);
1885 35 : sContext.nTotalPixelsOfSources +=
1886 35 : static_cast<uint64_t>(poSimpleSourceBand->GetXSize()) *
1887 35 : poSimpleSourceBand->GetYSize();
1888 : }
1889 :
1890 22 : if (poThreadPool)
1891 : {
1892 12 : CPLDebugOnly("VRT", "ComputeStatistics(): use optimized "
1893 : "multi-threaded code path for mosaic");
1894 24 : std::vector<Job> asJobs(nSources);
1895 24 : auto poQueue = poThreadPool->CreateJobQueue();
1896 36 : for (int i = 0; i < nSources; ++i)
1897 : {
1898 24 : auto poSimpleSource =
1899 24 : static_cast<VRTSimpleSource *>(papoSources[i]);
1900 24 : assert(poSimpleSource);
1901 24 : auto poSimpleSourceBand = poSimpleSource->GetRasterBand();
1902 24 : assert(poSimpleSourceBand);
1903 24 : asJobs[i].psContext = &sContext;
1904 24 : asJobs[i].poRasterBand = poSimpleSourceBand;
1905 24 : if (!poQueue->SubmitJob(JobRunner, &asJobs[i]))
1906 : {
1907 0 : sContext.bFailure = true;
1908 0 : break;
1909 : }
1910 : }
1911 12 : poQueue->WaitCompletion();
1912 12 : if (!(sContext.bFailure || sContext.bFallbackToBase))
1913 : {
1914 33 : for (int i = 0; i < nSources; ++i)
1915 : {
1916 22 : Job::UpdateStats(&asJobs[i]);
1917 : }
1918 : }
1919 : }
1920 : else
1921 : {
1922 10 : CPLDebugOnly(
1923 : "VRT",
1924 : "ComputeStatistics(): use optimized code path for mosaic");
1925 20 : for (int i = 0; i < nSources; ++i)
1926 : {
1927 11 : auto poSimpleSource =
1928 11 : static_cast<VRTSimpleSource *>(papoSources[i]);
1929 11 : assert(poSimpleSource);
1930 11 : auto poSimpleSourceBand = poSimpleSource->GetRasterBand();
1931 11 : assert(poSimpleSourceBand);
1932 11 : Job sJob;
1933 11 : sJob.psContext = &sContext;
1934 11 : sJob.poRasterBand = poSimpleSourceBand;
1935 11 : JobRunner(&sJob);
1936 11 : if (sContext.bFailure || sContext.bFallbackToBase)
1937 : break;
1938 10 : Job::UpdateStats(&sJob);
1939 : }
1940 : }
1941 :
1942 22 : if (sContext.bFailure)
1943 0 : return CE_Failure;
1944 22 : if (sContext.bFallbackToBase)
1945 : {
1946 : // If the VRT band nodata value is in the [min, max] range
1947 : // of the source and that the source has no nodata value set,
1948 : // then we can't use the optimization.
1949 2 : CPLDebugOnly("VRT", "ComputeStatistics(): revert back to "
1950 : "generic case because of nodata value in range "
1951 : "of source raster");
1952 2 : return GDALRasterBand::ComputeStatistics(
1953 : bApproxOK, pdfMin, pdfMax, pdfMean, pdfStdDev, pfnProgress,
1954 2 : pProgressData);
1955 : }
1956 :
1957 20 : const uint64_t nTotalPixels =
1958 20 : static_cast<uint64_t>(nRasterXSize) * nRasterYSize;
1959 7 : if (m_bNoDataValueSet && m_bHideNoDataValue &&
1960 27 : !std::isnan(m_dfNoDataValue) && IsNoDataValueInDataTypeRange())
1961 : {
1962 1 : UpdateStatsWithConstantValue(sContext, m_dfNoDataValue,
1963 : nTotalPixels -
1964 1 : sContext.nGlobalValidPixels);
1965 : }
1966 19 : else if (!m_bNoDataValueSet)
1967 : {
1968 13 : sContext.nGlobalValidPixels = nTotalPixels;
1969 : }
1970 :
1971 : #ifdef naive_update_not_used
1972 : const double dfGlobalMean =
1973 : sContext.nGlobalValidPixels > 0
1974 : ? sContext.dfGlobalSum / sContext.nGlobalValidPixels
1975 : : 0;
1976 : const double dfGlobalStdDev =
1977 : sContext.nGlobalValidPixels > 0
1978 : ? sqrt(sContext.dfGlobalSumSquare /
1979 : sContext.nGlobalValidPixels -
1980 : dfGlobalMean * dfGlobalMean)
1981 : : 0;
1982 : #else
1983 20 : const double dfGlobalMean = sContext.dfGlobalMean;
1984 : const double dfGlobalStdDev =
1985 20 : sContext.nGlobalValidPixels > 0
1986 20 : ? sqrt(sContext.dfGlobalM2 / sContext.nGlobalValidPixels)
1987 20 : : 0;
1988 : #endif
1989 20 : if (sContext.nGlobalValidPixels > 0)
1990 : {
1991 18 : if (bApproxOK)
1992 : {
1993 1 : SetMetadataItem("STATISTICS_APPROXIMATE", "YES");
1994 : }
1995 17 : else if (GetMetadataItem("STATISTICS_APPROXIMATE"))
1996 : {
1997 0 : SetMetadataItem("STATISTICS_APPROXIMATE", nullptr);
1998 : }
1999 18 : SetStatistics(sContext.dfGlobalMin, sContext.dfGlobalMax,
2000 18 : dfGlobalMean, dfGlobalStdDev);
2001 : }
2002 : else
2003 : {
2004 2 : sContext.dfGlobalMin = 0.0;
2005 2 : sContext.dfGlobalMax = 0.0;
2006 : }
2007 :
2008 20 : SetValidPercent(nTotalPixels, sContext.nGlobalValidPixels);
2009 :
2010 20 : if (pdfMin)
2011 16 : *pdfMin = sContext.dfGlobalMin;
2012 20 : if (pdfMax)
2013 16 : *pdfMax = sContext.dfGlobalMax;
2014 20 : if (pdfMean)
2015 12 : *pdfMean = dfGlobalMean;
2016 20 : if (pdfStdDev)
2017 12 : *pdfStdDev = dfGlobalStdDev;
2018 :
2019 20 : if (sContext.nGlobalValidPixels == 0)
2020 : {
2021 2 : ReportError(CE_Failure, CPLE_AppDefined,
2022 : "Failed to compute statistics, no valid pixels found "
2023 : "in sampling.");
2024 : }
2025 :
2026 20 : return sContext.nGlobalValidPixels > 0 ? CE_None : CE_Failure;
2027 : }
2028 : else
2029 : {
2030 6 : return GDALRasterBand::ComputeStatistics(bApproxOK, pdfMin, pdfMax,
2031 : pdfMean, pdfStdDev,
2032 6 : pfnProgress, pProgressData);
2033 : }
2034 : }
2035 :
2036 : /************************************************************************/
2037 : /* GetHistogram() */
2038 : /************************************************************************/
2039 :
2040 7 : CPLErr VRTSourcedRasterBand::GetHistogram(double dfMin, double dfMax,
2041 : int nBuckets, GUIntBig *panHistogram,
2042 : int bIncludeOutOfRange, int bApproxOK,
2043 : GDALProgressFunc pfnProgress,
2044 : void *pProgressData)
2045 :
2046 : {
2047 : /* -------------------------------------------------------------------- */
2048 : /* If we have overviews, use them for the histogram. */
2049 : /* -------------------------------------------------------------------- */
2050 7 : if (bApproxOK && GetOverviewCount() > 0 && !HasArbitraryOverviews())
2051 : {
2052 : // FIXME: Should we use the most reduced overview here or use some
2053 : // minimum number of samples like GDALRasterBand::ComputeStatistics()
2054 : // does?
2055 1 : GDALRasterBand *poBand = GetRasterSampleOverview(0);
2056 :
2057 1 : if (poBand != nullptr && poBand != this)
2058 : {
2059 1 : auto l_poDS = dynamic_cast<VRTDataset *>(poDS);
2060 2 : if (l_poDS && !l_poDS->m_apoOverviews.empty() &&
2061 1 : dynamic_cast<VRTSourcedRasterBand *>(poBand) != nullptr)
2062 : {
2063 1 : auto apoTmpOverviews = std::move(l_poDS->m_apoOverviews);
2064 1 : l_poDS->m_apoOverviews.clear();
2065 1 : auto eErr = poBand->GDALRasterBand::GetHistogram(
2066 : dfMin, dfMax, nBuckets, panHistogram, bIncludeOutOfRange,
2067 : bApproxOK, pfnProgress, pProgressData);
2068 1 : l_poDS->m_apoOverviews = std::move(apoTmpOverviews);
2069 1 : return eErr;
2070 : }
2071 : else
2072 : {
2073 0 : return poBand->GetHistogram(
2074 : dfMin, dfMax, nBuckets, panHistogram, bIncludeOutOfRange,
2075 0 : bApproxOK, pfnProgress, pProgressData);
2076 : }
2077 : }
2078 : }
2079 :
2080 6 : if (nSources != 1)
2081 2 : return VRTRasterBand::GetHistogram(dfMin, dfMax, nBuckets, panHistogram,
2082 : bIncludeOutOfRange, bApproxOK,
2083 2 : pfnProgress, pProgressData);
2084 :
2085 4 : if (pfnProgress == nullptr)
2086 4 : pfnProgress = GDALDummyProgress;
2087 :
2088 8 : const std::string osFctId("VRTSourcedRasterBand::GetHistogram");
2089 8 : GDALAntiRecursionGuard oGuard(osFctId);
2090 4 : if (oGuard.GetCallDepth() >= 32)
2091 : {
2092 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
2093 0 : return CE_Failure;
2094 : }
2095 :
2096 12 : GDALAntiRecursionGuard oGuard2(oGuard, poDS->GetDescription());
2097 4 : if (oGuard2.GetCallDepth() >= 2)
2098 : {
2099 0 : CPLError(CE_Failure, CPLE_AppDefined, "Recursion detected");
2100 0 : return CE_Failure;
2101 : }
2102 :
2103 : /* -------------------------------------------------------------------- */
2104 : /* Try with source bands. */
2105 : /* -------------------------------------------------------------------- */
2106 4 : const CPLErr eErr = papoSources[0]->GetHistogram(
2107 : GetXSize(), GetYSize(), dfMin, dfMax, nBuckets, panHistogram,
2108 4 : bIncludeOutOfRange, bApproxOK, pfnProgress, pProgressData);
2109 4 : if (eErr != CE_None)
2110 : {
2111 0 : const CPLErr eErr2 = GDALRasterBand::GetHistogram(
2112 : dfMin, dfMax, nBuckets, panHistogram, bIncludeOutOfRange, bApproxOK,
2113 : pfnProgress, pProgressData);
2114 0 : return eErr2;
2115 : }
2116 :
2117 4 : SetDefaultHistogram(dfMin, dfMax, nBuckets, panHistogram);
2118 :
2119 4 : return CE_None;
2120 : }
2121 :
2122 : /************************************************************************/
2123 : /* AddSource() */
2124 : /************************************************************************/
2125 :
2126 12948 : CPLErr VRTSourcedRasterBand::AddSource(VRTSource *poNewSource)
2127 :
2128 : {
2129 12948 : nSources++;
2130 :
2131 12948 : papoSources = static_cast<VRTSource **>(
2132 12948 : CPLRealloc(papoSources, sizeof(void *) * nSources));
2133 12948 : papoSources[nSources - 1] = poNewSource;
2134 :
2135 12948 : auto l_poDS = static_cast<VRTDataset *>(poDS);
2136 12948 : l_poDS->SetNeedsFlush();
2137 12948 : l_poDS->SourceAdded();
2138 :
2139 12948 : if (poNewSource->IsSimpleSource())
2140 : {
2141 12922 : VRTSimpleSource *poSS = static_cast<VRTSimpleSource *>(poNewSource);
2142 12922 : if (GetMetadataItem("NBITS", "IMAGE_STRUCTURE") != nullptr)
2143 : {
2144 47 : int nBits = atoi(GetMetadataItem("NBITS", "IMAGE_STRUCTURE"));
2145 47 : if (nBits >= 1 && nBits <= 31)
2146 : {
2147 47 : poSS->SetMaxValue(static_cast<int>((1U << nBits) - 1));
2148 : }
2149 : }
2150 : }
2151 :
2152 12948 : return CE_None;
2153 : }
2154 :
2155 : /*! @endcond */
2156 :
2157 : /************************************************************************/
2158 : /* VRTAddSource() */
2159 : /************************************************************************/
2160 :
2161 : /**
2162 : * @see VRTSourcedRasterBand::AddSource().
2163 : */
2164 :
2165 0 : CPLErr CPL_STDCALL VRTAddSource(VRTSourcedRasterBandH hVRTBand,
2166 : VRTSourceH hNewSource)
2167 : {
2168 0 : VALIDATE_POINTER1(hVRTBand, "VRTAddSource", CE_Failure);
2169 :
2170 0 : return reinterpret_cast<VRTSourcedRasterBand *>(hVRTBand)->AddSource(
2171 0 : reinterpret_cast<VRTSource *>(hNewSource));
2172 : }
2173 :
2174 : /*! @cond Doxygen_Suppress */
2175 :
2176 : /************************************************************************/
2177 : /* XMLInit() */
2178 : /************************************************************************/
2179 :
2180 1206 : CPLErr VRTSourcedRasterBand::XMLInit(const CPLXMLNode *psTree,
2181 : const char *pszVRTPath,
2182 : VRTMapSharedResources &oMapSharedSources)
2183 :
2184 : {
2185 : {
2186 : const CPLErr eErr =
2187 1206 : VRTRasterBand::XMLInit(psTree, pszVRTPath, oMapSharedSources);
2188 1206 : if (eErr != CE_None)
2189 0 : return eErr;
2190 : }
2191 :
2192 : /* -------------------------------------------------------------------- */
2193 : /* Process sources. */
2194 : /* -------------------------------------------------------------------- */
2195 : VRTDriver *const poDriver =
2196 1206 : static_cast<VRTDriver *>(GDALGetDriverByName("VRT"));
2197 :
2198 1206 : for (const CPLXMLNode *psChild = psTree->psChild;
2199 8751 : psChild != nullptr && poDriver != nullptr; psChild = psChild->psNext)
2200 : {
2201 7566 : if (psChild->eType != CXT_Element)
2202 2858 : continue;
2203 :
2204 4708 : CPLErrorReset();
2205 : VRTSource *const poSource =
2206 4708 : poDriver->ParseSource(psChild, pszVRTPath, oMapSharedSources);
2207 4708 : if (poSource != nullptr)
2208 3098 : AddSource(poSource);
2209 1610 : else if (CPLGetLastErrorType() != CE_None)
2210 21 : return CE_Failure;
2211 : }
2212 :
2213 : /* -------------------------------------------------------------------- */
2214 : /* Done. */
2215 : /* -------------------------------------------------------------------- */
2216 : const char *pszSubclass =
2217 1185 : CPLGetXMLValue(psTree, "subclass", "VRTSourcedRasterBand");
2218 1185 : if (nSources == 0 && !EQUAL(pszSubclass, "VRTDerivedRasterBand"))
2219 147 : CPLDebug("VRT", "No valid sources found for band in VRT file %s",
2220 147 : GetDataset() ? GetDataset()->GetDescription() : "");
2221 :
2222 1185 : return CE_None;
2223 : }
2224 :
2225 : /************************************************************************/
2226 : /* SerializeToXML() */
2227 : /************************************************************************/
2228 :
2229 605 : CPLXMLNode *VRTSourcedRasterBand::SerializeToXML(const char *pszVRTPath,
2230 : bool &bHasWarnedAboutRAMUsage,
2231 : size_t &nAccRAMUsage)
2232 :
2233 : {
2234 605 : CPLXMLNode *psTree = VRTRasterBand::SerializeToXML(
2235 : pszVRTPath, bHasWarnedAboutRAMUsage, nAccRAMUsage);
2236 605 : CPLXMLNode *psLastChild = psTree->psChild;
2237 2320 : while (psLastChild != nullptr && psLastChild->psNext != nullptr)
2238 1715 : psLastChild = psLastChild->psNext;
2239 :
2240 : /* -------------------------------------------------------------------- */
2241 : /* Process Sources. */
2242 : /* -------------------------------------------------------------------- */
2243 :
2244 605 : GIntBig nUsableRAM = -1;
2245 :
2246 3143 : for (int iSource = 0; iSource < nSources; iSource++)
2247 : {
2248 : CPLXMLNode *const psXMLSrc =
2249 2538 : papoSources[iSource]->SerializeToXML(pszVRTPath);
2250 :
2251 2538 : if (psXMLSrc == nullptr)
2252 0 : break;
2253 :
2254 : // Creating the CPLXMLNode tree representation of a VRT can easily
2255 : // take several times RAM usage than its string serialization, or its
2256 : // internal representation in the driver.
2257 : // We multiply the estimate by a factor of 2, experimentally found to
2258 : // be more realistic than the conservative raw estimate.
2259 2538 : nAccRAMUsage += 2 * CPLXMLNodeGetRAMUsageEstimate(psXMLSrc);
2260 2538 : if (!bHasWarnedAboutRAMUsage && nAccRAMUsage > 512 * 1024 * 1024)
2261 : {
2262 0 : if (nUsableRAM < 0)
2263 0 : nUsableRAM = CPLGetUsablePhysicalRAM();
2264 0 : if (nUsableRAM > 0 &&
2265 0 : nAccRAMUsage > static_cast<uint64_t>(nUsableRAM) / 10 * 8)
2266 : {
2267 0 : bHasWarnedAboutRAMUsage = true;
2268 0 : CPLError(CE_Warning, CPLE_AppDefined,
2269 : "Serialization of this VRT file has already consumed "
2270 : "at least %.02f GB of RAM over a total of %.02f. This "
2271 : "process may abort",
2272 0 : double(nAccRAMUsage) / (1024 * 1024 * 1024),
2273 0 : double(nUsableRAM) / (1024 * 1024 * 1024));
2274 : }
2275 : }
2276 :
2277 2538 : if (psLastChild == nullptr)
2278 0 : psTree->psChild = psXMLSrc;
2279 : else
2280 2538 : psLastChild->psNext = psXMLSrc;
2281 2538 : psLastChild = psXMLSrc;
2282 : }
2283 :
2284 605 : return psTree;
2285 : }
2286 :
2287 : /************************************************************************/
2288 : /* SkipBufferInitialization() */
2289 : /************************************************************************/
2290 :
2291 48879 : bool VRTSourcedRasterBand::SkipBufferInitialization()
2292 : {
2293 48879 : if (m_nSkipBufferInitialization >= 0)
2294 44317 : return m_nSkipBufferInitialization != 0;
2295 : /* -------------------------------------------------------------------- */
2296 : /* Check if we can avoid buffer initialization. */
2297 : /* -------------------------------------------------------------------- */
2298 :
2299 : // Note: if one day we do alpha compositing, we will need to check that.
2300 4562 : m_nSkipBufferInitialization = FALSE;
2301 4562 : if (nSources != 1 || !papoSources[0]->IsSimpleSource())
2302 : {
2303 659 : return false;
2304 : }
2305 3903 : VRTSimpleSource *poSS = static_cast<VRTSimpleSource *>(papoSources[0]);
2306 3903 : if (poSS->GetType() == VRTSimpleSource::GetTypeStatic())
2307 : {
2308 3628 : auto l_poBand = poSS->GetRasterBand();
2309 3616 : if (l_poBand != nullptr && poSS->m_dfSrcXOff >= 0.0 &&
2310 3578 : poSS->m_dfSrcYOff >= 0.0 &&
2311 3578 : poSS->m_dfSrcXOff + poSS->m_dfSrcXSize <= l_poBand->GetXSize() &&
2312 3572 : poSS->m_dfSrcYOff + poSS->m_dfSrcYSize <= l_poBand->GetYSize() &&
2313 3572 : poSS->m_dfDstXOff <= 0.0 && poSS->m_dfDstYOff <= 0.0 &&
2314 10799 : poSS->m_dfDstXOff + poSS->m_dfDstXSize >= nRasterXSize &&
2315 3555 : poSS->m_dfDstYOff + poSS->m_dfDstYSize >= nRasterYSize)
2316 : {
2317 3552 : m_nSkipBufferInitialization = TRUE;
2318 : }
2319 : }
2320 3903 : return m_nSkipBufferInitialization != 0;
2321 : }
2322 :
2323 : /************************************************************************/
2324 : /* ConfigureSource() */
2325 : /************************************************************************/
2326 :
2327 9310 : void VRTSourcedRasterBand::ConfigureSource(VRTSimpleSource *poSimpleSource,
2328 : GDALRasterBand *poSrcBand,
2329 : int bAddAsMaskBand, double dfSrcXOff,
2330 : double dfSrcYOff, double dfSrcXSize,
2331 : double dfSrcYSize, double dfDstXOff,
2332 : double dfDstYOff, double dfDstXSize,
2333 : double dfDstYSize)
2334 : {
2335 : /* -------------------------------------------------------------------- */
2336 : /* Default source and dest rectangles. */
2337 : /* -------------------------------------------------------------------- */
2338 9310 : if (dfSrcYSize == -1)
2339 : {
2340 294 : dfSrcXOff = 0;
2341 294 : dfSrcYOff = 0;
2342 294 : dfSrcXSize = poSrcBand->GetXSize();
2343 294 : dfSrcYSize = poSrcBand->GetYSize();
2344 : }
2345 :
2346 9310 : if (dfDstYSize == -1)
2347 : {
2348 294 : dfDstXOff = 0;
2349 294 : dfDstYOff = 0;
2350 294 : dfDstXSize = nRasterXSize;
2351 294 : dfDstYSize = nRasterYSize;
2352 : }
2353 :
2354 9310 : if (bAddAsMaskBand)
2355 32 : poSimpleSource->SetSrcMaskBand(poSrcBand);
2356 : else
2357 9278 : poSimpleSource->SetSrcBand(poSrcBand);
2358 :
2359 9310 : poSimpleSource->SetSrcWindow(dfSrcXOff, dfSrcYOff, dfSrcXSize, dfSrcYSize);
2360 9310 : poSimpleSource->SetDstWindow(dfDstXOff, dfDstYOff, dfDstXSize, dfDstYSize);
2361 :
2362 : /* -------------------------------------------------------------------- */
2363 : /* If we can get the associated GDALDataset, add a reference to it.*/
2364 : /* -------------------------------------------------------------------- */
2365 9310 : GDALDataset *poSrcBandDataset = poSrcBand->GetDataset();
2366 9310 : if (poSrcBandDataset != nullptr)
2367 : {
2368 : VRTDataset *poVRTSrcBandDataset =
2369 9310 : dynamic_cast<VRTDataset *>(poSrcBandDataset);
2370 9310 : if (poVRTSrcBandDataset && !poVRTSrcBandDataset->m_bCanTakeRef)
2371 : {
2372 : // Situation triggered by VRTDataset::AddVirtualOverview()
2373 : // We create an overview dataset that is a VRT of a reduction of
2374 : // ourselves. But we don't want to take a reference on ourselves,
2375 : // otherwise this will prevent us to be closed in number of
2376 : // circumstances
2377 60 : poSimpleSource->m_bDropRefOnSrcBand = false;
2378 : }
2379 : else
2380 : {
2381 9250 : poSrcBandDataset->Reference();
2382 : }
2383 : }
2384 9310 : }
2385 :
2386 : /************************************************************************/
2387 : /* AddSimpleSource() */
2388 : /************************************************************************/
2389 :
2390 360 : CPLErr VRTSourcedRasterBand::AddSimpleSource(
2391 : const char *pszFilename, int nBandIn, double dfSrcXOff, double dfSrcYOff,
2392 : double dfSrcXSize, double dfSrcYSize, double dfDstXOff, double dfDstYOff,
2393 : double dfDstXSize, double dfDstYSize, const char *pszResampling,
2394 : double dfNoDataValueIn)
2395 :
2396 : {
2397 : /* -------------------------------------------------------------------- */
2398 : /* Create source. */
2399 : /* -------------------------------------------------------------------- */
2400 360 : VRTSimpleSource *poSimpleSource = nullptr;
2401 :
2402 360 : if (pszResampling != nullptr && STARTS_WITH_CI(pszResampling, "aver"))
2403 : {
2404 0 : auto poAveragedSource = new VRTAveragedSource();
2405 0 : poSimpleSource = poAveragedSource;
2406 0 : if (dfNoDataValueIn != VRT_NODATA_UNSET)
2407 0 : poAveragedSource->SetNoDataValue(dfNoDataValueIn);
2408 : }
2409 : else
2410 : {
2411 360 : poSimpleSource = new VRTSimpleSource();
2412 360 : if (dfNoDataValueIn != VRT_NODATA_UNSET)
2413 0 : CPLError(
2414 : CE_Warning, CPLE_AppDefined,
2415 : "NODATA setting not currently supported for nearest "
2416 : "neighbour sampled simple sources on Virtual Datasources.");
2417 : }
2418 :
2419 360 : poSimpleSource->SetSrcBand(pszFilename, nBandIn);
2420 :
2421 360 : poSimpleSource->SetSrcWindow(dfSrcXOff, dfSrcYOff, dfSrcXSize, dfSrcYSize);
2422 360 : poSimpleSource->SetDstWindow(dfDstXOff, dfDstYOff, dfDstXSize, dfDstYSize);
2423 :
2424 : /* -------------------------------------------------------------------- */
2425 : /* add to list. */
2426 : /* -------------------------------------------------------------------- */
2427 360 : return AddSource(poSimpleSource);
2428 : }
2429 :
2430 : /************************************************************************/
2431 : /* AddSimpleSource() */
2432 : /************************************************************************/
2433 :
2434 1554 : CPLErr VRTSourcedRasterBand::AddSimpleSource(
2435 : GDALRasterBand *poSrcBand, double dfSrcXOff, double dfSrcYOff,
2436 : double dfSrcXSize, double dfSrcYSize, double dfDstXOff, double dfDstYOff,
2437 : double dfDstXSize, double dfDstYSize, const char *pszResampling,
2438 : double dfNoDataValueIn)
2439 :
2440 : {
2441 : /* -------------------------------------------------------------------- */
2442 : /* Create source. */
2443 : /* -------------------------------------------------------------------- */
2444 1554 : VRTSimpleSource *poSimpleSource = nullptr;
2445 :
2446 1554 : if (pszResampling != nullptr && STARTS_WITH_CI(pszResampling, "aver"))
2447 : {
2448 0 : auto poAveragedSource = new VRTAveragedSource();
2449 0 : poSimpleSource = poAveragedSource;
2450 0 : if (dfNoDataValueIn != VRT_NODATA_UNSET)
2451 0 : poAveragedSource->SetNoDataValue(dfNoDataValueIn);
2452 : }
2453 : else
2454 : {
2455 1554 : poSimpleSource = new VRTSimpleSource();
2456 1554 : if (dfNoDataValueIn != VRT_NODATA_UNSET)
2457 0 : CPLError(
2458 : CE_Warning, CPLE_AppDefined,
2459 : "NODATA setting not currently supported for "
2460 : "neighbour sampled simple sources on Virtual Datasources.");
2461 : }
2462 :
2463 1554 : ConfigureSource(poSimpleSource, poSrcBand, FALSE, dfSrcXOff, dfSrcYOff,
2464 : dfSrcXSize, dfSrcYSize, dfDstXOff, dfDstYOff, dfDstXSize,
2465 : dfDstYSize);
2466 :
2467 : /* -------------------------------------------------------------------- */
2468 : /* add to list. */
2469 : /* -------------------------------------------------------------------- */
2470 1554 : return AddSource(poSimpleSource);
2471 : }
2472 :
2473 : /************************************************************************/
2474 : /* AddMaskBandSource() */
2475 : /************************************************************************/
2476 :
2477 : // poSrcBand is not the mask band, but the band from which the mask band is
2478 : // taken.
2479 18 : CPLErr VRTSourcedRasterBand::AddMaskBandSource(
2480 : GDALRasterBand *poSrcBand, double dfSrcXOff, double dfSrcYOff,
2481 : double dfSrcXSize, double dfSrcYSize, double dfDstXOff, double dfDstYOff,
2482 : double dfDstXSize, double dfDstYSize)
2483 : {
2484 : /* -------------------------------------------------------------------- */
2485 : /* Create source. */
2486 : /* -------------------------------------------------------------------- */
2487 18 : VRTSimpleSource *poSimpleSource = new VRTSimpleSource();
2488 :
2489 18 : ConfigureSource(poSimpleSource, poSrcBand, TRUE, dfSrcXOff, dfSrcYOff,
2490 : dfSrcXSize, dfSrcYSize, dfDstXOff, dfDstYOff, dfDstXSize,
2491 : dfDstYSize);
2492 :
2493 : /* -------------------------------------------------------------------- */
2494 : /* add to list. */
2495 : /* -------------------------------------------------------------------- */
2496 18 : return AddSource(poSimpleSource);
2497 : }
2498 :
2499 : /*! @endcond */
2500 :
2501 : /************************************************************************/
2502 : /* VRTAddSimpleSource() */
2503 : /************************************************************************/
2504 :
2505 : /**
2506 : * @see VRTSourcedRasterBand::AddSimpleSource().
2507 : */
2508 :
2509 960 : CPLErr CPL_STDCALL VRTAddSimpleSource(VRTSourcedRasterBandH hVRTBand,
2510 : GDALRasterBandH hSrcBand, int nSrcXOff,
2511 : int nSrcYOff, int nSrcXSize,
2512 : int nSrcYSize, int nDstXOff, int nDstYOff,
2513 : int nDstXSize, int nDstYSize,
2514 : const char *pszResampling,
2515 : double dfNoDataValue)
2516 : {
2517 960 : VALIDATE_POINTER1(hVRTBand, "VRTAddSimpleSource", CE_Failure);
2518 :
2519 960 : return reinterpret_cast<VRTSourcedRasterBand *>(hVRTBand)->AddSimpleSource(
2520 : reinterpret_cast<GDALRasterBand *>(hSrcBand), nSrcXOff, nSrcYOff,
2521 : nSrcXSize, nSrcYSize, nDstXOff, nDstYOff, nDstXSize, nDstYSize,
2522 960 : pszResampling, dfNoDataValue);
2523 : }
2524 :
2525 : /*! @cond Doxygen_Suppress */
2526 :
2527 : /************************************************************************/
2528 : /* AddComplexSource() */
2529 : /************************************************************************/
2530 :
2531 29 : CPLErr VRTSourcedRasterBand::AddComplexSource(
2532 : const char *pszFilename, int nBandIn, double dfSrcXOff, double dfSrcYOff,
2533 : double dfSrcXSize, double dfSrcYSize, double dfDstXOff, double dfDstYOff,
2534 : double dfDstXSize, double dfDstYSize, double dfScaleOff,
2535 : double dfScaleRatio, double dfNoDataValueIn, int nColorTableComponent)
2536 :
2537 : {
2538 : /* -------------------------------------------------------------------- */
2539 : /* Create source. */
2540 : /* -------------------------------------------------------------------- */
2541 29 : VRTComplexSource *const poSource = new VRTComplexSource();
2542 :
2543 29 : poSource->SetSrcBand(pszFilename, nBandIn);
2544 :
2545 29 : poSource->SetSrcWindow(dfSrcXOff, dfSrcYOff, dfSrcXSize, dfSrcYSize);
2546 29 : poSource->SetDstWindow(dfDstXOff, dfDstYOff, dfDstXSize, dfDstYSize);
2547 :
2548 : /* -------------------------------------------------------------------- */
2549 : /* Set complex parameters. */
2550 : /* -------------------------------------------------------------------- */
2551 29 : if (dfNoDataValueIn != VRT_NODATA_UNSET)
2552 8 : poSource->SetNoDataValue(dfNoDataValueIn);
2553 :
2554 29 : if (dfScaleOff != 0.0 || dfScaleRatio != 1.0)
2555 4 : poSource->SetLinearScaling(dfScaleOff, dfScaleRatio);
2556 :
2557 29 : poSource->SetColorTableComponent(nColorTableComponent);
2558 :
2559 : /* -------------------------------------------------------------------- */
2560 : /* add to list. */
2561 : /* -------------------------------------------------------------------- */
2562 29 : return AddSource(poSource);
2563 : }
2564 :
2565 : /************************************************************************/
2566 : /* AddComplexSource() */
2567 : /************************************************************************/
2568 :
2569 12 : CPLErr VRTSourcedRasterBand::AddComplexSource(
2570 : GDALRasterBand *poSrcBand, double dfSrcXOff, double dfSrcYOff,
2571 : double dfSrcXSize, double dfSrcYSize, double dfDstXOff, double dfDstYOff,
2572 : double dfDstXSize, double dfDstYSize, double dfScaleOff,
2573 : double dfScaleRatio, double dfNoDataValueIn, int nColorTableComponent)
2574 :
2575 : {
2576 : /* -------------------------------------------------------------------- */
2577 : /* Create source. */
2578 : /* -------------------------------------------------------------------- */
2579 12 : VRTComplexSource *const poSource = new VRTComplexSource();
2580 :
2581 12 : ConfigureSource(poSource, poSrcBand, FALSE, dfSrcXOff, dfSrcYOff,
2582 : dfSrcXSize, dfSrcYSize, dfDstXOff, dfDstYOff, dfDstXSize,
2583 : dfDstYSize);
2584 :
2585 : /* -------------------------------------------------------------------- */
2586 : /* Set complex parameters. */
2587 : /* -------------------------------------------------------------------- */
2588 12 : if (dfNoDataValueIn != VRT_NODATA_UNSET)
2589 0 : poSource->SetNoDataValue(dfNoDataValueIn);
2590 :
2591 12 : if (dfScaleOff != 0.0 || dfScaleRatio != 1.0)
2592 12 : poSource->SetLinearScaling(dfScaleOff, dfScaleRatio);
2593 :
2594 12 : poSource->SetColorTableComponent(nColorTableComponent);
2595 :
2596 : /* -------------------------------------------------------------------- */
2597 : /* add to list. */
2598 : /* -------------------------------------------------------------------- */
2599 12 : return AddSource(poSource);
2600 : }
2601 :
2602 : /*! @endcond */
2603 :
2604 : /************************************************************************/
2605 : /* VRTAddComplexSource() */
2606 : /************************************************************************/
2607 :
2608 : /**
2609 : * @see VRTSourcedRasterBand::AddComplexSource().
2610 : */
2611 :
2612 0 : CPLErr CPL_STDCALL VRTAddComplexSource(
2613 : VRTSourcedRasterBandH hVRTBand, GDALRasterBandH hSrcBand, int nSrcXOff,
2614 : int nSrcYOff, int nSrcXSize, int nSrcYSize, int nDstXOff, int nDstYOff,
2615 : int nDstXSize, int nDstYSize, double dfScaleOff, double dfScaleRatio,
2616 : double dfNoDataValue)
2617 : {
2618 0 : VALIDATE_POINTER1(hVRTBand, "VRTAddComplexSource", CE_Failure);
2619 :
2620 0 : return reinterpret_cast<VRTSourcedRasterBand *>(hVRTBand)->AddComplexSource(
2621 : reinterpret_cast<GDALRasterBand *>(hSrcBand), nSrcXOff, nSrcYOff,
2622 : nSrcXSize, nSrcYSize, nDstXOff, nDstYOff, nDstXSize, nDstYSize,
2623 0 : dfScaleOff, dfScaleRatio, dfNoDataValue);
2624 : }
2625 :
2626 : /*! @cond Doxygen_Suppress */
2627 :
2628 : /************************************************************************/
2629 : /* AddFuncSource() */
2630 : /************************************************************************/
2631 :
2632 12 : CPLErr VRTSourcedRasterBand::AddFuncSource(VRTImageReadFunc pfnReadFunc,
2633 : void *pCBData,
2634 : double dfNoDataValueIn)
2635 :
2636 : {
2637 : /* -------------------------------------------------------------------- */
2638 : /* Create source. */
2639 : /* -------------------------------------------------------------------- */
2640 12 : VRTFuncSource *const poFuncSource = new VRTFuncSource;
2641 :
2642 12 : poFuncSource->fNoDataValue = static_cast<float>(dfNoDataValueIn);
2643 12 : poFuncSource->pfnReadFunc = pfnReadFunc;
2644 12 : poFuncSource->pCBData = pCBData;
2645 12 : poFuncSource->eType = GetRasterDataType();
2646 :
2647 : /* -------------------------------------------------------------------- */
2648 : /* add to list. */
2649 : /* -------------------------------------------------------------------- */
2650 12 : return AddSource(poFuncSource);
2651 : }
2652 :
2653 : /*! @endcond */
2654 :
2655 : /************************************************************************/
2656 : /* VRTAddFuncSource() */
2657 : /************************************************************************/
2658 :
2659 : /**
2660 : * @see VRTSourcedRasterBand::AddFuncSource().
2661 : */
2662 :
2663 0 : CPLErr CPL_STDCALL VRTAddFuncSource(VRTSourcedRasterBandH hVRTBand,
2664 : VRTImageReadFunc pfnReadFunc, void *pCBData,
2665 : double dfNoDataValue)
2666 : {
2667 0 : VALIDATE_POINTER1(hVRTBand, "VRTAddFuncSource", CE_Failure);
2668 :
2669 0 : return reinterpret_cast<VRTSourcedRasterBand *>(hVRTBand)->AddFuncSource(
2670 0 : pfnReadFunc, pCBData, dfNoDataValue);
2671 : }
2672 :
2673 : /*! @cond Doxygen_Suppress */
2674 :
2675 : /************************************************************************/
2676 : /* GetMetadataDomainList() */
2677 : /************************************************************************/
2678 :
2679 0 : char **VRTSourcedRasterBand::GetMetadataDomainList()
2680 : {
2681 0 : return CSLAddString(GDALRasterBand::GetMetadataDomainList(),
2682 0 : "LocationInfo");
2683 : }
2684 :
2685 : /************************************************************************/
2686 : /* GetMetadataItem() */
2687 : /************************************************************************/
2688 :
2689 20414 : const char *VRTSourcedRasterBand::GetMetadataItem(const char *pszName,
2690 : const char *pszDomain)
2691 :
2692 : {
2693 : /* ==================================================================== */
2694 : /* LocationInfo handling. */
2695 : /* ==================================================================== */
2696 20414 : if (pszDomain != nullptr && EQUAL(pszDomain, "LocationInfo") &&
2697 1 : (STARTS_WITH_CI(pszName, "Pixel_") ||
2698 0 : STARTS_WITH_CI(pszName, "GeoPixel_")))
2699 : {
2700 : /* --------------------------------------------------------------------
2701 : */
2702 : /* What pixel are we aiming at? */
2703 : /* --------------------------------------------------------------------
2704 : */
2705 1 : int iPixel = 0;
2706 1 : int iLine = 0;
2707 :
2708 1 : if (STARTS_WITH_CI(pszName, "Pixel_"))
2709 : {
2710 : // TODO(schwehr): Replace sscanf.
2711 1 : if (sscanf(pszName + 6, "%d_%d", &iPixel, &iLine) != 2)
2712 0 : return nullptr;
2713 : }
2714 0 : else if (STARTS_WITH_CI(pszName, "GeoPixel_"))
2715 : {
2716 0 : const double dfGeoX = CPLAtof(pszName + 9);
2717 0 : const char *const pszUnderscore = strchr(pszName + 9, '_');
2718 0 : if (!pszUnderscore)
2719 0 : return nullptr;
2720 0 : const double dfGeoY = CPLAtof(pszUnderscore + 1);
2721 :
2722 0 : if (GetDataset() == nullptr)
2723 0 : return nullptr;
2724 :
2725 0 : double adfGeoTransform[6] = {0.0};
2726 0 : if (GetDataset()->GetGeoTransform(adfGeoTransform) != CE_None)
2727 0 : return nullptr;
2728 :
2729 0 : double adfInvGeoTransform[6] = {0.0};
2730 0 : if (!GDALInvGeoTransform(adfGeoTransform, adfInvGeoTransform))
2731 0 : return nullptr;
2732 :
2733 0 : iPixel = static_cast<int>(floor(adfInvGeoTransform[0] +
2734 0 : adfInvGeoTransform[1] * dfGeoX +
2735 0 : adfInvGeoTransform[2] * dfGeoY));
2736 0 : iLine = static_cast<int>(floor(adfInvGeoTransform[3] +
2737 0 : adfInvGeoTransform[4] * dfGeoX +
2738 0 : adfInvGeoTransform[5] * dfGeoY));
2739 : }
2740 : else
2741 : {
2742 0 : return nullptr;
2743 : }
2744 :
2745 2 : if (iPixel < 0 || iLine < 0 || iPixel >= GetXSize() ||
2746 1 : iLine >= GetYSize())
2747 0 : return nullptr;
2748 :
2749 : /* --------------------------------------------------------------------
2750 : */
2751 : /* Find the file(s) at this location. */
2752 : /* --------------------------------------------------------------------
2753 : */
2754 1 : char **papszFileList = nullptr;
2755 1 : int nListSize = 0; // keep it in this scope
2756 1 : int nListMaxSize = 0; // keep it in this scope
2757 : CPLHashSet *const hSetFiles =
2758 1 : CPLHashSetNew(CPLHashSetHashStr, CPLHashSetEqualStr, nullptr);
2759 :
2760 2 : for (int iSource = 0; iSource < nSources; iSource++)
2761 : {
2762 1 : if (!papoSources[iSource]->IsSimpleSource())
2763 0 : continue;
2764 :
2765 1 : VRTSimpleSource *const poSrc =
2766 1 : static_cast<VRTSimpleSource *>(papoSources[iSource]);
2767 :
2768 1 : double dfReqXOff = 0.0;
2769 1 : double dfReqYOff = 0.0;
2770 1 : double dfReqXSize = 0.0;
2771 1 : double dfReqYSize = 0.0;
2772 1 : int nReqXOff = 0;
2773 1 : int nReqYOff = 0;
2774 1 : int nReqXSize = 0;
2775 1 : int nReqYSize = 0;
2776 1 : int nOutXOff = 0;
2777 1 : int nOutYOff = 0;
2778 1 : int nOutXSize = 0;
2779 1 : int nOutYSize = 0;
2780 :
2781 1 : bool bError = false;
2782 1 : if (!poSrc->GetSrcDstWindow(iPixel, iLine, 1, 1, 1, 1, &dfReqXOff,
2783 : &dfReqYOff, &dfReqXSize, &dfReqYSize,
2784 : &nReqXOff, &nReqYOff, &nReqXSize,
2785 : &nReqYSize, &nOutXOff, &nOutYOff,
2786 : &nOutXSize, &nOutYSize, bError))
2787 : {
2788 0 : if (bError)
2789 : {
2790 0 : CSLDestroy(papszFileList);
2791 0 : CPLHashSetDestroy(hSetFiles);
2792 0 : return nullptr;
2793 : }
2794 0 : continue;
2795 : }
2796 :
2797 1 : poSrc->GetFileList(&papszFileList, &nListSize, &nListMaxSize,
2798 1 : hSetFiles);
2799 : }
2800 :
2801 : /* --------------------------------------------------------------------
2802 : */
2803 : /* Format into XML. */
2804 : /* --------------------------------------------------------------------
2805 : */
2806 1 : m_osLastLocationInfo = "<LocationInfo>";
2807 2 : for (int i = 0; i < nListSize && papszFileList[i] != nullptr; i++)
2808 : {
2809 1 : m_osLastLocationInfo += "<File>";
2810 : char *const pszXMLEscaped =
2811 1 : CPLEscapeString(papszFileList[i], -1, CPLES_XML);
2812 1 : m_osLastLocationInfo += pszXMLEscaped;
2813 1 : CPLFree(pszXMLEscaped);
2814 1 : m_osLastLocationInfo += "</File>";
2815 : }
2816 1 : m_osLastLocationInfo += "</LocationInfo>";
2817 :
2818 1 : CSLDestroy(papszFileList);
2819 1 : CPLHashSetDestroy(hSetFiles);
2820 :
2821 1 : return m_osLastLocationInfo.c_str();
2822 : }
2823 :
2824 : /* ==================================================================== */
2825 : /* Other domains. */
2826 : /* ==================================================================== */
2827 :
2828 20413 : return GDALRasterBand::GetMetadataItem(pszName, pszDomain);
2829 : }
2830 :
2831 : /************************************************************************/
2832 : /* GetMetadata() */
2833 : /************************************************************************/
2834 :
2835 13381 : char **VRTSourcedRasterBand::GetMetadata(const char *pszDomain)
2836 :
2837 : {
2838 : /* ==================================================================== */
2839 : /* vrt_sources domain handling. */
2840 : /* ==================================================================== */
2841 13381 : if (pszDomain != nullptr && EQUAL(pszDomain, "vrt_sources"))
2842 : {
2843 1 : CSLDestroy(m_papszSourceList);
2844 1 : m_papszSourceList = nullptr;
2845 :
2846 : /* --------------------------------------------------------------------
2847 : */
2848 : /* Process SimpleSources. */
2849 : /* --------------------------------------------------------------------
2850 : */
2851 2 : for (int iSource = 0; iSource < nSources; iSource++)
2852 : {
2853 : CPLXMLNode *const psXMLSrc =
2854 1 : papoSources[iSource]->SerializeToXML(nullptr);
2855 1 : if (psXMLSrc == nullptr)
2856 0 : continue;
2857 :
2858 1 : char *const pszXML = CPLSerializeXMLTree(psXMLSrc);
2859 :
2860 1 : m_papszSourceList = CSLSetNameValue(
2861 : m_papszSourceList, CPLSPrintf("source_%d", iSource), pszXML);
2862 1 : CPLFree(pszXML);
2863 1 : CPLDestroyXMLNode(psXMLSrc);
2864 : }
2865 :
2866 1 : return m_papszSourceList;
2867 : }
2868 :
2869 : /* ==================================================================== */
2870 : /* Other domains. */
2871 : /* ==================================================================== */
2872 :
2873 13380 : return GDALRasterBand::GetMetadata(pszDomain);
2874 : }
2875 :
2876 : /************************************************************************/
2877 : /* SetMetadataItem() */
2878 : /************************************************************************/
2879 :
2880 2340 : CPLErr VRTSourcedRasterBand::SetMetadataItem(const char *pszName,
2881 : const char *pszValue,
2882 : const char *pszDomain)
2883 :
2884 : {
2885 : #if DEBUG_VERBOSE
2886 : CPLDebug("VRT", "VRTSourcedRasterBand::SetMetadataItem(%s,%s,%s)\n",
2887 : pszName, pszValue ? pszValue : "(null)",
2888 : pszDomain ? pszDomain : "(null)");
2889 : #endif
2890 :
2891 2340 : if (pszDomain != nullptr && EQUAL(pszDomain, "new_vrt_sources"))
2892 : {
2893 : VRTDriver *const poDriver =
2894 1 : static_cast<VRTDriver *>(GDALGetDriverByName("VRT"));
2895 :
2896 1 : CPLXMLNode *const psTree = CPLParseXMLString(pszValue);
2897 1 : if (psTree == nullptr)
2898 0 : return CE_Failure;
2899 :
2900 1 : auto l_poDS = dynamic_cast<VRTDataset *>(GetDataset());
2901 1 : if (l_poDS == nullptr)
2902 : {
2903 0 : CPLDestroyXMLNode(psTree);
2904 0 : return CE_Failure;
2905 : }
2906 : VRTSource *const poSource =
2907 1 : poDriver->ParseSource(psTree, nullptr, l_poDS->m_oMapSharedSources);
2908 1 : CPLDestroyXMLNode(psTree);
2909 :
2910 1 : if (poSource != nullptr)
2911 1 : return AddSource(poSource);
2912 :
2913 0 : return CE_Failure;
2914 : }
2915 2339 : else if (pszDomain != nullptr && EQUAL(pszDomain, "vrt_sources"))
2916 : {
2917 1 : int iSource = 0;
2918 : // TODO(schwehr): Replace sscanf.
2919 2 : if (sscanf(pszName, "source_%d", &iSource) != 1 || iSource < 0 ||
2920 1 : iSource >= nSources)
2921 : {
2922 0 : CPLError(CE_Failure, CPLE_AppDefined,
2923 : "%s metadata item name is not recognized. "
2924 : "Should be between source_0 and source_%d",
2925 0 : pszName, nSources - 1);
2926 0 : return CE_Failure;
2927 : }
2928 :
2929 : VRTDriver *const poDriver =
2930 1 : static_cast<VRTDriver *>(GDALGetDriverByName("VRT"));
2931 :
2932 1 : CPLXMLNode *const psTree = CPLParseXMLString(pszValue);
2933 1 : if (psTree == nullptr)
2934 0 : return CE_Failure;
2935 :
2936 1 : auto l_poDS = dynamic_cast<VRTDataset *>(GetDataset());
2937 1 : if (l_poDS == nullptr)
2938 : {
2939 0 : CPLDestroyXMLNode(psTree);
2940 0 : return CE_Failure;
2941 : }
2942 : VRTSource *const poSource =
2943 1 : poDriver->ParseSource(psTree, nullptr, l_poDS->m_oMapSharedSources);
2944 1 : CPLDestroyXMLNode(psTree);
2945 :
2946 1 : if (poSource != nullptr)
2947 : {
2948 1 : delete papoSources[iSource];
2949 1 : papoSources[iSource] = poSource;
2950 1 : static_cast<VRTDataset *>(poDS)->SetNeedsFlush();
2951 1 : return CE_None;
2952 : }
2953 :
2954 0 : return CE_Failure;
2955 : }
2956 :
2957 2338 : return VRTRasterBand::SetMetadataItem(pszName, pszValue, pszDomain);
2958 : }
2959 :
2960 : /************************************************************************/
2961 : /* SetMetadata() */
2962 : /************************************************************************/
2963 :
2964 4381 : CPLErr VRTSourcedRasterBand::SetMetadata(char **papszNewMD,
2965 : const char *pszDomain)
2966 :
2967 : {
2968 4381 : if (pszDomain != nullptr && (EQUAL(pszDomain, "new_vrt_sources") ||
2969 4381 : EQUAL(pszDomain, "vrt_sources")))
2970 : {
2971 : VRTDriver *const poDriver =
2972 8 : static_cast<VRTDriver *>(GDALGetDriverByName("VRT"));
2973 :
2974 8 : if (EQUAL(pszDomain, "vrt_sources"))
2975 : {
2976 8 : for (int i = 0; i < nSources; i++)
2977 0 : delete papoSources[i];
2978 8 : CPLFree(papoSources);
2979 8 : papoSources = nullptr;
2980 8 : nSources = 0;
2981 : }
2982 :
2983 6 : for (const char *const pszMDItem :
2984 20 : cpl::Iterate(CSLConstList(papszNewMD)))
2985 : {
2986 8 : const char *const pszXML = CPLParseNameValue(pszMDItem, nullptr);
2987 8 : CPLXMLTreeCloser psTree(CPLParseXMLString(pszXML));
2988 8 : if (!psTree)
2989 0 : return CE_Failure;
2990 :
2991 8 : auto l_poDS = dynamic_cast<VRTDataset *>(GetDataset());
2992 8 : if (l_poDS == nullptr)
2993 : {
2994 0 : return CE_Failure;
2995 : }
2996 8 : VRTSource *const poSource = poDriver->ParseSource(
2997 8 : psTree.get(), nullptr, l_poDS->m_oMapSharedSources);
2998 8 : if (poSource == nullptr)
2999 2 : return CE_Failure;
3000 :
3001 6 : const CPLErr eErr = AddSource(poSource);
3002 : // cppcheck-suppress knownConditionTrueFalse
3003 6 : if (eErr != CE_None)
3004 0 : return eErr;
3005 : }
3006 :
3007 6 : return CE_None;
3008 : }
3009 :
3010 4373 : return VRTRasterBand::SetMetadata(papszNewMD, pszDomain);
3011 : }
3012 :
3013 : /************************************************************************/
3014 : /* GetFileList() */
3015 : /************************************************************************/
3016 :
3017 67 : void VRTSourcedRasterBand::GetFileList(char ***ppapszFileList, int *pnSize,
3018 : int *pnMaxSize, CPLHashSet *hSetFiles)
3019 : {
3020 178 : for (int i = 0; i < nSources; i++)
3021 : {
3022 111 : papoSources[i]->GetFileList(ppapszFileList, pnSize, pnMaxSize,
3023 111 : hSetFiles);
3024 : }
3025 :
3026 67 : VRTRasterBand::GetFileList(ppapszFileList, pnSize, pnMaxSize, hSetFiles);
3027 67 : }
3028 :
3029 : /************************************************************************/
3030 : /* CloseDependentDatasets() */
3031 : /************************************************************************/
3032 :
3033 7714 : int VRTSourcedRasterBand::CloseDependentDatasets()
3034 : {
3035 7714 : int ret = VRTRasterBand::CloseDependentDatasets();
3036 :
3037 7714 : if (nSources == 0)
3038 373 : return ret;
3039 :
3040 20282 : for (int i = 0; i < nSources; i++)
3041 12941 : delete papoSources[i];
3042 :
3043 7341 : CPLFree(papoSources);
3044 7341 : papoSources = nullptr;
3045 7341 : nSources = 0;
3046 :
3047 7341 : return TRUE;
3048 : }
3049 :
3050 : /************************************************************************/
3051 : /* FlushCache() */
3052 : /************************************************************************/
3053 :
3054 8119 : CPLErr VRTSourcedRasterBand::FlushCache(bool bAtClosing)
3055 : {
3056 8119 : CPLErr eErr = VRTRasterBand::FlushCache(bAtClosing);
3057 21494 : for (int i = 0; i < nSources && eErr == CE_None; i++)
3058 : {
3059 13375 : eErr = papoSources[i]->FlushCache(bAtClosing);
3060 : }
3061 8119 : return eErr;
3062 : }
3063 :
3064 : /************************************************************************/
3065 : /* RemoveCoveredSources() */
3066 : /************************************************************************/
3067 :
3068 : /** Remove sources that are covered by other sources.
3069 : *
3070 : * This method removes sources that are covered entirely by (one or several)
3071 : * sources of higher priority (even if they declare a nodata setting).
3072 : * This optimizes the size of the VRT and the rendering time.
3073 : */
3074 13 : void VRTSourcedRasterBand::RemoveCoveredSources(CSLConstList papszOptions)
3075 : {
3076 : #ifndef HAVE_GEOS
3077 : if (CPLTestBool(CSLFetchNameValueDef(
3078 : papszOptions, "EMIT_ERROR_IF_GEOS_NOT_AVAILABLE", "TRUE")))
3079 : {
3080 : CPLError(CE_Failure, CPLE_NotSupported,
3081 : "RemoveCoveredSources() not implemented in builds "
3082 : "without GEOS support");
3083 : }
3084 : #else
3085 : (void)papszOptions;
3086 :
3087 : CPLRectObj globalBounds;
3088 13 : globalBounds.minx = 0;
3089 13 : globalBounds.miny = 0;
3090 13 : globalBounds.maxx = nRasterXSize;
3091 13 : globalBounds.maxy = nRasterYSize;
3092 :
3093 : // Create an index with the bbox of all sources
3094 13 : CPLQuadTree *hTree = CPLQuadTreeCreate(&globalBounds, nullptr);
3095 36 : for (int i = 0; i < nSources; i++)
3096 : {
3097 23 : if (papoSources[i]->IsSimpleSource())
3098 : {
3099 : VRTSimpleSource *poSS =
3100 23 : cpl::down_cast<VRTSimpleSource *>(papoSources[i]);
3101 23 : void *hFeature =
3102 23 : reinterpret_cast<void *>(static_cast<uintptr_t>(i));
3103 : CPLRectObj rect;
3104 23 : rect.minx = std::max(0.0, poSS->m_dfDstXOff);
3105 23 : rect.miny = std::max(0.0, poSS->m_dfDstYOff);
3106 46 : rect.maxx = std::min(double(nRasterXSize),
3107 23 : poSS->m_dfDstXOff + poSS->m_dfDstXSize);
3108 46 : rect.maxy = std::min(double(nRasterYSize),
3109 23 : poSS->m_dfDstYOff + poSS->m_dfDstYSize);
3110 23 : CPLQuadTreeInsertWithBounds(hTree, hFeature, &rect);
3111 : }
3112 : }
3113 :
3114 36 : for (int i = 0; i < nSources; i++)
3115 : {
3116 23 : if (papoSources[i]->IsSimpleSource())
3117 : {
3118 : VRTSimpleSource *poSS =
3119 23 : cpl::down_cast<VRTSimpleSource *>(papoSources[i]);
3120 : CPLRectObj rect;
3121 23 : rect.minx = std::max(0.0, poSS->m_dfDstXOff);
3122 23 : rect.miny = std::max(0.0, poSS->m_dfDstYOff);
3123 46 : rect.maxx = std::min(double(nRasterXSize),
3124 23 : poSS->m_dfDstXOff + poSS->m_dfDstXSize);
3125 46 : rect.maxy = std::min(double(nRasterYSize),
3126 23 : poSS->m_dfDstYOff + poSS->m_dfDstYSize);
3127 :
3128 : // Find sources whose extent intersect with the current one
3129 23 : int nFeatureCount = 0;
3130 : void **pahFeatures =
3131 23 : CPLQuadTreeSearch(hTree, &rect, &nFeatureCount);
3132 :
3133 : // Compute the bounding box of those sources, only if they are
3134 : // on top of the current one
3135 : CPLRectObj rectIntersecting;
3136 23 : rectIntersecting.minx = std::numeric_limits<double>::max();
3137 23 : rectIntersecting.miny = std::numeric_limits<double>::max();
3138 23 : rectIntersecting.maxx = -std::numeric_limits<double>::max();
3139 23 : rectIntersecting.maxy = -std::numeric_limits<double>::max();
3140 59 : for (int j = 0; j < nFeatureCount; j++)
3141 : {
3142 36 : const int curFeature = static_cast<int>(
3143 36 : reinterpret_cast<uintptr_t>(pahFeatures[j]));
3144 36 : if (curFeature > i)
3145 : {
3146 : VRTSimpleSource *poOtherSS =
3147 26 : cpl::down_cast<VRTSimpleSource *>(
3148 13 : papoSources[curFeature]);
3149 13 : rectIntersecting.minx =
3150 13 : std::min(rectIntersecting.minx, poOtherSS->m_dfDstXOff);
3151 13 : rectIntersecting.miny =
3152 13 : std::min(rectIntersecting.miny, poOtherSS->m_dfDstYOff);
3153 13 : rectIntersecting.maxx = std::max(
3154 : rectIntersecting.maxx,
3155 13 : poOtherSS->m_dfDstXOff + poOtherSS->m_dfDstXSize);
3156 13 : rectIntersecting.maxy = std::max(
3157 : rectIntersecting.maxy,
3158 13 : poOtherSS->m_dfDstYOff + poOtherSS->m_dfDstXSize);
3159 : }
3160 : }
3161 :
3162 : // If the boundinx box of those sources overlap the current one,
3163 : // then compute their union, and check if it contains the current
3164 : // source
3165 23 : if (rectIntersecting.minx <= rect.minx &&
3166 7 : rectIntersecting.miny <= rect.miny &&
3167 7 : rectIntersecting.maxx >= rect.maxx &&
3168 7 : rectIntersecting.maxy >= rect.maxy)
3169 : {
3170 14 : OGRPolygon oPoly;
3171 : {
3172 7 : auto poLR = new OGRLinearRing();
3173 7 : poLR->addPoint(rect.minx, rect.miny);
3174 7 : poLR->addPoint(rect.minx, rect.maxy);
3175 7 : poLR->addPoint(rect.maxx, rect.maxy);
3176 7 : poLR->addPoint(rect.maxx, rect.miny);
3177 7 : poLR->addPoint(rect.minx, rect.miny);
3178 7 : oPoly.addRingDirectly(poLR);
3179 : }
3180 :
3181 7 : std::unique_ptr<OGRGeometry> poUnion;
3182 24 : for (int j = 0; j < nFeatureCount; j++)
3183 : {
3184 17 : const int curFeature = static_cast<int>(
3185 17 : reinterpret_cast<uintptr_t>(pahFeatures[j]));
3186 17 : if (curFeature > i)
3187 : {
3188 : VRTSimpleSource *poOtherSS =
3189 20 : cpl::down_cast<VRTSimpleSource *>(
3190 10 : papoSources[curFeature]);
3191 : CPLRectObj otherRect;
3192 10 : otherRect.minx = std::max(0.0, poOtherSS->m_dfDstXOff);
3193 10 : otherRect.miny = std::max(0.0, poOtherSS->m_dfDstYOff);
3194 20 : otherRect.maxx = std::min(double(nRasterXSize),
3195 20 : poOtherSS->m_dfDstXOff +
3196 10 : poOtherSS->m_dfDstXSize);
3197 20 : otherRect.maxy = std::min(double(nRasterYSize),
3198 20 : poOtherSS->m_dfDstYOff +
3199 10 : poOtherSS->m_dfDstYSize);
3200 20 : OGRPolygon oOtherPoly;
3201 : {
3202 10 : auto poLR = new OGRLinearRing();
3203 10 : poLR->addPoint(otherRect.minx, otherRect.miny);
3204 10 : poLR->addPoint(otherRect.minx, otherRect.maxy);
3205 10 : poLR->addPoint(otherRect.maxx, otherRect.maxy);
3206 10 : poLR->addPoint(otherRect.maxx, otherRect.miny);
3207 10 : poLR->addPoint(otherRect.minx, otherRect.miny);
3208 10 : oOtherPoly.addRingDirectly(poLR);
3209 : }
3210 10 : if (poUnion == nullptr)
3211 7 : poUnion.reset(oOtherPoly.clone());
3212 : else
3213 3 : poUnion.reset(oOtherPoly.Union(poUnion.get()));
3214 : }
3215 : }
3216 :
3217 7 : if (poUnion != nullptr && poUnion->Contains(&oPoly))
3218 : {
3219 : // We can remove the current source
3220 7 : delete papoSources[i];
3221 7 : papoSources[i] = nullptr;
3222 : }
3223 : }
3224 23 : CPLFree(pahFeatures);
3225 :
3226 23 : void *hFeature =
3227 23 : reinterpret_cast<void *>(static_cast<uintptr_t>(i));
3228 23 : CPLQuadTreeRemove(hTree, hFeature, &rect);
3229 : }
3230 : }
3231 :
3232 : // Compact the papoSources array
3233 13 : int iDst = 0;
3234 36 : for (int iSrc = 0; iSrc < nSources; iSrc++)
3235 : {
3236 23 : if (papoSources[iSrc])
3237 16 : papoSources[iDst++] = papoSources[iSrc];
3238 : }
3239 13 : nSources = iDst;
3240 :
3241 13 : CPLQuadTreeDestroy(hTree);
3242 : #endif
3243 13 : }
3244 :
3245 : /*! @endcond */
|
