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