Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Virtual GDAL Datasets
4 : * Purpose: Implementation of VRTWarpedRasterBand *and VRTWarpedDataset.
5 : * Author: Frank Warmerdam <warmerdam@pobox.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2004, 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 "vrtdataset.h"
16 :
17 : #include <stddef.h>
18 : #include <stdio.h>
19 : #include <stdlib.h>
20 : #include <string.h>
21 : #include <algorithm>
22 : #include <map>
23 :
24 : // Suppress deprecation warning for GDALOpenVerticalShiftGrid and
25 : // GDALApplyVerticalShiftGrid
26 : #ifndef CPL_WARN_DEPRECATED_GDALOpenVerticalShiftGrid
27 : #define CPL_WARN_DEPRECATED_GDALOpenVerticalShiftGrid(x)
28 : #define CPL_WARN_DEPRECATED_GDALApplyVerticalShiftGrid(x)
29 : #endif
30 :
31 : #include "cpl_conv.h"
32 : #include "cpl_error.h"
33 : #include "cpl_minixml.h"
34 : #include "cpl_progress.h"
35 : #include "cpl_string.h"
36 : #include "cpl_vsi.h"
37 : #include "gdal.h"
38 : #include "gdal_alg.h"
39 : #include "gdal_alg_priv.h"
40 : #include "gdal_priv.h"
41 : #include "gdalwarper.h"
42 : #include "ogr_geometry.h"
43 :
44 : /************************************************************************/
45 : /* GDALAutoCreateWarpedVRT() */
46 : /************************************************************************/
47 :
48 : /**
49 : * Create virtual warped dataset automatically.
50 : *
51 : * This function will create a warped virtual file representing the
52 : * input image warped into the target coordinate system. A GenImgProj
53 : * transformation is created to accomplish any required GCP/Geotransform
54 : * warp and reprojection to the target coordinate system. The output virtual
55 : * dataset will be "northup" in the target coordinate system. The
56 : * GDALSuggestedWarpOutput() function is used to determine the bounds and
57 : * resolution of the output virtual file which should be large enough to
58 : * include all the input image
59 : *
60 : * If you want to create an alpha band if the source dataset has none, set
61 : * psOptionsIn->nDstAlphaBand = GDALGetRasterCount(hSrcDS) + 1.
62 : *
63 : * Note that the constructed GDALDatasetH will acquire one or more references
64 : * to the passed in hSrcDS. Reference counting semantics on the source
65 : * dataset should be honoured. That is, don't just GDALClose() it unless it
66 : * was opened with GDALOpenShared().
67 : *
68 : * It is possible to "transfer" the ownership of the source dataset
69 : * to the warped dataset in the following way:
70 : *
71 : * \code{.c}
72 : * GDALDatasetH src_ds = GDALOpen("source.tif");
73 : * GDALDatasetH warped_ds = GDALAutoCreateWarpedVRT( src_ds, ... );
74 : * GDALReleaseDataset(src_ds); // src_ds is not "owned" fully by warped_ds.
75 : * Do NOT use GDALClose(src_ds) here
76 : * ...
77 : * ...
78 : * GDALReleaseDataset(warped_ds); // or GDALClose(warped_ds);
79 : * \endcode
80 : *
81 : * Traditional nested calls are also possible of course:
82 : *
83 : * \code{.c}
84 : * GDALDatasetH src_ds = GDALOpen("source.tif");
85 : * GDALDatasetH warped_ds = GDALAutoCreateWarpedVRT( src_ds, ... );
86 : * ...
87 : * ...
88 : * GDALReleaseDataset(warped_ds); // or GDALClose(warped_ds);
89 : * GDALReleaseDataset(src_ds); // or GDALClose(src_ds);
90 : * \endcode
91 : *
92 : * The returned dataset will have no associated filename for itself. If you
93 : * want to write the virtual dataset description to a file, use the
94 : * GDALSetDescription() function (or SetDescription() method) on the dataset
95 : * to assign a filename before it is closed.
96 : *
97 : * @param hSrcDS The source dataset.
98 : *
99 : * @param pszSrcWKT The coordinate system of the source image. If NULL, it
100 : * will be read from the source image.
101 : *
102 : * @param pszDstWKT The coordinate system to convert to. If NULL no change
103 : * of coordinate system will take place.
104 : *
105 : * @param eResampleAlg One of GRA_NearestNeighbour, GRA_Bilinear, GRA_Cubic,
106 : * GRA_CubicSpline, GRA_Lanczos, GRA_Average, GRA_RMS or GRA_Mode.
107 : * Controls the sampling method used.
108 : *
109 : * @param dfMaxError Maximum error measured in input pixels that is allowed in
110 : * approximating the transformation (0.0 for exact calculations).
111 : *
112 : * @param psOptionsIn Additional warp options, normally NULL.
113 : *
114 : * @return NULL on failure, or a new virtual dataset handle on success.
115 : */
116 :
117 : GDALDatasetH CPL_STDCALL
118 23 : GDALAutoCreateWarpedVRT(GDALDatasetH hSrcDS, const char *pszSrcWKT,
119 : const char *pszDstWKT, GDALResampleAlg eResampleAlg,
120 : double dfMaxError, const GDALWarpOptions *psOptionsIn)
121 : {
122 23 : return GDALAutoCreateWarpedVRTEx(hSrcDS, pszSrcWKT, pszDstWKT, eResampleAlg,
123 23 : dfMaxError, psOptionsIn, nullptr);
124 : }
125 :
126 : /************************************************************************/
127 : /* GDALAutoCreateWarpedVRTEx() */
128 : /************************************************************************/
129 :
130 : /**
131 : * Create virtual warped dataset automatically.
132 : *
133 : * Compared to GDALAutoCreateWarpedVRT() this function adds one extra
134 : * argument: options to be passed to GDALCreateGenImgProjTransformer2().
135 : *
136 : * @since 3.2
137 : */
138 :
139 23 : GDALDatasetH CPL_STDCALL GDALAutoCreateWarpedVRTEx(
140 : GDALDatasetH hSrcDS, const char *pszSrcWKT, const char *pszDstWKT,
141 : GDALResampleAlg eResampleAlg, double dfMaxError,
142 : const GDALWarpOptions *psOptionsIn, CSLConstList papszTransformerOptions)
143 : {
144 23 : VALIDATE_POINTER1(hSrcDS, "GDALAutoCreateWarpedVRT", nullptr);
145 :
146 : /* -------------------------------------------------------------------- */
147 : /* Populate the warp options. */
148 : /* -------------------------------------------------------------------- */
149 23 : GDALWarpOptions *psWO = nullptr;
150 23 : if (psOptionsIn != nullptr)
151 1 : psWO = GDALCloneWarpOptions(psOptionsIn);
152 : else
153 22 : psWO = GDALCreateWarpOptions();
154 :
155 23 : psWO->eResampleAlg = eResampleAlg;
156 :
157 23 : psWO->hSrcDS = hSrcDS;
158 :
159 23 : GDALWarpInitDefaultBandMapping(psWO, GDALGetRasterCount(hSrcDS));
160 :
161 : /* -------------------------------------------------------------------- */
162 : /* Setup no data values (if not done in psOptionsIn) */
163 : /* -------------------------------------------------------------------- */
164 23 : if (psWO->padfSrcNoDataReal == nullptr &&
165 23 : psWO->padfDstNoDataReal == nullptr && psWO->nSrcAlphaBand == 0)
166 : {
167 : // If none of the provided input nodata values can be represented in the
168 : // data type of the corresponding source band, ignore them.
169 23 : int nCountInvalidSrcNoDataReal = 0;
170 48 : for (int i = 0; i < psWO->nBandCount; i++)
171 : {
172 : GDALRasterBandH rasterBand =
173 25 : GDALGetRasterBand(psWO->hSrcDS, psWO->panSrcBands[i]);
174 :
175 : int hasNoDataValue;
176 : double noDataValue =
177 25 : GDALGetRasterNoDataValue(rasterBand, &hasNoDataValue);
178 :
179 25 : if (hasNoDataValue &&
180 3 : !GDALIsValueExactAs(noDataValue,
181 25 : GDALGetRasterDataType(rasterBand)))
182 : {
183 2 : nCountInvalidSrcNoDataReal++;
184 : }
185 : }
186 :
187 23 : if (nCountInvalidSrcNoDataReal != psWO->nBandCount)
188 : {
189 44 : for (int i = 0; i < psWO->nBandCount; i++)
190 : {
191 : GDALRasterBandH rasterBand =
192 23 : GDALGetRasterBand(psWO->hSrcDS, psWO->panSrcBands[i]);
193 :
194 : int hasNoDataValue;
195 : double noDataValue =
196 23 : GDALGetRasterNoDataValue(rasterBand, &hasNoDataValue);
197 :
198 23 : if (hasNoDataValue)
199 : {
200 : // Check if the nodata value is out of range
201 1 : int bClamped = FALSE;
202 1 : int bRounded = FALSE;
203 1 : CPL_IGNORE_RET_VAL(GDALAdjustValueToDataType(
204 : GDALGetRasterDataType(rasterBand), noDataValue,
205 : &bClamped, &bRounded));
206 1 : if (!bClamped)
207 : {
208 1 : GDALWarpInitNoDataReal(psWO, -1e10);
209 1 : if (psWO->padfSrcNoDataReal != nullptr &&
210 1 : psWO->padfDstNoDataReal != nullptr)
211 : {
212 1 : psWO->padfSrcNoDataReal[i] = noDataValue;
213 1 : psWO->padfDstNoDataReal[i] = noDataValue;
214 : }
215 : }
216 : }
217 : }
218 : }
219 :
220 23 : if (psWO->padfDstNoDataReal != nullptr)
221 : {
222 1 : if (CSLFetchNameValue(psWO->papszWarpOptions, "INIT_DEST") ==
223 : nullptr)
224 : {
225 1 : psWO->papszWarpOptions = CSLSetNameValue(
226 : psWO->papszWarpOptions, "INIT_DEST", "NO_DATA");
227 : }
228 : }
229 : }
230 :
231 : /* -------------------------------------------------------------------- */
232 : /* Create the transformer. */
233 : /* -------------------------------------------------------------------- */
234 23 : psWO->pfnTransformer = GDALGenImgProjTransform;
235 :
236 23 : char **papszOptions = nullptr;
237 23 : if (pszSrcWKT != nullptr)
238 2 : papszOptions = CSLSetNameValue(papszOptions, "SRC_SRS", pszSrcWKT);
239 23 : if (pszDstWKT != nullptr)
240 6 : papszOptions = CSLSetNameValue(papszOptions, "DST_SRS", pszDstWKT);
241 23 : papszOptions = CSLMerge(papszOptions, papszTransformerOptions);
242 23 : psWO->pTransformerArg =
243 23 : GDALCreateGenImgProjTransformer2(psWO->hSrcDS, nullptr, papszOptions);
244 23 : CSLDestroy(papszOptions);
245 :
246 23 : if (psWO->pTransformerArg == nullptr)
247 : {
248 1 : GDALDestroyWarpOptions(psWO);
249 1 : return nullptr;
250 : }
251 :
252 : /* -------------------------------------------------------------------- */
253 : /* Figure out the desired output bounds and resolution. */
254 : /* -------------------------------------------------------------------- */
255 22 : double adfDstGeoTransform[6] = {0.0};
256 22 : int nDstPixels = 0;
257 22 : int nDstLines = 0;
258 22 : CPLErr eErr = GDALSuggestedWarpOutput(
259 : hSrcDS, psWO->pfnTransformer, psWO->pTransformerArg, adfDstGeoTransform,
260 : &nDstPixels, &nDstLines);
261 22 : if (eErr != CE_None)
262 : {
263 0 : GDALDestroyTransformer(psWO->pTransformerArg);
264 0 : GDALDestroyWarpOptions(psWO);
265 0 : return nullptr;
266 : }
267 :
268 22 : if (nDstPixels <= 0 || nDstLines <= 0)
269 : {
270 : // Can happen when the extent of the source in the target CRS is
271 : // degenerate, for example when the source crosses the antimeridian
272 : // and the suggested longitude span covers the whole globe, so that
273 : // the suggested latitude span rounds down to zero line.
274 1 : CPLError(CE_Failure, CPLE_AppDefined,
275 : "Computed output size is %d x %d. The extent of the source "
276 : "in the target CRS is degenerate.",
277 : nDstPixels, nDstLines);
278 1 : GDALDestroyTransformer(psWO->pTransformerArg);
279 1 : GDALDestroyWarpOptions(psWO);
280 1 : return nullptr;
281 : }
282 :
283 : /* -------------------------------------------------------------------- */
284 : /* Update the transformer to include an output geotransform */
285 : /* back to pixel/line coordinates. */
286 : /* */
287 : /* -------------------------------------------------------------------- */
288 21 : GDALSetGenImgProjTransformerDstGeoTransform(psWO->pTransformerArg,
289 : adfDstGeoTransform);
290 :
291 : /* -------------------------------------------------------------------- */
292 : /* Do we want to apply an approximating transformation? */
293 : /* -------------------------------------------------------------------- */
294 21 : if (dfMaxError > 0.0)
295 : {
296 3 : psWO->pTransformerArg = GDALCreateApproxTransformer(
297 : psWO->pfnTransformer, psWO->pTransformerArg, dfMaxError);
298 3 : psWO->pfnTransformer = GDALApproxTransform;
299 3 : GDALApproxTransformerOwnsSubtransformer(psWO->pTransformerArg, TRUE);
300 : }
301 :
302 : /* -------------------------------------------------------------------- */
303 : /* Create the VRT file. */
304 : /* -------------------------------------------------------------------- */
305 21 : GDALDatasetH hDstDS = GDALCreateWarpedVRT(hSrcDS, nDstPixels, nDstLines,
306 : adfDstGeoTransform, psWO);
307 :
308 21 : GDALDestroyWarpOptions(psWO);
309 :
310 21 : if (hDstDS != nullptr)
311 : {
312 21 : if (pszDstWKT != nullptr)
313 4 : GDALSetProjection(hDstDS, pszDstWKT);
314 17 : else if (pszSrcWKT != nullptr)
315 0 : GDALSetProjection(hDstDS, pszSrcWKT);
316 17 : else if (GDALGetGCPCount(hSrcDS) > 0)
317 2 : GDALSetProjection(hDstDS, GDALGetGCPProjection(hSrcDS));
318 : else
319 15 : GDALSetProjection(hDstDS, GDALGetProjectionRef(hSrcDS));
320 : }
321 :
322 21 : return hDstDS;
323 : }
324 :
325 : /************************************************************************/
326 : /* GDALCreateWarpedVRT() */
327 : /************************************************************************/
328 :
329 : /**
330 : * Create virtual warped dataset.
331 : *
332 : * This function will create a warped virtual file representing the
333 : * input image warped based on a provided transformation. Output bounds
334 : * and resolution are provided explicitly.
335 : *
336 : * If you want to create an alpha band if the source dataset has none, set
337 : * psOptions->nDstAlphaBand = GDALGetRasterCount(hSrcDS) + 1.
338 : *
339 : * Note that the constructed GDALDatasetH will acquire one or more references
340 : * to the passed in hSrcDS. Reference counting semantics on the source
341 : * dataset should be honoured. That is, don't just GDALClose() it unless it
342 : * was opened with GDALOpenShared().
343 : *
344 : * It is possible to "transfer" the ownership of the source dataset
345 : * to the warped dataset in the following way:
346 : *
347 : * \code{.c}
348 : * GDALDatasetH src_ds = GDALOpen("source.tif");
349 : * GDALDatasetH warped_ds = GDALAutoCreateWarpedVRT( src_ds, ... );
350 : * GDALReleaseDataset(src_ds); // src_ds is not "owned" fully by warped_ds.
351 : * Do NOT use GDALClose(src_ds) here
352 : * ...
353 : * ...
354 : * GDALReleaseDataset(warped_ds); // or GDALClose(warped_ds);
355 : * \endcode
356 : *
357 : * Traditional nested calls are also possible of course:
358 : *
359 : * \code{.c}
360 : * GDALDatasetH src_ds = GDALOpen("source.tif");
361 : * GDALDatasetH warped_ds = GDALAutoCreateWarpedVRT( src_ds, ... );
362 : * ...
363 : * ...
364 : * GDALReleaseDataset(warped_ds); // or GDALClose(warped_ds);
365 : * GDALReleaseDataset(src_ds); // or GDALClose(src_ds);
366 : * \endcode
367 : *
368 : * The returned dataset will have no associated filename for itself. If you
369 : * want to write the virtual dataset description to a file, use the
370 : * GDALSetDescription() function (or SetDescription() method) on the dataset
371 : * to assign a filename before it is closed.
372 : *
373 : * @param hSrcDS The source dataset.
374 : *
375 : * @param nPixels Width of the virtual warped dataset to create
376 : *
377 : * @param nLines Height of the virtual warped dataset to create
378 : *
379 : * @param padfGeoTransform Geotransform matrix of the virtual warped dataset
380 : * to create
381 : *
382 : * @param psOptions Warp options. Must be different from NULL.
383 : *
384 : * @return NULL on failure, or a new virtual dataset handle on success.
385 : */
386 :
387 44 : GDALDatasetH CPL_STDCALL GDALCreateWarpedVRT(GDALDatasetH hSrcDS, int nPixels,
388 : int nLines,
389 : const double *padfGeoTransform,
390 : GDALWarpOptions *psOptions)
391 :
392 : {
393 44 : VALIDATE_POINTER1(hSrcDS, "GDALCreateWarpedVRT", nullptr);
394 44 : VALIDATE_POINTER1(psOptions, "GDALCreateWarpedVRT", nullptr);
395 :
396 44 : if (nPixels <= 0 || nLines <= 0)
397 : {
398 0 : CPLError(CE_Failure, CPLE_IllegalArg,
399 : "Attempt to create %dx%d warped VRT is illegal, sizes must "
400 : "be larger than zero.",
401 : nPixels, nLines);
402 0 : return nullptr;
403 : }
404 :
405 : /* -------------------------------------------------------------------- */
406 : /* Create the VRTDataset and populate it with bands. */
407 : /* -------------------------------------------------------------------- */
408 44 : VRTWarpedDataset *poDS = new VRTWarpedDataset(nPixels, nLines);
409 :
410 : // Call this before assigning hDstDS
411 44 : GDALWarpResolveWorkingDataType(psOptions);
412 :
413 44 : psOptions->hDstDS = poDS;
414 44 : poDS->SetGeoTransform(GDALGeoTransform(padfGeoTransform));
415 :
416 94 : for (int i = 0; i < psOptions->nBandCount; i++)
417 : {
418 50 : int nDstBand = psOptions->panDstBands[i];
419 100 : while (poDS->GetRasterCount() < nDstBand)
420 : {
421 50 : poDS->AddBand(psOptions->eWorkingDataType, nullptr);
422 : }
423 :
424 : VRTWarpedRasterBand *poBand =
425 50 : static_cast<VRTWarpedRasterBand *>(poDS->GetRasterBand(nDstBand));
426 : GDALRasterBand *poSrcBand = static_cast<GDALRasterBand *>(
427 50 : GDALGetRasterBand(hSrcDS, psOptions->panSrcBands[i]));
428 :
429 50 : poBand->CopyCommonInfoFrom(poSrcBand);
430 : }
431 :
432 45 : while (poDS->GetRasterCount() < psOptions->nDstAlphaBand)
433 : {
434 1 : poDS->AddBand(psOptions->eWorkingDataType, nullptr);
435 : }
436 44 : if (psOptions->nDstAlphaBand)
437 : {
438 1 : poDS->GetRasterBand(psOptions->nDstAlphaBand)
439 1 : ->SetColorInterpretation(GCI_AlphaBand);
440 : }
441 :
442 : /* -------------------------------------------------------------------- */
443 : /* Initialize the warp on the VRTWarpedDataset. */
444 : /* -------------------------------------------------------------------- */
445 44 : const CPLErr eErr = poDS->Initialize(psOptions);
446 44 : if (eErr == CE_Failure)
447 : {
448 0 : psOptions->hDstDS = nullptr;
449 0 : delete poDS;
450 0 : return nullptr;
451 : }
452 :
453 44 : return poDS;
454 : }
455 :
456 : /*! @cond Doxygen_Suppress */
457 :
458 : /************************************************************************/
459 : /* ==================================================================== */
460 : /* VRTWarpedDataset */
461 : /* ==================================================================== */
462 : /************************************************************************/
463 :
464 : /************************************************************************/
465 : /* VRTWarpedDataset() */
466 : /************************************************************************/
467 :
468 392 : VRTWarpedDataset::VRTWarpedDataset(int nXSize, int nYSize, int nBlockXSize,
469 392 : int nBlockYSize)
470 : : VRTDataset(nXSize, nYSize,
471 391 : nBlockXSize > 0 ? nBlockXSize : std::min(nXSize, 512),
472 391 : nBlockYSize > 0 ? nBlockYSize : std::min(nYSize, 128)),
473 1174 : m_poWarper(nullptr), m_nSrcOvrLevel(-2)
474 : {
475 392 : eAccess = GA_Update;
476 392 : DisableReadWriteMutex();
477 392 : }
478 :
479 : /************************************************************************/
480 : /* ~VRTWarpedDataset() */
481 : /************************************************************************/
482 :
483 784 : VRTWarpedDataset::~VRTWarpedDataset()
484 :
485 : {
486 392 : VRTWarpedDataset::FlushCache(true);
487 392 : VRTWarpedDataset::CloseDependentDatasets();
488 784 : }
489 :
490 : /************************************************************************/
491 : /* CloseDependentDatasets() */
492 : /************************************************************************/
493 :
494 405 : int VRTWarpedDataset::CloseDependentDatasets()
495 : {
496 405 : bool bHasDroppedRef = CPL_TO_BOOL(VRTDataset::CloseDependentDatasets());
497 :
498 : /* -------------------------------------------------------------------- */
499 : /* Cleanup overviews. */
500 : /* -------------------------------------------------------------------- */
501 430 : for (auto &poDS : m_apoOverviews)
502 : {
503 25 : if (poDS && poDS->Release())
504 : {
505 0 : bHasDroppedRef = true;
506 : }
507 : }
508 :
509 405 : m_apoOverviews.clear();
510 :
511 : /* -------------------------------------------------------------------- */
512 : /* Cleanup warper if one is in effect. */
513 : /* -------------------------------------------------------------------- */
514 405 : if (m_poWarper != nullptr)
515 : {
516 390 : const GDALWarpOptions *psWO = m_poWarper->GetOptions();
517 :
518 : /* --------------------------------------------------------------------
519 : */
520 : /* We take care to only call GDALClose() on psWO->hSrcDS if the */
521 : /* reference count drops to zero. This is makes it so that we */
522 : /* can operate reference counting semantics more-or-less */
523 : /* properly even if the dataset isn't open in shared mode, */
524 : /* though we require that the caller also honour the reference */
525 : /* counting semantics even though it isn't a shared dataset. */
526 : /* --------------------------------------------------------------------
527 : */
528 390 : if (psWO != nullptr && psWO->hSrcDS != nullptr)
529 : {
530 389 : if (GDALReleaseDataset(psWO->hSrcDS))
531 : {
532 264 : bHasDroppedRef = true;
533 : }
534 : }
535 :
536 : /* --------------------------------------------------------------------
537 : */
538 : /* We are responsible for cleaning up the transformer ourselves. */
539 : /* --------------------------------------------------------------------
540 : */
541 390 : if (psWO != nullptr && psWO->pTransformerArg != nullptr)
542 389 : GDALDestroyTransformer(psWO->pTransformerArg);
543 :
544 390 : delete m_poWarper;
545 390 : m_poWarper = nullptr;
546 : }
547 :
548 : /* -------------------------------------------------------------------- */
549 : /* Destroy the raster bands if they exist. */
550 : /* -------------------------------------------------------------------- */
551 1060 : for (int iBand = 0; iBand < nBands; iBand++)
552 : {
553 655 : delete papoBands[iBand];
554 : }
555 405 : nBands = 0;
556 :
557 405 : return bHasDroppedRef;
558 : }
559 :
560 : /************************************************************************/
561 : /* SetSrcOverviewLevel() */
562 : /************************************************************************/
563 :
564 186 : CPLErr VRTWarpedDataset::SetMetadataItem(const char *pszName,
565 : const char *pszValue,
566 : const char *pszDomain)
567 :
568 : {
569 186 : if ((pszDomain == nullptr || EQUAL(pszDomain, "")) &&
570 186 : EQUAL(pszName, "SrcOvrLevel"))
571 : {
572 185 : const int nOldValue = m_nSrcOvrLevel;
573 185 : if (pszValue == nullptr || EQUAL(pszValue, "AUTO"))
574 1 : m_nSrcOvrLevel = -2;
575 184 : else if (STARTS_WITH_CI(pszValue, "AUTO-"))
576 1 : m_nSrcOvrLevel = -2 - atoi(pszValue + 5);
577 183 : else if (EQUAL(pszValue, "NONE"))
578 1 : m_nSrcOvrLevel = -1;
579 182 : else if (CPLGetValueType(pszValue) == CPL_VALUE_INTEGER)
580 182 : m_nSrcOvrLevel = atoi(pszValue);
581 185 : if (m_nSrcOvrLevel != nOldValue)
582 2 : SetNeedsFlush();
583 185 : return CE_None;
584 : }
585 1 : return VRTDataset::SetMetadataItem(pszName, pszValue, pszDomain);
586 : }
587 :
588 : /************************************************************************/
589 : /* VRTWarpedAddOptions() */
590 : /************************************************************************/
591 :
592 390 : static char **VRTWarpedAddOptions(char **papszWarpOptions)
593 : {
594 : /* Avoid errors when adding an alpha band, but source dataset has */
595 : /* no alpha band (#4571), and generally don't leave our buffer uninitialized
596 : */
597 390 : if (CSLFetchNameValue(papszWarpOptions, "INIT_DEST") == nullptr)
598 49 : papszWarpOptions = CSLSetNameValue(papszWarpOptions, "INIT_DEST", "0");
599 :
600 : /* For https://github.com/OSGeo/gdal/issues/1985 */
601 390 : if (CSLFetchNameValue(papszWarpOptions,
602 390 : "ERROR_OUT_IF_EMPTY_SOURCE_WINDOW") == nullptr)
603 : {
604 311 : papszWarpOptions = CSLSetNameValue(
605 : papszWarpOptions, "ERROR_OUT_IF_EMPTY_SOURCE_WINDOW", "FALSE");
606 : }
607 390 : return papszWarpOptions;
608 : }
609 :
610 : /************************************************************************/
611 : /* Initialize() */
612 : /* */
613 : /* Initialize a dataset from passed in warp options. */
614 : /************************************************************************/
615 :
616 243 : CPLErr VRTWarpedDataset::Initialize(void *psWO)
617 :
618 : {
619 243 : if (m_poWarper != nullptr)
620 0 : delete m_poWarper;
621 :
622 243 : m_poWarper = new GDALWarpOperation();
623 :
624 : GDALWarpOptions *psWO_Dup =
625 243 : GDALCloneWarpOptions(static_cast<GDALWarpOptions *>(psWO));
626 :
627 243 : psWO_Dup->papszWarpOptions =
628 243 : VRTWarpedAddOptions(psWO_Dup->papszWarpOptions);
629 :
630 243 : CPLErr eErr = m_poWarper->Initialize(psWO_Dup);
631 :
632 : // The act of initializing this warped dataset with this warp options
633 : // will result in our assuming ownership of a reference to the
634 : // hSrcDS.
635 :
636 243 : if (eErr == CE_None &&
637 242 : static_cast<GDALWarpOptions *>(psWO)->hSrcDS != nullptr)
638 : {
639 242 : GDALReferenceDataset(psWO_Dup->hSrcDS);
640 : }
641 :
642 243 : GDALDestroyWarpOptions(psWO_Dup);
643 :
644 243 : if (nBands > 1)
645 : {
646 60 : GDALDataset::SetMetadataItem(GDALMD_INTERLEAVE, "PIXEL",
647 : GDAL_MDD_IMAGE_STRUCTURE);
648 : }
649 :
650 243 : return eErr;
651 : }
652 :
653 : /************************************************************************/
654 : /* GDALWarpCoordRescaler */
655 : /************************************************************************/
656 :
657 1 : class GDALWarpCoordRescaler final : public OGRCoordinateTransformation
658 : {
659 : const double m_dfRatioX;
660 : const double m_dfRatioY;
661 :
662 : GDALWarpCoordRescaler &operator=(const GDALWarpCoordRescaler &) = delete;
663 : GDALWarpCoordRescaler(GDALWarpCoordRescaler &&) = delete;
664 : GDALWarpCoordRescaler &operator=(GDALWarpCoordRescaler &&) = delete;
665 :
666 : public:
667 1 : GDALWarpCoordRescaler(double dfRatioX, double dfRatioY)
668 1 : : m_dfRatioX(dfRatioX), m_dfRatioY(dfRatioY)
669 : {
670 1 : }
671 :
672 0 : GDALWarpCoordRescaler(const GDALWarpCoordRescaler &) = default;
673 :
674 : ~GDALWarpCoordRescaler() override;
675 :
676 0 : const OGRSpatialReference *GetSourceCS() const override
677 : {
678 0 : return nullptr;
679 : }
680 :
681 6 : const OGRSpatialReference *GetTargetCS() const override
682 : {
683 6 : return nullptr;
684 : }
685 :
686 3 : virtual int Transform(size_t nCount, double *x, double *y, double * /*z*/,
687 : double * /*t*/, int *pabSuccess) override
688 : {
689 17 : for (size_t i = 0; i < nCount; i++)
690 : {
691 14 : x[i] *= m_dfRatioX;
692 14 : y[i] *= m_dfRatioY;
693 14 : if (pabSuccess)
694 14 : pabSuccess[i] = TRUE;
695 : }
696 3 : return TRUE;
697 : }
698 :
699 0 : OGRCoordinateTransformation *Clone() const override
700 : {
701 0 : return new GDALWarpCoordRescaler(*this);
702 : }
703 :
704 0 : OGRCoordinateTransformation *GetInverse() const override
705 : {
706 0 : return nullptr;
707 : }
708 : };
709 :
710 : GDALWarpCoordRescaler::~GDALWarpCoordRescaler() = default;
711 :
712 : /************************************************************************/
713 : /* RescaleDstGeoTransform() */
714 : /************************************************************************/
715 :
716 23 : static void RescaleDstGeoTransform(double adfDstGeoTransform[6],
717 : int nRasterXSize, int nDstPixels,
718 : int nRasterYSize, int nDstLines)
719 : {
720 23 : adfDstGeoTransform[1] *= static_cast<double>(nRasterXSize) / nDstPixels;
721 23 : adfDstGeoTransform[2] *= static_cast<double>(nRasterXSize) / nDstPixels;
722 23 : adfDstGeoTransform[4] *= static_cast<double>(nRasterYSize) / nDstLines;
723 23 : adfDstGeoTransform[5] *= static_cast<double>(nRasterYSize) / nDstLines;
724 23 : }
725 :
726 : /************************************************************************/
727 : /* GetSrcOverviewLevel() */
728 : /************************************************************************/
729 :
730 81 : int VRTWarpedDataset::GetSrcOverviewLevel(int iOvr,
731 : bool &bThisLevelOnlyOut) const
732 : {
733 81 : bThisLevelOnlyOut = false;
734 81 : if (m_nSrcOvrLevel < -2)
735 : {
736 6 : if (iOvr + m_nSrcOvrLevel + 2 >= 0)
737 : {
738 3 : return iOvr + m_nSrcOvrLevel + 2;
739 : }
740 : }
741 75 : else if (m_nSrcOvrLevel == -2)
742 : {
743 69 : return iOvr;
744 : }
745 6 : else if (m_nSrcOvrLevel >= 0)
746 : {
747 0 : bThisLevelOnlyOut = true;
748 0 : return m_nSrcOvrLevel;
749 : }
750 9 : return -1;
751 : }
752 :
753 : /************************************************************************/
754 : /* GetOverviewSize() */
755 : /************************************************************************/
756 :
757 75 : bool VRTWarpedDataset::GetOverviewSize(GDALDataset *poSrcDS, int iOvr,
758 : int iSrcOvr, int &nOvrXSize,
759 : int &nOvrYSize, double &dfSrcRatioX,
760 : double &dfSrcRatioY) const
761 : {
762 : auto poSrcOvrBand = iSrcOvr >= 0
763 75 : ? poSrcDS->GetRasterBand(1)->GetOverview(iSrcOvr)
764 3 : : poSrcDS->GetRasterBand(1);
765 75 : if (!poSrcOvrBand)
766 : {
767 0 : return false;
768 : }
769 75 : dfSrcRatioX = static_cast<double>(poSrcDS->GetRasterXSize()) /
770 75 : poSrcOvrBand->GetXSize();
771 75 : dfSrcRatioY = static_cast<double>(poSrcDS->GetRasterYSize()) /
772 75 : poSrcOvrBand->GetYSize();
773 : const double dfTargetRatio =
774 75 : static_cast<double>(poSrcDS->GetRasterXSize()) /
775 75 : poSrcDS->GetRasterBand(1)->GetOverview(iOvr)->GetXSize();
776 :
777 75 : nOvrXSize = static_cast<int>(nRasterXSize / dfTargetRatio + 0.5);
778 75 : nOvrYSize = static_cast<int>(nRasterYSize / dfTargetRatio + 0.5);
779 75 : return nOvrXSize >= 1 && nOvrYSize >= 1;
780 : }
781 :
782 : /************************************************************************/
783 : /* CreateImplicitOverview() */
784 : /************************************************************************/
785 :
786 21 : VRTWarpedDataset *VRTWarpedDataset::CreateImplicitOverview(int iOvr) const
787 : {
788 21 : if (!m_poWarper)
789 0 : return nullptr;
790 21 : const GDALWarpOptions *psWO = m_poWarper->GetOptions();
791 21 : if (!psWO->hSrcDS || GDALGetRasterCount(psWO->hSrcDS) == 0)
792 0 : return nullptr;
793 21 : GDALDataset *poSrcDS = GDALDataset::FromHandle(psWO->hSrcDS);
794 21 : GDALDataset *poSrcOvrDS = poSrcDS;
795 21 : bool bThisLevelOnly = false;
796 21 : const int iSrcOvr = GetSrcOverviewLevel(iOvr, bThisLevelOnly);
797 21 : if (iSrcOvr >= 0)
798 : {
799 : poSrcOvrDS =
800 18 : GDALCreateOverviewDataset(poSrcDS, iSrcOvr, bThisLevelOnly);
801 : }
802 21 : if (poSrcOvrDS == nullptr)
803 0 : return nullptr;
804 21 : if (poSrcOvrDS == poSrcDS)
805 3 : poSrcOvrDS->Reference();
806 :
807 21 : int nDstPixels = 0;
808 21 : int nDstLines = 0;
809 21 : double dfSrcRatioX = 0;
810 21 : double dfSrcRatioY = 0;
811 : // Figure out the desired output bounds and resolution.
812 21 : if (!GetOverviewSize(poSrcDS, iOvr, iSrcOvr, nDstPixels, nDstLines,
813 : dfSrcRatioX, dfSrcRatioY))
814 : {
815 0 : poSrcOvrDS->ReleaseRef();
816 0 : return nullptr;
817 : }
818 :
819 : /* --------------------------------------------------------------------
820 : */
821 : /* Create transformer and warping options. */
822 : /* --------------------------------------------------------------------
823 : */
824 42 : void *pTransformerArg = GDALCreateSimilarTransformer(
825 21 : psWO->pTransformerArg, dfSrcRatioX, dfSrcRatioY);
826 21 : if (pTransformerArg == nullptr)
827 : {
828 0 : poSrcOvrDS->ReleaseRef();
829 0 : return nullptr;
830 : }
831 :
832 21 : GDALWarpOptions *psWOOvr = GDALCloneWarpOptions(psWO);
833 21 : psWOOvr->hSrcDS = poSrcOvrDS;
834 21 : psWOOvr->pfnTransformer = psWO->pfnTransformer;
835 21 : psWOOvr->pTransformerArg = pTransformerArg;
836 :
837 : /* --------------------------------------------------------------------
838 : */
839 : /* We need to rescale the potential CUTLINE */
840 : /* --------------------------------------------------------------------
841 : */
842 21 : if (psWOOvr->hCutline)
843 : {
844 2 : GDALWarpCoordRescaler oRescaler(1.0 / dfSrcRatioX, 1.0 / dfSrcRatioY);
845 1 : static_cast<OGRGeometry *>(psWOOvr->hCutline)->transform(&oRescaler);
846 : }
847 :
848 : /* --------------------------------------------------------------------
849 : */
850 : /* Rescale the output geotransform on the transformer. */
851 : /* --------------------------------------------------------------------
852 : */
853 21 : double adfDstGeoTransform[6] = {0.0};
854 21 : GDALGetTransformerDstGeoTransform(psWOOvr->pTransformerArg,
855 : adfDstGeoTransform);
856 21 : RescaleDstGeoTransform(adfDstGeoTransform, nRasterXSize, nDstPixels,
857 21 : nRasterYSize, nDstLines);
858 21 : GDALSetTransformerDstGeoTransform(psWOOvr->pTransformerArg,
859 : adfDstGeoTransform);
860 :
861 : /* --------------------------------------------------------------------
862 : */
863 : /* Create the VRT file. */
864 : /* --------------------------------------------------------------------
865 : */
866 21 : GDALDatasetH hDstDS = GDALCreateWarpedVRT(poSrcOvrDS, nDstPixels, nDstLines,
867 : adfDstGeoTransform, psWOOvr);
868 :
869 21 : poSrcOvrDS->ReleaseRef();
870 :
871 21 : GDALDestroyWarpOptions(psWOOvr);
872 :
873 21 : if (hDstDS == nullptr)
874 : {
875 0 : GDALDestroyTransformer(pTransformerArg);
876 0 : return nullptr;
877 : }
878 :
879 21 : auto poOvrDS = static_cast<VRTWarpedDataset *>(hDstDS);
880 21 : poOvrDS->m_bIsOverview = true;
881 21 : return poOvrDS;
882 : }
883 :
884 : /************************************************************************/
885 : /* GetOverviewCount() */
886 : /************************************************************************/
887 :
888 112 : int VRTWarpedDataset::GetOverviewCount() const
889 : {
890 112 : if (m_poWarper)
891 : {
892 112 : const GDALWarpOptions *psWO = m_poWarper->GetOptions();
893 112 : if (!m_bIsOverview && psWO->hSrcDS && GDALGetRasterCount(psWO->hSrcDS))
894 : {
895 112 : GDALDataset *poSrcDS = GDALDataset::FromHandle(psWO->hSrcDS);
896 : int nSrcOverviewCount =
897 112 : poSrcDS->GetRasterBand(1)->GetOverviewCount();
898 112 : int nCount = 0;
899 172 : for (int i = 0; i < nSrcOverviewCount; ++i)
900 : {
901 60 : bool bThisLevelOnly = false;
902 60 : const int iSrcOvr = GetSrcOverviewLevel(i, bThisLevelOnly);
903 60 : if (iSrcOvr >= 0)
904 : {
905 54 : int nDstPixels = 0;
906 54 : int nDstLines = 0;
907 54 : double dfSrcRatioX = 0;
908 54 : double dfSrcRatioY = 0;
909 54 : if (!GetOverviewSize(poSrcDS, i, iSrcOvr, nDstPixels,
910 : nDstLines, dfSrcRatioX, dfSrcRatioY))
911 : {
912 0 : break;
913 : }
914 : }
915 60 : ++nCount;
916 : }
917 112 : return nCount;
918 : }
919 : }
920 0 : return 0;
921 : }
922 :
923 : /************************************************************************/
924 : /* CreateImplicitOverviews() */
925 : /* */
926 : /* For each overview of the source dataset, create an overview */
927 : /* in the warped VRT dataset. */
928 : /************************************************************************/
929 :
930 8 : void VRTWarpedDataset::CreateImplicitOverviews()
931 : {
932 8 : if (m_bIsOverview)
933 0 : return;
934 8 : const int nOvrCount = GetOverviewCount();
935 8 : if (m_apoOverviews.empty())
936 4 : m_apoOverviews.resize(nOvrCount);
937 18 : for (int iOvr = 0; iOvr < nOvrCount; iOvr++)
938 : {
939 10 : if (!m_apoOverviews[iOvr])
940 : {
941 2 : m_apoOverviews[iOvr] = CreateImplicitOverview(iOvr);
942 : }
943 : }
944 : }
945 :
946 : /************************************************************************/
947 : /* GetFileList() */
948 : /************************************************************************/
949 :
950 20 : char **VRTWarpedDataset::GetFileList()
951 : {
952 20 : char **papszFileList = GDALDataset::GetFileList();
953 :
954 20 : if (m_poWarper != nullptr)
955 : {
956 20 : const GDALWarpOptions *psWO = m_poWarper->GetOptions();
957 20 : const char *pszFilename = nullptr;
958 :
959 40 : if (psWO->hSrcDS != nullptr &&
960 : (pszFilename =
961 20 : static_cast<GDALDataset *>(psWO->hSrcDS)->GetDescription()) !=
962 : nullptr)
963 : {
964 : VSIStatBufL sStat;
965 20 : if (VSIStatL(pszFilename, &sStat) == 0)
966 : {
967 20 : papszFileList = CSLAddString(papszFileList, pszFilename);
968 : }
969 : }
970 : }
971 :
972 20 : return papszFileList;
973 : }
974 :
975 : /************************************************************************/
976 : /* ==================================================================== */
977 : /* VRTWarpedOverviewTransformer */
978 : /* ==================================================================== */
979 : /************************************************************************/
980 :
981 : typedef struct
982 : {
983 : GDALTransformerInfo sTI;
984 :
985 : GDALTransformerFunc pfnBaseTransformer;
986 : void *pBaseTransformerArg;
987 : bool bOwnSubtransformer;
988 :
989 : double dfXOverviewFactor;
990 : double dfYOverviewFactor;
991 : } VWOTInfo;
992 :
993 : static void *VRTCreateWarpedOverviewTransformer(
994 : GDALTransformerFunc pfnBaseTransformer, void *pBaseTransformArg,
995 : double dfXOverviewFactor, double dfYOverviewFactor);
996 : static void VRTDestroyWarpedOverviewTransformer(void *pTransformArg);
997 :
998 : static int VRTWarpedOverviewTransform(void *pTransformArg, int bDstToSrc,
999 : int nPointCount, double *padfX,
1000 : double *padfY, double *padfZ,
1001 : int *panSuccess);
1002 :
1003 : #if 0 // TODO: Why?
1004 : /************************************************************************/
1005 : /* VRTSerializeWarpedOverviewTransformer() */
1006 : /************************************************************************/
1007 :
1008 : static CPLXMLNode *
1009 : VRTSerializeWarpedOverviewTransformer( void *pTransformArg )
1010 :
1011 : {
1012 : VWOTInfo *psInfo = static_cast<VWOTInfo *>( pTransformArg );
1013 :
1014 : CPLXMLNode *psTree
1015 : = CPLCreateXMLNode( NULL, CXT_Element, "WarpedOverviewTransformer" );
1016 :
1017 : CPLCreateXMLElementAndValue(
1018 : psTree, "XFactor",
1019 : CPLString().Printf("%g",psInfo->dfXOverviewFactor) );
1020 : CPLCreateXMLElementAndValue(
1021 : psTree, "YFactor",
1022 : CPLString().Printf("%g",psInfo->dfYOverviewFactor) );
1023 :
1024 : /* -------------------------------------------------------------------- */
1025 : /* Capture underlying transformer. */
1026 : /* -------------------------------------------------------------------- */
1027 : CPLXMLNode *psTransformerContainer
1028 : = CPLCreateXMLNode( psTree, CXT_Element, "BaseTransformer" );
1029 :
1030 : CPLXMLNode *psTransformer
1031 : = GDALSerializeTransformer( psInfo->pfnBaseTransformer,
1032 : psInfo->pBaseTransformerArg );
1033 : if( psTransformer != NULL )
1034 : CPLAddXMLChild( psTransformerContainer, psTransformer );
1035 :
1036 : return psTree;
1037 : }
1038 :
1039 : /************************************************************************/
1040 : /* VRTWarpedOverviewTransformerOwnsSubtransformer() */
1041 : /************************************************************************/
1042 :
1043 : static void VRTWarpedOverviewTransformerOwnsSubtransformer( void *pTransformArg,
1044 : bool bOwnFlag )
1045 : {
1046 : VWOTInfo *psInfo = static_cast<VWOTInfo *>( pTransformArg );
1047 :
1048 : psInfo->bOwnSubtransformer = bOwnFlag;
1049 : }
1050 :
1051 : /************************************************************************/
1052 : /* VRTDeserializeWarpedOverviewTransformer() */
1053 : /************************************************************************/
1054 :
1055 : void* VRTDeserializeWarpedOverviewTransformer( CPLXMLNode *psTree )
1056 :
1057 : {
1058 : const double dfXOverviewFactor =
1059 : CPLAtof(CPLGetXMLValue( psTree, "XFactor", "1" ));
1060 : const double dfYOverviewFactor =
1061 : CPLAtof(CPLGetXMLValue( psTree, "YFactor", "1" ));
1062 : GDALTransformerFunc pfnBaseTransform = NULL;
1063 : void *pBaseTransformerArg = NULL;
1064 :
1065 : CPLXMLNode *psContainer = CPLGetXMLNode( psTree, "BaseTransformer" );
1066 :
1067 : if( psContainer != NULL && psContainer->psChild != NULL )
1068 : {
1069 : GDALDeserializeTransformer( psContainer->psChild,
1070 : &pfnBaseTransform,
1071 : &pBaseTransformerArg );
1072 : }
1073 :
1074 : if( pfnBaseTransform == NULL )
1075 : {
1076 : CPLError( CE_Failure, CPLE_AppDefined,
1077 : "Cannot get base transform for scaled coord transformer." );
1078 : return NULL;
1079 : }
1080 : else
1081 : {
1082 : void *pApproxCBData =
1083 : VRTCreateWarpedOverviewTransformer( pfnBaseTransform,
1084 : pBaseTransformerArg,
1085 : dfXOverviewFactor,
1086 : dfYOverviewFactor );
1087 : VRTWarpedOverviewTransformerOwnsSubtransformer( pApproxCBData, true );
1088 :
1089 : return pApproxCBData;
1090 : }
1091 : }
1092 : #endif // TODO: Why disabled?
1093 :
1094 : /************************************************************************/
1095 : /* VRTCreateWarpedOverviewTransformer() */
1096 : /************************************************************************/
1097 :
1098 4 : static void *VRTCreateWarpedOverviewTransformer(
1099 : GDALTransformerFunc pfnBaseTransformer, void *pBaseTransformerArg,
1100 : double dfXOverviewFactor, double dfYOverviewFactor)
1101 :
1102 : {
1103 4 : if (pfnBaseTransformer == nullptr)
1104 0 : return nullptr;
1105 :
1106 4 : VWOTInfo *psSCTInfo = static_cast<VWOTInfo *>(CPLMalloc(sizeof(VWOTInfo)));
1107 4 : psSCTInfo->pfnBaseTransformer = pfnBaseTransformer;
1108 4 : psSCTInfo->pBaseTransformerArg = pBaseTransformerArg;
1109 4 : psSCTInfo->dfXOverviewFactor = dfXOverviewFactor;
1110 4 : psSCTInfo->dfYOverviewFactor = dfYOverviewFactor;
1111 4 : psSCTInfo->bOwnSubtransformer = false;
1112 :
1113 4 : memcpy(psSCTInfo->sTI.abySignature, GDAL_GTI2_SIGNATURE,
1114 : strlen(GDAL_GTI2_SIGNATURE));
1115 4 : psSCTInfo->sTI.pszClassName = "VRTWarpedOverviewTransformer";
1116 4 : psSCTInfo->sTI.pfnTransform = VRTWarpedOverviewTransform;
1117 4 : psSCTInfo->sTI.pfnCleanup = VRTDestroyWarpedOverviewTransformer;
1118 : #if 0
1119 : psSCTInfo->sTI.pfnSerialize = VRTSerializeWarpedOverviewTransformer;
1120 : #endif
1121 4 : return psSCTInfo;
1122 : }
1123 :
1124 : /************************************************************************/
1125 : /* VRTDestroyWarpedOverviewTransformer() */
1126 : /************************************************************************/
1127 :
1128 4 : static void VRTDestroyWarpedOverviewTransformer(void *pTransformArg)
1129 : {
1130 4 : VWOTInfo *psInfo = static_cast<VWOTInfo *>(pTransformArg);
1131 :
1132 4 : if (psInfo->bOwnSubtransformer)
1133 0 : GDALDestroyTransformer(psInfo->pBaseTransformerArg);
1134 :
1135 4 : CPLFree(psInfo);
1136 4 : }
1137 :
1138 : /************************************************************************/
1139 : /* VRTWarpedOverviewTransform() */
1140 : /************************************************************************/
1141 :
1142 262 : static int VRTWarpedOverviewTransform(void *pTransformArg, int bDstToSrc,
1143 : int nPointCount, double *padfX,
1144 : double *padfY, double *padfZ,
1145 : int *panSuccess)
1146 :
1147 : {
1148 262 : VWOTInfo *psInfo = static_cast<VWOTInfo *>(pTransformArg);
1149 :
1150 262 : if (bDstToSrc)
1151 : {
1152 63871 : for (int i = 0; i < nPointCount; i++)
1153 : {
1154 63609 : padfX[i] *= psInfo->dfXOverviewFactor;
1155 63609 : padfY[i] *= psInfo->dfYOverviewFactor;
1156 : }
1157 : }
1158 :
1159 262 : const int bSuccess = psInfo->pfnBaseTransformer(
1160 : psInfo->pBaseTransformerArg, bDstToSrc, nPointCount, padfX, padfY,
1161 : padfZ, panSuccess);
1162 :
1163 262 : if (!bDstToSrc)
1164 : {
1165 0 : for (int i = 0; i < nPointCount; i++)
1166 : {
1167 0 : padfX[i] /= psInfo->dfXOverviewFactor;
1168 0 : padfY[i] /= psInfo->dfYOverviewFactor;
1169 : }
1170 : }
1171 :
1172 262 : return bSuccess;
1173 : }
1174 :
1175 : /************************************************************************/
1176 : /* BuildOverviews() */
1177 : /* */
1178 : /* For overviews, we actually just build a whole new dataset */
1179 : /* with an extra layer of transformation on the warper used to */
1180 : /* accomplish downsampling by the desired factor. */
1181 : /************************************************************************/
1182 :
1183 6 : CPLErr VRTWarpedDataset::IBuildOverviews(
1184 : const char * /* pszResampling */, int nOverviews,
1185 : const int *panOverviewList, int /* nListBands */,
1186 : const int * /* panBandList */, GDALProgressFunc pfnProgress,
1187 : void *pProgressData, CSLConstList /*papszOptions*/)
1188 : {
1189 6 : if (m_poWarper == nullptr || m_bIsOverview)
1190 0 : return CE_Failure;
1191 :
1192 : /* -------------------------------------------------------------------- */
1193 : /* Initial progress result. */
1194 : /* -------------------------------------------------------------------- */
1195 6 : if (!pfnProgress(0.0, nullptr, pProgressData))
1196 : {
1197 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
1198 0 : return CE_Failure;
1199 : }
1200 :
1201 6 : CreateImplicitOverviews();
1202 :
1203 : /* -------------------------------------------------------------------- */
1204 : /* Establish which of the overview levels we already have, and */
1205 : /* which are new. */
1206 : /* -------------------------------------------------------------------- */
1207 6 : int nNewOverviews = 0;
1208 : int *panNewOverviewList =
1209 6 : static_cast<int *>(CPLCalloc(sizeof(int), nOverviews));
1210 6 : std::vector<bool> abFoundOverviewFactor(nOverviews);
1211 14 : for (int i = 0; i < nOverviews; i++)
1212 : {
1213 20 : for (GDALDataset *const poOverview : m_apoOverviews)
1214 : {
1215 12 : if (poOverview)
1216 : {
1217 12 : const int nOvFactor = GDALComputeOvFactor(
1218 : poOverview->GetRasterXSize(), GetRasterXSize(),
1219 : poOverview->GetRasterYSize(), GetRasterYSize());
1220 :
1221 20 : if (nOvFactor == panOverviewList[i] ||
1222 8 : nOvFactor == GDALOvLevelAdjust2(panOverviewList[i],
1223 : GetRasterXSize(),
1224 : GetRasterYSize()))
1225 4 : abFoundOverviewFactor[i] = true;
1226 : }
1227 : }
1228 :
1229 8 : if (!abFoundOverviewFactor[i])
1230 4 : panNewOverviewList[nNewOverviews++] = panOverviewList[i];
1231 : }
1232 :
1233 : /* -------------------------------------------------------------------- */
1234 : /* Create each missing overview (we don't need to do anything */
1235 : /* to update existing overviews). */
1236 : /* -------------------------------------------------------------------- */
1237 6 : CPLErr eErr = CE_None;
1238 10 : for (int i = 0; i < nNewOverviews; i++)
1239 : {
1240 : /* --------------------------------------------------------------------
1241 : */
1242 : /* What size should this overview be. */
1243 : /* --------------------------------------------------------------------
1244 : */
1245 4 : const int nOXSize = (GetRasterXSize() + panNewOverviewList[i] - 1) /
1246 4 : panNewOverviewList[i];
1247 :
1248 4 : const int nOYSize = (GetRasterYSize() + panNewOverviewList[i] - 1) /
1249 4 : panNewOverviewList[i];
1250 :
1251 : /* --------------------------------------------------------------------
1252 : */
1253 : /* Find the most appropriate base dataset onto which to build the
1254 : */
1255 : /* new one. The preference will be an overview dataset with a
1256 : * ratio*/
1257 : /* greater than ours, and which is not using */
1258 : /* VRTWarpedOverviewTransform, since those ones are slow. The
1259 : * other*/
1260 : /* ones are based on overviews of the source dataset. */
1261 : /* --------------------------------------------------------------------
1262 : */
1263 4 : VRTWarpedDataset *poBaseDataset = this;
1264 8 : for (auto *poOverview : m_apoOverviews)
1265 : {
1266 4 : if (poOverview && poOverview->GetRasterXSize() > nOXSize &&
1267 4 : poOverview->m_poWarper->GetOptions()->pfnTransformer !=
1268 8 : VRTWarpedOverviewTransform &&
1269 4 : poOverview->GetRasterXSize() < poBaseDataset->GetRasterXSize())
1270 : {
1271 4 : poBaseDataset = poOverview;
1272 : }
1273 : }
1274 :
1275 : /* --------------------------------------------------------------------
1276 : */
1277 : /* Create the overview dataset. */
1278 : /* --------------------------------------------------------------------
1279 : */
1280 4 : VRTWarpedDataset *poOverviewDS = new VRTWarpedDataset(nOXSize, nOYSize);
1281 :
1282 8 : for (int iBand = 0; iBand < GetRasterCount(); iBand++)
1283 : {
1284 4 : GDALRasterBand *const poOldBand = GetRasterBand(iBand + 1);
1285 : VRTWarpedRasterBand *const poNewBand = new VRTWarpedRasterBand(
1286 4 : poOverviewDS, iBand + 1, poOldBand->GetRasterDataType());
1287 :
1288 4 : poNewBand->CopyCommonInfoFrom(poOldBand);
1289 4 : poOverviewDS->SetBand(iBand + 1, poNewBand);
1290 : }
1291 :
1292 : /* --------------------------------------------------------------------
1293 : */
1294 : /* Prepare update transformation information that will apply */
1295 : /* the overview decimation. */
1296 : /* --------------------------------------------------------------------
1297 : */
1298 : GDALWarpOptions *psWO = const_cast<GDALWarpOptions *>(
1299 4 : poBaseDataset->m_poWarper->GetOptions());
1300 :
1301 : /* --------------------------------------------------------------------
1302 : */
1303 : /* Initialize the new dataset with adjusted warp options, and */
1304 : /* then restore to original condition. */
1305 : /* --------------------------------------------------------------------
1306 : */
1307 4 : GDALTransformerFunc pfnTransformerBase = psWO->pfnTransformer;
1308 4 : void *pTransformerBaseArg = psWO->pTransformerArg;
1309 :
1310 4 : psWO->pfnTransformer = VRTWarpedOverviewTransform;
1311 12 : psWO->pTransformerArg = VRTCreateWarpedOverviewTransformer(
1312 : pfnTransformerBase, pTransformerBaseArg,
1313 4 : poBaseDataset->GetRasterXSize() / static_cast<double>(nOXSize),
1314 4 : poBaseDataset->GetRasterYSize() / static_cast<double>(nOYSize));
1315 :
1316 4 : eErr = poOverviewDS->Initialize(psWO);
1317 :
1318 4 : psWO->pfnTransformer = pfnTransformerBase;
1319 4 : psWO->pTransformerArg = pTransformerBaseArg;
1320 :
1321 4 : if (eErr != CE_None)
1322 : {
1323 0 : delete poOverviewDS;
1324 0 : break;
1325 : }
1326 :
1327 4 : m_apoOverviews.push_back(poOverviewDS);
1328 : }
1329 :
1330 6 : CPLFree(panNewOverviewList);
1331 :
1332 : /* -------------------------------------------------------------------- */
1333 : /* Progress finished. */
1334 : /* -------------------------------------------------------------------- */
1335 6 : pfnProgress(1.0, nullptr, pProgressData);
1336 :
1337 6 : SetNeedsFlush();
1338 :
1339 6 : return eErr;
1340 : }
1341 :
1342 : /*! @endcond */
1343 :
1344 : /************************************************************************/
1345 : /* GDALInitializeWarpedVRT() */
1346 : /************************************************************************/
1347 :
1348 : /**
1349 : * Set warp info on virtual warped dataset.
1350 : *
1351 : * Initializes all the warping information for a virtual warped dataset.
1352 : *
1353 : * This method is the same as the C++ method VRTWarpedDataset::Initialize().
1354 : *
1355 : * @param hDS dataset previously created with the VRT driver, and a
1356 : * SUBCLASS of "VRTWarpedDataset".
1357 : *
1358 : * @param psWO the warp options to apply. Note that ownership of the
1359 : * transformation information is taken over by the function though everything
1360 : * else remains the property of the caller.
1361 : *
1362 : * @return CE_None on success or CE_Failure if an error occurs.
1363 : */
1364 :
1365 182 : CPLErr CPL_STDCALL GDALInitializeWarpedVRT(GDALDatasetH hDS,
1366 : GDALWarpOptions *psWO)
1367 :
1368 : {
1369 182 : VALIDATE_POINTER1(hDS, "GDALInitializeWarpedVRT", CE_Failure);
1370 :
1371 182 : return static_cast<VRTWarpedDataset *>(GDALDataset::FromHandle(hDS))
1372 182 : ->Initialize(psWO);
1373 : }
1374 :
1375 : /*! @cond Doxygen_Suppress */
1376 :
1377 : /************************************************************************/
1378 : /* XMLInit() */
1379 : /************************************************************************/
1380 :
1381 149 : CPLErr VRTWarpedDataset::XMLInit(const CPLXMLNode *psTree,
1382 : const char *pszVRTPathIn)
1383 :
1384 : {
1385 :
1386 : /* -------------------------------------------------------------------- */
1387 : /* Initialize blocksize before calling sub-init so that the */
1388 : /* band initializers can get it from the dataset object when */
1389 : /* they are created. */
1390 : /* -------------------------------------------------------------------- */
1391 149 : m_nBlockXSize = atoi(CPLGetXMLValue(psTree, "BlockXSize", "512"));
1392 149 : m_nBlockYSize = atoi(CPLGetXMLValue(psTree, "BlockYSize", "128"));
1393 :
1394 : /* -------------------------------------------------------------------- */
1395 : /* Initialize all the general VRT stuff. This will even */
1396 : /* create the VRTWarpedRasterBands and initialize them. */
1397 : /* -------------------------------------------------------------------- */
1398 : {
1399 149 : const CPLErr eErr = VRTDataset::XMLInit(psTree, pszVRTPathIn);
1400 :
1401 149 : if (eErr != CE_None)
1402 0 : return eErr;
1403 : }
1404 :
1405 : // Check that band block sizes didn't change the dataset block size.
1406 405 : for (int i = 1; i <= nBands; i++)
1407 : {
1408 258 : int nBlockXSize = 0;
1409 258 : int nBlockYSize = 0;
1410 258 : GetRasterBand(i)->GetBlockSize(&nBlockXSize, &nBlockYSize);
1411 258 : if (nBlockXSize != m_nBlockXSize || nBlockYSize != m_nBlockYSize)
1412 : {
1413 2 : CPLError(CE_Failure, CPLE_AppDefined,
1414 : "Block size specified on band %d not consistent with "
1415 : "dataset block size",
1416 : i);
1417 2 : return CE_Failure;
1418 : }
1419 : }
1420 :
1421 147 : if (nBands > 1)
1422 : {
1423 42 : GDALDataset::SetMetadataItem(GDALMD_INTERLEAVE, "PIXEL",
1424 : GDAL_MDD_IMAGE_STRUCTURE);
1425 : }
1426 :
1427 : /* -------------------------------------------------------------------- */
1428 : /* Find the GDALWarpOptions XML tree. */
1429 : /* -------------------------------------------------------------------- */
1430 : const CPLXMLNode *const psOptionsTree =
1431 147 : CPLGetXMLNode(psTree, "GDALWarpOptions");
1432 147 : if (psOptionsTree == nullptr)
1433 : {
1434 0 : CPLError(CE_Failure, CPLE_AppDefined,
1435 : "Count not find required GDALWarpOptions in XML.");
1436 0 : return CE_Failure;
1437 : }
1438 :
1439 : /* -------------------------------------------------------------------- */
1440 : /* Adjust the SourceDataset in the warp options to take into */
1441 : /* account that it is relative to the VRT if appropriate. */
1442 : /* -------------------------------------------------------------------- */
1443 147 : const bool bRelativeToVRT = CPL_TO_BOOL(atoi(
1444 : CPLGetXMLValue(psOptionsTree, "SourceDataset.relativeToVRT", "0")));
1445 :
1446 : const char *pszRelativePath =
1447 147 : CPLGetXMLValue(psOptionsTree, "SourceDataset", "");
1448 147 : char *pszAbsolutePath = nullptr;
1449 :
1450 147 : if (bRelativeToVRT)
1451 114 : pszAbsolutePath = CPLStrdup(
1452 228 : CPLProjectRelativeFilenameSafe(pszVRTPathIn, pszRelativePath)
1453 : .c_str());
1454 : else
1455 33 : pszAbsolutePath = CPLStrdup(pszRelativePath);
1456 :
1457 147 : CPLXMLNode *psOptionsTreeCloned = CPLCloneXMLTree(psOptionsTree);
1458 147 : CPLSetXMLValue(psOptionsTreeCloned, "SourceDataset", pszAbsolutePath);
1459 147 : CPLFree(pszAbsolutePath);
1460 :
1461 : /* -------------------------------------------------------------------- */
1462 : /* And instantiate the warp options, and corresponding warp */
1463 : /* operation. */
1464 : /* -------------------------------------------------------------------- */
1465 147 : GDALWarpOptions *psWO = GDALDeserializeWarpOptions(psOptionsTreeCloned);
1466 147 : CPLDestroyXMLNode(psOptionsTreeCloned);
1467 147 : if (psWO == nullptr)
1468 0 : return CE_Failure;
1469 :
1470 147 : psWO->papszWarpOptions = VRTWarpedAddOptions(psWO->papszWarpOptions);
1471 :
1472 147 : eAccess = GA_Update;
1473 :
1474 147 : if (psWO->hDstDS != nullptr)
1475 : {
1476 0 : GDALClose(psWO->hDstDS);
1477 0 : psWO->hDstDS = nullptr;
1478 : }
1479 :
1480 147 : psWO->hDstDS = this;
1481 :
1482 : /* -------------------------------------------------------------------- */
1483 : /* Deserialize vertical shift grids. */
1484 : /* -------------------------------------------------------------------- */
1485 147 : CPLXMLNode *psIter = psTree->psChild;
1486 1652 : for (; psWO->hSrcDS != nullptr && psIter != nullptr;
1487 1505 : psIter = psIter->psNext)
1488 : {
1489 1505 : if (psIter->eType != CXT_Element ||
1490 1064 : !EQUAL(psIter->pszValue, "VerticalShiftGrids"))
1491 : {
1492 1505 : continue;
1493 : }
1494 :
1495 0 : CPLError(CE_Warning, CPLE_AppDefined,
1496 : "The VerticalShiftGrids in a warped VRT is now deprecated, "
1497 : "and will no longer be handled in GDAL 4.0");
1498 :
1499 0 : const char *pszVGrids = CPLGetXMLValue(psIter, "Grids", nullptr);
1500 0 : if (pszVGrids)
1501 : {
1502 : int bInverse =
1503 0 : CSLTestBoolean(CPLGetXMLValue(psIter, "Inverse", "FALSE"));
1504 : double dfToMeterSrc =
1505 0 : CPLAtof(CPLGetXMLValue(psIter, "ToMeterSrc", "1.0"));
1506 : double dfToMeterDest =
1507 0 : CPLAtof(CPLGetXMLValue(psIter, "ToMeterDest", "1.0"));
1508 0 : char **papszOptions = nullptr;
1509 0 : CPLXMLNode *psIter2 = psIter->psChild;
1510 0 : for (; psIter2 != nullptr; psIter2 = psIter2->psNext)
1511 : {
1512 0 : if (psIter2->eType != CXT_Element ||
1513 0 : !EQUAL(psIter2->pszValue, "Option"))
1514 : {
1515 0 : continue;
1516 : }
1517 0 : const char *pszName = CPLGetXMLValue(psIter2, "name", nullptr);
1518 : const char *pszValue =
1519 0 : CPLGetXMLValue(psIter2, nullptr, nullptr);
1520 0 : if (pszName && pszValue)
1521 : {
1522 : papszOptions =
1523 0 : CSLSetNameValue(papszOptions, pszName, pszValue);
1524 : }
1525 : }
1526 :
1527 0 : int bError = FALSE;
1528 : GDALDatasetH hGridDataset =
1529 0 : GDALOpenVerticalShiftGrid(pszVGrids, &bError);
1530 0 : if (bError && hGridDataset == nullptr)
1531 : {
1532 0 : CPLError(CE_Warning, CPLE_AppDefined,
1533 : "Cannot open %s. Source dataset will no "
1534 : "be vertically adjusted regarding "
1535 : "vertical datum",
1536 : pszVGrids);
1537 : }
1538 0 : else if (hGridDataset != nullptr)
1539 : {
1540 : // Transform from source vertical datum to WGS84
1541 0 : GDALDatasetH hTmpDS = GDALApplyVerticalShiftGrid(
1542 : psWO->hSrcDS, hGridDataset, bInverse, dfToMeterSrc,
1543 : dfToMeterDest, papszOptions);
1544 0 : GDALReleaseDataset(hGridDataset);
1545 0 : if (hTmpDS == nullptr)
1546 : {
1547 0 : CPLError(CE_Warning, CPLE_AppDefined,
1548 : "Source dataset will no "
1549 : "be vertically adjusted regarding "
1550 : "vertical datum %s",
1551 : pszVGrids);
1552 : }
1553 : else
1554 : {
1555 0 : CPLDebug("GDALWARP",
1556 : "Adjusting source dataset "
1557 : "with vertical datum using %s",
1558 : pszVGrids);
1559 0 : GDALReleaseDataset(psWO->hSrcDS);
1560 0 : psWO->hSrcDS = hTmpDS;
1561 : }
1562 : }
1563 :
1564 0 : CSLDestroy(papszOptions);
1565 : }
1566 : }
1567 :
1568 : /* -------------------------------------------------------------------- */
1569 : /* Instantiate the warp operation. */
1570 : /* -------------------------------------------------------------------- */
1571 147 : m_poWarper = new GDALWarpOperation();
1572 :
1573 147 : const CPLErr eErr = m_poWarper->Initialize(psWO);
1574 147 : if (eErr != CE_None)
1575 : {
1576 : /* --------------------------------------------------------------------
1577 : */
1578 : /* We are responsible for cleaning up the transformer ourselves. */
1579 : /* --------------------------------------------------------------------
1580 : */
1581 0 : if (psWO->pTransformerArg != nullptr)
1582 : {
1583 0 : GDALDestroyTransformer(psWO->pTransformerArg);
1584 0 : psWO->pTransformerArg = nullptr;
1585 : }
1586 :
1587 0 : if (psWO->hSrcDS != nullptr)
1588 : {
1589 0 : GDALClose(psWO->hSrcDS);
1590 0 : psWO->hSrcDS = nullptr;
1591 : }
1592 : }
1593 :
1594 147 : GDALDestroyWarpOptions(psWO);
1595 147 : if (eErr != CE_None)
1596 : {
1597 0 : delete m_poWarper;
1598 0 : m_poWarper = nullptr;
1599 : }
1600 :
1601 : /* -------------------------------------------------------------------- */
1602 : /* Deserialize SrcOvrLevel */
1603 : /* -------------------------------------------------------------------- */
1604 147 : const char *pszSrcOvrLevel = CPLGetXMLValue(psTree, "SrcOvrLevel", nullptr);
1605 147 : if (pszSrcOvrLevel != nullptr)
1606 : {
1607 0 : SetMetadataItem("SrcOvrLevel", pszSrcOvrLevel);
1608 : }
1609 :
1610 : /* -------------------------------------------------------------------- */
1611 : /* Generate overviews, if appropriate. */
1612 : /* -------------------------------------------------------------------- */
1613 :
1614 : // OverviewList is historical, and quite inefficient, since it uses
1615 : // the full resolution source dataset, so only build it afterwards.
1616 : const CPLStringList aosOverviews(
1617 147 : CSLTokenizeString(CPLGetXMLValue(psTree, "OverviewList", "")));
1618 147 : if (!aosOverviews.empty())
1619 2 : CreateImplicitOverviews();
1620 :
1621 151 : for (int iOverview = 0; iOverview < aosOverviews.size(); ++iOverview)
1622 : {
1623 4 : int nOvFactor = atoi(aosOverviews[iOverview]);
1624 :
1625 4 : if (nOvFactor > 0)
1626 4 : BuildOverviews("NEAREST", 1, &nOvFactor, 0, nullptr, nullptr,
1627 : nullptr,
1628 : /*papszOptions=*/nullptr);
1629 : else
1630 0 : CPLError(CE_Failure, CPLE_AppDefined,
1631 : "Bad value for overview factor : %s",
1632 : aosOverviews[iOverview]);
1633 : }
1634 :
1635 147 : return eErr;
1636 : }
1637 :
1638 : /************************************************************************/
1639 : /* SerializeToXML() */
1640 : /************************************************************************/
1641 :
1642 42 : CPLXMLNode *VRTWarpedDataset::SerializeToXML(const char *pszVRTPathIn)
1643 :
1644 : {
1645 42 : CPLXMLNode *psTree = VRTDataset::SerializeToXML(pszVRTPathIn);
1646 :
1647 42 : if (psTree == nullptr)
1648 0 : return psTree;
1649 :
1650 : /* -------------------------------------------------------------------- */
1651 : /* Set subclass. */
1652 : /* -------------------------------------------------------------------- */
1653 42 : CPLCreateXMLNode(CPLCreateXMLNode(psTree, CXT_Attribute, "subClass"),
1654 : CXT_Text, "VRTWarpedDataset");
1655 :
1656 : /* -------------------------------------------------------------------- */
1657 : /* Serialize the block size. */
1658 : /* -------------------------------------------------------------------- */
1659 42 : CPLCreateXMLElementAndValue(psTree, "BlockXSize",
1660 : CPLSPrintf("%d", m_nBlockXSize));
1661 42 : CPLCreateXMLElementAndValue(psTree, "BlockYSize",
1662 : CPLSPrintf("%d", m_nBlockYSize));
1663 :
1664 : /* -------------------------------------------------------------------- */
1665 : /* Serialize the overview list (only for non implicit overviews) */
1666 : /* -------------------------------------------------------------------- */
1667 42 : if (!m_apoOverviews.empty())
1668 : {
1669 3 : int nSrcDSOvrCount = 0;
1670 3 : if (m_poWarper != nullptr && m_poWarper->GetOptions() != nullptr &&
1671 9 : m_poWarper->GetOptions()->hSrcDS != nullptr &&
1672 3 : GDALGetRasterCount(m_poWarper->GetOptions()->hSrcDS) > 0)
1673 : {
1674 : nSrcDSOvrCount =
1675 3 : static_cast<GDALDataset *>(m_poWarper->GetOptions()->hSrcDS)
1676 3 : ->GetRasterBand(1)
1677 3 : ->GetOverviewCount();
1678 : }
1679 :
1680 3 : if (static_cast<int>(m_apoOverviews.size()) != nSrcDSOvrCount)
1681 : {
1682 2 : const size_t nLen = m_apoOverviews.size() * 8 + 10;
1683 2 : char *pszOverviewList = static_cast<char *>(CPLMalloc(nLen));
1684 2 : pszOverviewList[0] = '\0';
1685 6 : for (auto *poOverviewDS : m_apoOverviews)
1686 : {
1687 4 : if (poOverviewDS)
1688 : {
1689 : const int nOvFactor = static_cast<int>(
1690 4 : 0.5 +
1691 4 : GetRasterXSize() / static_cast<double>(
1692 4 : poOverviewDS->GetRasterXSize()));
1693 :
1694 4 : snprintf(pszOverviewList + strlen(pszOverviewList),
1695 4 : nLen - strlen(pszOverviewList), "%d ", nOvFactor);
1696 : }
1697 : }
1698 :
1699 2 : CPLCreateXMLElementAndValue(psTree, "OverviewList",
1700 : pszOverviewList);
1701 :
1702 2 : CPLFree(pszOverviewList);
1703 : }
1704 : }
1705 :
1706 : /* -------------------------------------------------------------------- */
1707 : /* Serialize source overview level. */
1708 : /* -------------------------------------------------------------------- */
1709 42 : if (m_nSrcOvrLevel != -2)
1710 : {
1711 0 : if (m_nSrcOvrLevel < -2)
1712 0 : CPLCreateXMLElementAndValue(
1713 : psTree, "SrcOvrLevel",
1714 0 : CPLSPrintf("AUTO%d", m_nSrcOvrLevel + 2));
1715 0 : else if (m_nSrcOvrLevel == -1)
1716 0 : CPLCreateXMLElementAndValue(psTree, "SrcOvrLevel", "NONE");
1717 : else
1718 0 : CPLCreateXMLElementAndValue(psTree, "SrcOvrLevel",
1719 : CPLSPrintf("%d", m_nSrcOvrLevel));
1720 : }
1721 :
1722 : /* ==================================================================== */
1723 : /* Serialize the warp options. */
1724 : /* ==================================================================== */
1725 42 : if (m_poWarper != nullptr)
1726 : {
1727 : /* --------------------------------------------------------------------
1728 : */
1729 : /* We reset the destination dataset name so it doesn't get */
1730 : /* written out in the serialize warp options. */
1731 : /* --------------------------------------------------------------------
1732 : */
1733 42 : char *const pszSavedName = CPLStrdup(GetDescription());
1734 42 : SetDescription("");
1735 :
1736 : CPLXMLNode *const psWOTree =
1737 42 : GDALSerializeWarpOptions(m_poWarper->GetOptions());
1738 42 : CPLAddXMLChild(psTree, psWOTree);
1739 :
1740 42 : SetDescription(pszSavedName);
1741 42 : CPLFree(pszSavedName);
1742 :
1743 : /* --------------------------------------------------------------------
1744 : */
1745 : /* We need to consider making the source dataset relative to */
1746 : /* the VRT file if possible. Adjust accordingly. */
1747 : /* --------------------------------------------------------------------
1748 : */
1749 42 : CPLXMLNode *psSDS = CPLGetXMLNode(psWOTree, "SourceDataset");
1750 42 : int bRelativeToVRT = FALSE;
1751 : VSIStatBufL sStat;
1752 :
1753 42 : if (VSIStatExL(psSDS->psChild->pszValue, &sStat,
1754 42 : VSI_STAT_EXISTS_FLAG) == 0)
1755 : {
1756 80 : std::string osVRTFilename = pszVRTPathIn;
1757 40 : std::string osSourceDataset = psSDS->psChild->pszValue;
1758 40 : char *pszCurDir = CPLGetCurrentDir();
1759 40 : if (CPLIsFilenameRelative(osSourceDataset.c_str()) &&
1760 40 : !CPLIsFilenameRelative(osVRTFilename.c_str()) &&
1761 : pszCurDir != nullptr)
1762 : {
1763 6 : osSourceDataset = CPLFormFilenameSafe(
1764 3 : pszCurDir, osSourceDataset.c_str(), nullptr);
1765 : }
1766 37 : else if (!CPLIsFilenameRelative(osSourceDataset.c_str()) &&
1767 37 : CPLIsFilenameRelative(osVRTFilename.c_str()) &&
1768 : pszCurDir != nullptr)
1769 : {
1770 8 : osVRTFilename = CPLFormFilenameSafe(
1771 4 : pszCurDir, osVRTFilename.c_str(), nullptr);
1772 : }
1773 40 : CPLFree(pszCurDir);
1774 40 : char *pszRelativePath = CPLStrdup(CPLExtractRelativePath(
1775 : osVRTFilename.c_str(), osSourceDataset.c_str(),
1776 : &bRelativeToVRT));
1777 :
1778 40 : CPLFree(psSDS->psChild->pszValue);
1779 40 : psSDS->psChild->pszValue = pszRelativePath;
1780 : }
1781 :
1782 42 : CPLCreateXMLNode(
1783 : CPLCreateXMLNode(psSDS, CXT_Attribute, "relativeToVRT"), CXT_Text,
1784 42 : bRelativeToVRT ? "1" : "0");
1785 : }
1786 :
1787 42 : return psTree;
1788 : }
1789 :
1790 : /************************************************************************/
1791 : /* GetBlockSize() */
1792 : /************************************************************************/
1793 :
1794 655 : void VRTWarpedDataset::GetBlockSize(int *pnBlockXSize, int *pnBlockYSize) const
1795 :
1796 : {
1797 655 : CPLAssert(nullptr != pnBlockXSize);
1798 655 : CPLAssert(nullptr != pnBlockYSize);
1799 :
1800 655 : *pnBlockXSize = m_nBlockXSize;
1801 655 : *pnBlockYSize = m_nBlockYSize;
1802 655 : }
1803 :
1804 : /************************************************************************/
1805 : /* ProcessBlock() */
1806 : /* */
1807 : /* Warp a single requested block, and then push each band of */
1808 : /* the result into the block cache. */
1809 : /************************************************************************/
1810 :
1811 963 : CPLErr VRTWarpedDataset::ProcessBlock(int iBlockX, int iBlockY)
1812 :
1813 : {
1814 963 : if (m_poWarper == nullptr)
1815 0 : return CE_Failure;
1816 :
1817 963 : const int nReqXOff = iBlockX * m_nBlockXSize;
1818 963 : const int nReqYOff = iBlockY * m_nBlockYSize;
1819 963 : const int nReqXSize = std::min(m_nBlockXSize, nRasterXSize - nReqXOff);
1820 963 : const int nReqYSize = std::min(m_nBlockYSize, nRasterYSize - nReqYOff);
1821 :
1822 : GByte *pabyDstBuffer = static_cast<GByte *>(
1823 963 : m_poWarper->CreateDestinationBuffer(nReqXSize, nReqYSize));
1824 :
1825 963 : if (pabyDstBuffer == nullptr)
1826 : {
1827 0 : return CE_Failure;
1828 : }
1829 :
1830 : /* -------------------------------------------------------------------- */
1831 : /* Warp into this buffer. */
1832 : /* -------------------------------------------------------------------- */
1833 :
1834 963 : const GDALWarpOptions *psWO = m_poWarper->GetOptions();
1835 : const CPLErr eErr =
1836 1926 : m_poWarper->WarpRegionToBuffer(nReqXOff, nReqYOff, nReqXSize, nReqYSize,
1837 963 : pabyDstBuffer, psWO->eWorkingDataType);
1838 :
1839 963 : if (eErr != CE_None)
1840 : {
1841 0 : m_poWarper->DestroyDestinationBuffer(pabyDstBuffer);
1842 0 : return eErr;
1843 : }
1844 :
1845 : /* -------------------------------------------------------------------- */
1846 : /* Copy out into cache blocks for each band. */
1847 : /* -------------------------------------------------------------------- */
1848 963 : const int nWordSize = GDALGetDataTypeSizeBytes(psWO->eWorkingDataType);
1849 3721 : for (int i = 0; i < psWO->nBandCount; i++)
1850 : {
1851 2758 : int nDstBand = psWO->panDstBands[i];
1852 2758 : if (GetRasterCount() < nDstBand)
1853 : {
1854 0 : continue;
1855 : }
1856 :
1857 2758 : GDALRasterBand *poBand = GetRasterBand(nDstBand);
1858 : GDALRasterBlock *poBlock =
1859 2758 : poBand->GetLockedBlockRef(iBlockX, iBlockY, TRUE);
1860 :
1861 2758 : const GByte *pabyDstBandBuffer =
1862 : pabyDstBuffer +
1863 2758 : static_cast<GPtrDiff_t>(i) * nReqXSize * nReqYSize * nWordSize;
1864 :
1865 2758 : if (poBlock != nullptr)
1866 : {
1867 2758 : if (poBlock->GetDataRef() != nullptr)
1868 : {
1869 2758 : if (nReqXSize == m_nBlockXSize && nReqYSize == m_nBlockYSize)
1870 : {
1871 2586 : GDALCopyWords64(
1872 2586 : pabyDstBandBuffer, psWO->eWorkingDataType, nWordSize,
1873 : poBlock->GetDataRef(), poBlock->GetDataType(),
1874 : GDALGetDataTypeSizeBytes(poBlock->GetDataType()),
1875 2586 : static_cast<GPtrDiff_t>(m_nBlockXSize) * m_nBlockYSize);
1876 : }
1877 : else
1878 : {
1879 : GByte *pabyBlock =
1880 172 : static_cast<GByte *>(poBlock->GetDataRef());
1881 : const int nDTSize =
1882 172 : GDALGetDataTypeSizeBytes(poBlock->GetDataType());
1883 16981 : for (int iY = 0; iY < nReqYSize; iY++)
1884 : {
1885 16809 : GDALCopyWords(
1886 16809 : pabyDstBandBuffer + static_cast<GPtrDiff_t>(iY) *
1887 16809 : nReqXSize * nWordSize,
1888 16809 : psWO->eWorkingDataType, nWordSize,
1889 16809 : pabyBlock + static_cast<GPtrDiff_t>(iY) *
1890 16809 : m_nBlockXSize * nDTSize,
1891 : poBlock->GetDataType(), nDTSize, nReqXSize);
1892 : }
1893 : }
1894 : }
1895 :
1896 2758 : poBlock->DropLock();
1897 : }
1898 : }
1899 :
1900 963 : m_poWarper->DestroyDestinationBuffer(pabyDstBuffer);
1901 :
1902 963 : return CE_None;
1903 : }
1904 :
1905 : /************************************************************************/
1906 : /* IRasterIO() */
1907 : /************************************************************************/
1908 :
1909 : // Specialized implementation of IRasterIO() that will be faster than
1910 : // using the VRTWarpedRasterBand::IReadBlock() method in situations where
1911 : // - a large enough chunk of data is requested at once
1912 : // - and multi-threaded warping is enabled (it only kicks in if the warped
1913 : // chunk is large enough) and/or when reading the source dataset is
1914 : // multi-threaded (e.g JP2KAK or JP2OpenJPEG driver).
1915 2545 : CPLErr VRTWarpedDataset::IRasterIO(
1916 : GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
1917 : void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
1918 : int nBandCount, BANDMAP_TYPE panBandMap, GSpacing nPixelSpace,
1919 : GSpacing nLineSpace, GSpacing nBandSpace, GDALRasterIOExtraArg *psExtraArg)
1920 : {
1921 2521 : const bool bWholeImage = nXOff == 0 && nYOff == 0 &&
1922 5066 : nXSize == nRasterXSize && nYSize == nRasterYSize;
1923 :
1924 5090 : if (eRWFlag == GF_Write ||
1925 : // For too small request fall back to the block-based approach to
1926 : // benefit from caching
1927 2545 : (!bWholeImage &&
1928 2314 : (nBufXSize <= m_nBlockXSize || nBufYSize <= m_nBlockYSize)) ||
1929 : // Or if we don't request all bands at once
1930 5331 : nBandCount < nBands ||
1931 241 : !CPLTestBool(
1932 : CPLGetConfigOption("GDAL_VRT_WARP_USE_DATASET_RASTERIO", "YES")))
1933 : {
1934 2312 : return GDALDataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
1935 : pData, nBufXSize, nBufYSize, eBufType,
1936 : nBandCount, panBandMap, nPixelSpace,
1937 2312 : nLineSpace, nBandSpace, psExtraArg);
1938 : }
1939 :
1940 : // Try overviews for sub-sampled requests
1941 233 : if (nBufXSize < nXSize || nBufYSize < nYSize)
1942 : {
1943 5 : int bTried = FALSE;
1944 5 : const CPLErr eErr = TryOverviewRasterIO(
1945 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
1946 : eBufType, nBandCount, panBandMap, nPixelSpace, nLineSpace,
1947 : nBandSpace, psExtraArg, &bTried);
1948 :
1949 5 : if (bTried)
1950 : {
1951 3 : return eErr;
1952 : }
1953 : }
1954 :
1955 230 : if (m_poWarper == nullptr)
1956 0 : return CE_Failure;
1957 :
1958 230 : const GDALWarpOptions *psWO = m_poWarper->GetOptions();
1959 :
1960 230 : if (nBufXSize != nXSize || nBufYSize != nYSize)
1961 : {
1962 2 : if (!bWholeImage || !GDALTransformHasFastClone(psWO->pTransformerArg))
1963 : {
1964 0 : return GDALDataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
1965 : pData, nBufXSize, nBufYSize, eBufType,
1966 : nBandCount, panBandMap, nPixelSpace,
1967 0 : nLineSpace, nBandSpace, psExtraArg);
1968 : }
1969 :
1970 : // Build a temporary dataset taking into account the rescaling
1971 2 : void *pTransformerArg = GDALCloneTransformer(psWO->pTransformerArg);
1972 2 : if (pTransformerArg == nullptr)
1973 : {
1974 0 : return GDALDataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
1975 : pData, nBufXSize, nBufYSize, eBufType,
1976 : nBandCount, panBandMap, nPixelSpace,
1977 0 : nLineSpace, nBandSpace, psExtraArg);
1978 : }
1979 :
1980 2 : GDALWarpOptions *psRescaledWO = GDALCloneWarpOptions(psWO);
1981 2 : psRescaledWO->hSrcDS = psWO->hSrcDS;
1982 2 : psRescaledWO->pfnTransformer = psWO->pfnTransformer;
1983 2 : psRescaledWO->pTransformerArg = pTransformerArg;
1984 :
1985 : // Rescale the output geotransform on the transformer.
1986 2 : double adfDstGeoTransform[6] = {0.0};
1987 2 : GDALGetTransformerDstGeoTransform(psRescaledWO->pTransformerArg,
1988 : adfDstGeoTransform);
1989 2 : RescaleDstGeoTransform(adfDstGeoTransform, nRasterXSize, nBufXSize,
1990 : nRasterYSize, nBufYSize);
1991 2 : GDALSetTransformerDstGeoTransform(psRescaledWO->pTransformerArg,
1992 : adfDstGeoTransform);
1993 :
1994 : GDALDatasetH hDstDS =
1995 2 : GDALCreateWarpedVRT(psWO->hSrcDS, nBufXSize, nBufYSize,
1996 : adfDstGeoTransform, psRescaledWO);
1997 :
1998 2 : GDALDestroyWarpOptions(psRescaledWO);
1999 :
2000 2 : if (hDstDS == nullptr)
2001 : {
2002 : // Not supposed to happen in nominal circumstances. Could perhaps
2003 : // happen if some memory allocation error occurred in code called
2004 : // by GDALCreateWarpedVRT()
2005 0 : GDALDestroyTransformer(pTransformerArg);
2006 0 : return GDALDataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
2007 : pData, nBufXSize, nBufYSize, eBufType,
2008 : nBandCount, panBandMap, nPixelSpace,
2009 0 : nLineSpace, nBandSpace, psExtraArg);
2010 : }
2011 :
2012 2 : auto poOvrDS = static_cast<VRTWarpedDataset *>(hDstDS);
2013 2 : poOvrDS->m_bIsOverview = true;
2014 :
2015 : GDALRasterIOExtraArg sExtraArg;
2016 2 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
2017 2 : CPLErr eErr = poOvrDS->IRasterIO(GF_Read, 0, 0, nBufXSize, nBufYSize,
2018 : pData, nBufXSize, nBufYSize, eBufType,
2019 : nBandCount, panBandMap, nPixelSpace,
2020 : nLineSpace, nBandSpace, &sExtraArg);
2021 :
2022 2 : poOvrDS->ReleaseRef();
2023 2 : return eErr;
2024 : }
2025 :
2026 : // Build a map from warped output bands to their index
2027 456 : std::map<int, int> oMapBandToWarpingBandIndex;
2028 228 : bool bAllBandsIncreasingOrder =
2029 228 : (psWO->nBandCount == nBands && nBands == nBandCount);
2030 520 : for (int i = 0; i < psWO->nBandCount; ++i)
2031 : {
2032 292 : oMapBandToWarpingBandIndex[psWO->panDstBands[i]] = i;
2033 292 : if (psWO->panDstBands[i] != i + 1 || panBandMap[i] != i + 1)
2034 : {
2035 3 : bAllBandsIncreasingOrder = false;
2036 : }
2037 : }
2038 :
2039 : // Check that all requested bands are actually warped output bands.
2040 520 : for (int i = 0; i < nBandCount; ++i)
2041 : {
2042 316 : const int nRasterIOBand = panBandMap[i];
2043 316 : if (oMapBandToWarpingBandIndex.find(nRasterIOBand) ==
2044 632 : oMapBandToWarpingBandIndex.end())
2045 : {
2046 : // Not sure if that can happen...
2047 : // but if that does, that will likely later fail in ProcessBlock()
2048 24 : return GDALDataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
2049 : pData, nBufXSize, nBufYSize, eBufType,
2050 : nBandCount, panBandMap, nPixelSpace,
2051 24 : nLineSpace, nBandSpace, psExtraArg);
2052 : }
2053 : }
2054 :
2055 204 : int nSrcXOff = 0;
2056 204 : int nSrcYOff = 0;
2057 204 : int nSrcXSize = 0;
2058 204 : int nSrcYSize = 0;
2059 204 : double dfSrcXExtraSize = 0;
2060 204 : double dfSrcYExtraSize = 0;
2061 204 : double dfSrcFillRatio = 0;
2062 : // Find the source window that corresponds to our target window
2063 204 : if (m_poWarper->ComputeSourceWindow(nXOff, nYOff, nXSize, nYSize, &nSrcXOff,
2064 : &nSrcYOff, &nSrcXSize, &nSrcYSize,
2065 : &dfSrcXExtraSize, &dfSrcYExtraSize,
2066 204 : &dfSrcFillRatio) != CE_None)
2067 : {
2068 0 : return CE_Failure;
2069 : }
2070 :
2071 204 : GByte *const pabyDst = static_cast<GByte *>(pData);
2072 204 : const int nWarpDTSize = GDALGetDataTypeSizeBytes(psWO->eWorkingDataType);
2073 :
2074 204 : const double dfMemRequired = m_poWarper->GetWorkingMemoryForWindow(
2075 : nSrcXSize, nSrcYSize, nXSize, nYSize);
2076 : // If we need more warp working memory than allowed, we have to use a
2077 : // splitting strategy until we get below the limit.
2078 204 : if (dfMemRequired > psWO->dfWarpMemoryLimit && nXSize >= 2 && nYSize >= 2)
2079 : {
2080 3 : CPLDebugOnly("VRT", "VRTWarpedDataset::IRasterIO(): exceeding warp "
2081 : "memory. Splitting region");
2082 :
2083 : GDALRasterIOExtraArg sExtraArg;
2084 3 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
2085 :
2086 : bool bOK;
2087 : // Split along the longest dimension
2088 3 : if (nXSize >= nYSize)
2089 : {
2090 2 : const int nHalfXSize = nXSize / 2;
2091 4 : bOK = IRasterIO(GF_Read, nXOff, nYOff, nHalfXSize, nYSize, pabyDst,
2092 : nHalfXSize, nYSize, eBufType, nBandCount,
2093 : panBandMap, nPixelSpace, nLineSpace, nBandSpace,
2094 4 : &sExtraArg) == CE_None &&
2095 2 : IRasterIO(GF_Read, nXOff + nHalfXSize, nYOff,
2096 : nXSize - nHalfXSize, nYSize,
2097 2 : pabyDst + nHalfXSize * nPixelSpace,
2098 : nXSize - nHalfXSize, nYSize, eBufType, nBandCount,
2099 : panBandMap, nPixelSpace, nLineSpace, nBandSpace,
2100 : &sExtraArg) == CE_None;
2101 : }
2102 : else
2103 : {
2104 1 : const int nHalfYSize = nYSize / 2;
2105 2 : bOK = IRasterIO(GF_Read, nXOff, nYOff, nXSize, nHalfYSize, pabyDst,
2106 : nXSize, nHalfYSize, eBufType, nBandCount,
2107 : panBandMap, nPixelSpace, nLineSpace, nBandSpace,
2108 2 : &sExtraArg) == CE_None &&
2109 1 : IRasterIO(GF_Read, nXOff, nYOff + nHalfYSize, nXSize,
2110 : nYSize - nHalfYSize,
2111 1 : pabyDst + nHalfYSize * nLineSpace, nXSize,
2112 : nYSize - nHalfYSize, eBufType, nBandCount,
2113 : panBandMap, nPixelSpace, nLineSpace, nBandSpace,
2114 : &sExtraArg) == CE_None;
2115 : }
2116 3 : return bOK ? CE_None : CE_Failure;
2117 : }
2118 :
2119 201 : CPLDebugOnly("VRT",
2120 : "Using optimized VRTWarpedDataset::IRasterIO() code path");
2121 :
2122 : // Allocate a warping destination buffer if needed.
2123 : // We can use directly the output buffer pData if:
2124 : // - we request exactly all warped bands, and that there are as many
2125 : // warped bands as dataset bands (no alpha)
2126 : // - the output buffer data atype is the warping working data type
2127 : // - the output buffer has a band-sequential layout.
2128 : GByte *pabyWarpBuffer;
2129 :
2130 199 : if (bAllBandsIncreasingOrder && psWO->eWorkingDataType == eBufType &&
2131 112 : nPixelSpace == GDALGetDataTypeSizeBytes(eBufType) &&
2132 500 : nLineSpace == nPixelSpace * nXSize &&
2133 100 : (nBands == 1 || nBandSpace == nLineSpace * nYSize))
2134 : {
2135 100 : pabyWarpBuffer = static_cast<GByte *>(pData);
2136 100 : m_poWarper->InitializeDestinationBuffer(pabyWarpBuffer, nXSize, nYSize);
2137 : }
2138 : else
2139 : {
2140 : pabyWarpBuffer = static_cast<GByte *>(
2141 101 : m_poWarper->CreateDestinationBuffer(nXSize, nYSize));
2142 :
2143 101 : if (pabyWarpBuffer == nullptr)
2144 : {
2145 0 : return CE_Failure;
2146 : }
2147 : }
2148 :
2149 402 : const CPLErr eErr = m_poWarper->WarpRegionToBuffer(
2150 201 : nXOff, nYOff, nXSize, nYSize, pabyWarpBuffer, psWO->eWorkingDataType,
2151 : nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize, dfSrcXExtraSize,
2152 : dfSrcYExtraSize);
2153 :
2154 201 : if (pabyWarpBuffer != pData)
2155 : {
2156 101 : if (eErr == CE_None)
2157 : {
2158 : // Copy warping buffer into user destination buffer
2159 210 : for (int i = 0; i < nBandCount; i++)
2160 : {
2161 109 : const int nRasterIOBand = panBandMap[i];
2162 : const auto oIterToWarpingBandIndex =
2163 109 : oMapBandToWarpingBandIndex.find(nRasterIOBand);
2164 : // cannot happen due to earlier check
2165 109 : CPLAssert(oIterToWarpingBandIndex !=
2166 : oMapBandToWarpingBandIndex.end());
2167 :
2168 : const GByte *const pabyWarpBandBuffer =
2169 : pabyWarpBuffer +
2170 109 : static_cast<GPtrDiff_t>(oIterToWarpingBandIndex->second) *
2171 109 : nXSize * nYSize * nWarpDTSize;
2172 109 : GByte *const pabyDstBand = pabyDst + i * nBandSpace;
2173 :
2174 17342 : for (int iY = 0; iY < nYSize; iY++)
2175 : {
2176 17233 : GDALCopyWords(pabyWarpBandBuffer +
2177 17233 : static_cast<GPtrDiff_t>(iY) * nXSize *
2178 17233 : nWarpDTSize,
2179 17233 : psWO->eWorkingDataType, nWarpDTSize,
2180 17233 : pabyDstBand + iY * nLineSpace, eBufType,
2181 : static_cast<int>(nPixelSpace), nXSize);
2182 : }
2183 : }
2184 : }
2185 :
2186 101 : m_poWarper->DestroyDestinationBuffer(pabyWarpBuffer);
2187 : }
2188 :
2189 201 : return eErr;
2190 : }
2191 :
2192 : /************************************************************************/
2193 : /* AddBand() */
2194 : /************************************************************************/
2195 :
2196 394 : CPLErr VRTWarpedDataset::AddBand(GDALDataType eType,
2197 : CSLConstList /* papszOptions */)
2198 :
2199 : {
2200 394 : if (eType == GDT_Unknown || eType == GDT_TypeCount)
2201 : {
2202 1 : ReportError(CE_Failure, CPLE_IllegalArg,
2203 : "Illegal GDT_Unknown/GDT_TypeCount argument");
2204 1 : return CE_Failure;
2205 : }
2206 :
2207 393 : SetBand(GetRasterCount() + 1,
2208 393 : new VRTWarpedRasterBand(this, GetRasterCount() + 1, eType));
2209 :
2210 393 : return CE_None;
2211 : }
2212 :
2213 : /************************************************************************/
2214 : /* ==================================================================== */
2215 : /* VRTWarpedRasterBand */
2216 : /* ==================================================================== */
2217 : /************************************************************************/
2218 :
2219 : /************************************************************************/
2220 : /* VRTWarpedRasterBand() */
2221 : /************************************************************************/
2222 :
2223 655 : VRTWarpedRasterBand::VRTWarpedRasterBand(GDALDataset *poDSIn, int nBandIn,
2224 655 : GDALDataType eType)
2225 :
2226 : {
2227 655 : Initialize(poDSIn->GetRasterXSize(), poDSIn->GetRasterYSize());
2228 :
2229 655 : poDS = poDSIn;
2230 655 : nBand = nBandIn;
2231 655 : eAccess = GA_Update;
2232 :
2233 655 : static_cast<VRTWarpedDataset *>(poDS)->GetBlockSize(&nBlockXSize,
2234 : &nBlockYSize);
2235 :
2236 655 : if (eType != GDT_Unknown)
2237 397 : eDataType = eType;
2238 655 : }
2239 :
2240 : /************************************************************************/
2241 : /* ~VRTWarpedRasterBand() */
2242 : /************************************************************************/
2243 :
2244 1310 : VRTWarpedRasterBand::~VRTWarpedRasterBand()
2245 :
2246 : {
2247 655 : FlushCache(true);
2248 1310 : }
2249 :
2250 : /************************************************************************/
2251 : /* IReadBlock() */
2252 : /************************************************************************/
2253 :
2254 963 : CPLErr VRTWarpedRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
2255 : void *pImage)
2256 :
2257 : {
2258 963 : VRTWarpedDataset *poWDS = static_cast<VRTWarpedDataset *>(poDS);
2259 : const GPtrDiff_t nDataBytes =
2260 963 : static_cast<GPtrDiff_t>(GDALGetDataTypeSizeBytes(eDataType)) *
2261 963 : nBlockXSize * nBlockYSize;
2262 :
2263 963 : GDALRasterBlock *poBlock = GetLockedBlockRef(nBlockXOff, nBlockYOff, TRUE);
2264 963 : if (poBlock == nullptr)
2265 0 : return CE_Failure;
2266 :
2267 963 : if (poWDS->m_poWarper)
2268 : {
2269 963 : const GDALWarpOptions *psWO = poWDS->m_poWarper->GetOptions();
2270 963 : if (nBand == psWO->nDstAlphaBand)
2271 : {
2272 : // For a reader starting by asking on band 1, we should normally
2273 : // not reach here, because ProcessBlock() on band 1 will have
2274 : // populated the block cache for the regular bands and the alpha
2275 : // band.
2276 : // But if there's no source window corresponding to the block,
2277 : // the alpha band block will not be written through RasterIO(),
2278 : // so we nee to initialize it.
2279 34 : memset(poBlock->GetDataRef(), 0, nDataBytes);
2280 : }
2281 : }
2282 :
2283 963 : const CPLErr eErr = poWDS->ProcessBlock(nBlockXOff, nBlockYOff);
2284 :
2285 963 : if (eErr == CE_None && pImage != poBlock->GetDataRef())
2286 : {
2287 0 : memcpy(pImage, poBlock->GetDataRef(), nDataBytes);
2288 : }
2289 :
2290 963 : poBlock->DropLock();
2291 :
2292 963 : return eErr;
2293 : }
2294 :
2295 : /************************************************************************/
2296 : /* IWriteBlock() */
2297 : /************************************************************************/
2298 :
2299 806 : CPLErr VRTWarpedRasterBand::IWriteBlock(int nBlockXOff, int nBlockYOff,
2300 : void *pImage)
2301 :
2302 : {
2303 806 : VRTWarpedDataset *poWDS = static_cast<VRTWarpedDataset *>(poDS);
2304 :
2305 : // This is a bit tricky. In the case we are warping a VRTWarpedDataset
2306 : // with a destination alpha band, IWriteBlock can be called on that alpha
2307 : // band by GDALWarpDstAlphaMasker
2308 : // We don't need to do anything since the data will have hopefully been
2309 : // read from the block cache before if the reader processes all the bands
2310 : // of a same block.
2311 806 : if (poWDS->m_poWarper->GetOptions()->nDstAlphaBand != nBand)
2312 : {
2313 : /* Otherwise, call the superclass method, that will fail of course */
2314 0 : return VRTRasterBand::IWriteBlock(nBlockXOff, nBlockYOff, pImage);
2315 : }
2316 :
2317 806 : return CE_None;
2318 : }
2319 :
2320 : /************************************************************************/
2321 : /* EmitErrorMessageIfWriteNotSupported() */
2322 : /************************************************************************/
2323 :
2324 808 : bool VRTWarpedRasterBand::EmitErrorMessageIfWriteNotSupported(
2325 : const char *pszCaller) const
2326 : {
2327 808 : VRTWarpedDataset *poWDS = static_cast<VRTWarpedDataset *>(poDS);
2328 : // Cf comment in IWriteBlock()
2329 808 : if (poWDS->m_poWarper->GetOptions()->nDstAlphaBand != nBand)
2330 : {
2331 2 : ReportError(CE_Failure, CPLE_NoWriteAccess,
2332 : "%s: attempt to write to a VRTWarpedRasterBand.",
2333 : pszCaller);
2334 :
2335 2 : return true;
2336 : }
2337 806 : return false;
2338 : }
2339 :
2340 : /************************************************************************/
2341 : /* GetBestOverviewLevel() */
2342 : /************************************************************************/
2343 :
2344 0 : int VRTWarpedRasterBand::GetBestOverviewLevel(
2345 : int &nXOff, int &nYOff, int &nXSize, int &nYSize, int nBufXSize,
2346 : int nBufYSize, GDALRasterIOExtraArg *psExtraArg) const
2347 : {
2348 0 : VRTWarpedDataset *poWDS = static_cast<VRTWarpedDataset *>(poDS);
2349 :
2350 : /* -------------------------------------------------------------------- */
2351 : /* Compute the desired downsampling factor. It is */
2352 : /* based on the least reduced axis, and represents the number */
2353 : /* of source pixels to one destination pixel. */
2354 : /* -------------------------------------------------------------------- */
2355 0 : const double dfDesiredDownsamplingFactor =
2356 0 : ((nXSize / static_cast<double>(nBufXSize)) <
2357 0 : (nYSize / static_cast<double>(nBufYSize)) ||
2358 : nBufYSize == 1)
2359 0 : ? nXSize / static_cast<double>(nBufXSize)
2360 0 : : nYSize / static_cast<double>(nBufYSize);
2361 :
2362 : /* -------------------------------------------------------------------- */
2363 : /* Find the overview level that largest downsampling factor (most */
2364 : /* downsampled) that is still less than (or only a little more) */
2365 : /* downsampled than the request. */
2366 : /* -------------------------------------------------------------------- */
2367 0 : const GDALWarpOptions *psWO = poWDS->m_poWarper->GetOptions();
2368 0 : GDALDataset *poSrcDS = GDALDataset::FromHandle(psWO->hSrcDS);
2369 0 : const int nOverviewCount = poSrcDS->GetRasterBand(1)->GetOverviewCount();
2370 :
2371 0 : int nBestOverviewXSize = 1;
2372 0 : int nBestOverviewYSize = 1;
2373 0 : double dfBestDownsamplingFactor = 0;
2374 0 : int nBestOverviewLevel = -1;
2375 :
2376 : const char *pszOversampligThreshold =
2377 0 : CPLGetConfigOption("GDAL_OVERVIEW_OVERSAMPLING_THRESHOLD", nullptr);
2378 :
2379 : // Cf https://github.com/OSGeo/gdal/pull/9040#issuecomment-1898524693
2380 : // Do not exactly use a oversampling threshold of 1.0 because of numerical
2381 : // instability.
2382 0 : const auto AdjustThreshold = [](double x)
2383 : {
2384 0 : constexpr double EPS = 1e-2;
2385 0 : return x == 1.0 ? x + EPS : x;
2386 : };
2387 0 : const double dfOversamplingThreshold = AdjustThreshold(
2388 0 : pszOversampligThreshold ? CPLAtof(pszOversampligThreshold)
2389 0 : : psExtraArg && psExtraArg->eResampleAlg != GRIORA_NearestNeighbour
2390 0 : ? 1.0
2391 : : 1.2);
2392 0 : for (int iOverview = 0; iOverview < nOverviewCount; iOverview++)
2393 : {
2394 0 : const GDALRasterBand *poSrcOvrBand = this;
2395 0 : bool bThisLevelOnly = false;
2396 : const int iSrcOvr =
2397 0 : poWDS->GetSrcOverviewLevel(iOverview, bThisLevelOnly);
2398 0 : if (iSrcOvr >= 0)
2399 : {
2400 0 : poSrcOvrBand = poSrcDS->GetRasterBand(1)->GetOverview(iSrcOvr);
2401 : }
2402 0 : if (poSrcOvrBand == nullptr)
2403 0 : break;
2404 :
2405 0 : int nDstPixels = 0;
2406 0 : int nDstLines = 0;
2407 0 : double dfSrcRatioX = 0;
2408 0 : double dfSrcRatioY = 0;
2409 0 : if (!poWDS->GetOverviewSize(poSrcDS, iOverview, iSrcOvr, nDstPixels,
2410 : nDstLines, dfSrcRatioX, dfSrcRatioY))
2411 : {
2412 0 : break;
2413 : }
2414 :
2415 : // Compute downsampling factor of this overview
2416 : const double dfDownsamplingFactor =
2417 0 : std::min(nRasterXSize / static_cast<double>(nDstPixels),
2418 0 : nRasterYSize / static_cast<double>(nDstLines));
2419 :
2420 : // Is it nearly the requested factor and better (lower) than
2421 : // the current best factor?
2422 0 : if (dfDownsamplingFactor >=
2423 0 : dfDesiredDownsamplingFactor * dfOversamplingThreshold ||
2424 : dfDownsamplingFactor <= dfBestDownsamplingFactor)
2425 : {
2426 0 : continue;
2427 : }
2428 :
2429 : // Ignore AVERAGE_BIT2GRAYSCALE overviews for RasterIO purposes.
2430 : const char *pszResampling = const_cast<GDALRasterBand *>(poSrcOvrBand)
2431 0 : ->GetMetadataItem("RESAMPLING");
2432 :
2433 0 : if (pszResampling != nullptr &&
2434 0 : STARTS_WITH_CI(pszResampling, "AVERAGE_BIT2"))
2435 0 : continue;
2436 :
2437 : // OK, this is our new best overview.
2438 0 : nBestOverviewXSize = nDstPixels;
2439 0 : nBestOverviewYSize = nDstLines;
2440 0 : nBestOverviewLevel = iOverview;
2441 0 : dfBestDownsamplingFactor = dfDownsamplingFactor;
2442 : }
2443 :
2444 : /* -------------------------------------------------------------------- */
2445 : /* If we didn't find an overview that helps us, just return */
2446 : /* indicating failure and the full resolution image will be used. */
2447 : /* -------------------------------------------------------------------- */
2448 0 : if (nBestOverviewLevel < 0)
2449 0 : return -1;
2450 :
2451 : /* -------------------------------------------------------------------- */
2452 : /* Recompute the source window in terms of the selected */
2453 : /* overview. */
2454 : /* -------------------------------------------------------------------- */
2455 0 : const double dfXFactor =
2456 0 : nRasterXSize / static_cast<double>(nBestOverviewXSize);
2457 0 : const double dfYFactor =
2458 0 : nRasterYSize / static_cast<double>(nBestOverviewYSize);
2459 0 : CPLDebug("GDAL", "Selecting overview %d x %d", nBestOverviewXSize,
2460 : nBestOverviewYSize);
2461 :
2462 0 : const int nOXOff = std::min(nBestOverviewXSize - 1,
2463 0 : static_cast<int>(nXOff / dfXFactor + 0.5));
2464 0 : const int nOYOff = std::min(nBestOverviewYSize - 1,
2465 0 : static_cast<int>(nYOff / dfYFactor + 0.5));
2466 0 : int nOXSize = std::max(1, static_cast<int>(nXSize / dfXFactor + 0.5));
2467 0 : int nOYSize = std::max(1, static_cast<int>(nYSize / dfYFactor + 0.5));
2468 0 : if (nOXOff + nOXSize > nBestOverviewXSize)
2469 0 : nOXSize = nBestOverviewXSize - nOXOff;
2470 0 : if (nOYOff + nOYSize > nBestOverviewYSize)
2471 0 : nOYSize = nBestOverviewYSize - nOYOff;
2472 :
2473 0 : if (psExtraArg)
2474 : {
2475 0 : if (psExtraArg->bFloatingPointWindowValidity)
2476 : {
2477 0 : psExtraArg->dfXOff /= dfXFactor;
2478 0 : psExtraArg->dfXSize /= dfXFactor;
2479 0 : psExtraArg->dfYOff /= dfYFactor;
2480 0 : psExtraArg->dfYSize /= dfYFactor;
2481 : }
2482 0 : else if (psExtraArg->eResampleAlg != GRIORA_NearestNeighbour)
2483 : {
2484 0 : psExtraArg->bFloatingPointWindowValidity = true;
2485 0 : psExtraArg->dfXOff = nXOff / dfXFactor;
2486 0 : psExtraArg->dfXSize = nXSize / dfXFactor;
2487 0 : psExtraArg->dfYOff = nYOff / dfYFactor;
2488 0 : psExtraArg->dfYSize = nYSize / dfYFactor;
2489 : }
2490 : }
2491 :
2492 0 : nXOff = nOXOff;
2493 0 : nYOff = nOYOff;
2494 0 : nXSize = nOXSize;
2495 0 : nYSize = nOYSize;
2496 :
2497 0 : return nBestOverviewLevel;
2498 : }
2499 :
2500 : /************************************************************************/
2501 : /* IRasterIO() */
2502 : /************************************************************************/
2503 :
2504 12327 : CPLErr VRTWarpedRasterBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
2505 : int nXSize, int nYSize, void *pData,
2506 : int nBufXSize, int nBufYSize,
2507 : GDALDataType eBufType,
2508 : GSpacing nPixelSpace, GSpacing nLineSpace,
2509 : GDALRasterIOExtraArg *psExtraArg)
2510 : {
2511 12327 : VRTWarpedDataset *poWDS = static_cast<VRTWarpedDataset *>(poDS);
2512 12327 : if (m_nIRasterIOCounter == 0 && poWDS->GetRasterCount() == 1)
2513 : {
2514 916 : int anBandMap[] = {nBand};
2515 916 : ++m_nIRasterIOCounter;
2516 916 : const CPLErr eErr = poWDS->IRasterIO(
2517 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
2518 : eBufType, 1, anBandMap, nPixelSpace, nLineSpace, 0, psExtraArg);
2519 916 : --m_nIRasterIOCounter;
2520 916 : return eErr;
2521 : }
2522 :
2523 : /* ==================================================================== */
2524 : /* Do we have overviews that would be appropriate to satisfy */
2525 : /* this request? */
2526 : /* ==================================================================== */
2527 11411 : if ((nBufXSize < nXSize || nBufYSize < nYSize) && GetOverviewCount() &&
2528 : eRWFlag == GF_Read)
2529 : {
2530 : GDALRasterIOExtraArg sExtraArg;
2531 0 : GDALCopyRasterIOExtraArg(&sExtraArg, psExtraArg);
2532 :
2533 0 : const int nOverview = GetBestOverviewLevel(
2534 : nXOff, nYOff, nXSize, nYSize, nBufXSize, nBufYSize, &sExtraArg);
2535 0 : if (nOverview >= 0)
2536 : {
2537 0 : auto poOvrBand = GetOverview(nOverview);
2538 0 : if (!poOvrBand)
2539 0 : return CE_Failure;
2540 :
2541 0 : return poOvrBand->RasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
2542 : pData, nBufXSize, nBufYSize, eBufType,
2543 0 : nPixelSpace, nLineSpace, &sExtraArg);
2544 : }
2545 : }
2546 :
2547 11411 : return GDALRasterBand::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
2548 : pData, nBufXSize, nBufYSize, eBufType,
2549 11411 : nPixelSpace, nLineSpace, psExtraArg);
2550 : }
2551 :
2552 : /************************************************************************/
2553 : /* SerializeToXML() */
2554 : /************************************************************************/
2555 :
2556 71 : CPLXMLNode *VRTWarpedRasterBand::SerializeToXML(const char *pszVRTPathIn,
2557 : bool &bHasWarnedAboutRAMUsage,
2558 : size_t &nAccRAMUsage)
2559 :
2560 : {
2561 71 : CPLXMLNode *const psTree = VRTRasterBand::SerializeToXML(
2562 : pszVRTPathIn, bHasWarnedAboutRAMUsage, nAccRAMUsage);
2563 :
2564 : /* -------------------------------------------------------------------- */
2565 : /* Set subclass. */
2566 : /* -------------------------------------------------------------------- */
2567 71 : CPLCreateXMLNode(CPLCreateXMLNode(psTree, CXT_Attribute, "subClass"),
2568 : CXT_Text, "VRTWarpedRasterBand");
2569 :
2570 71 : return psTree;
2571 : }
2572 :
2573 : /************************************************************************/
2574 : /* GetOverviewCount() */
2575 : /************************************************************************/
2576 :
2577 136 : int VRTWarpedRasterBand::GetOverviewCount()
2578 :
2579 : {
2580 136 : VRTWarpedDataset *const poWDS = static_cast<VRTWarpedDataset *>(poDS);
2581 136 : if (poWDS->m_bIsOverview)
2582 1 : return 0;
2583 :
2584 135 : if (poWDS->m_apoOverviews.empty())
2585 : {
2586 104 : return poWDS->GetOverviewCount();
2587 : }
2588 :
2589 31 : return static_cast<int>(poWDS->m_apoOverviews.size());
2590 : }
2591 :
2592 : /************************************************************************/
2593 : /* GetOverview() */
2594 : /************************************************************************/
2595 :
2596 45 : GDALRasterBand *VRTWarpedRasterBand::GetOverview(int iOverview)
2597 :
2598 : {
2599 45 : VRTWarpedDataset *const poWDS = static_cast<VRTWarpedDataset *>(poDS);
2600 :
2601 45 : const int nOvrCount = GetOverviewCount();
2602 45 : if (iOverview < 0 || iOverview >= nOvrCount)
2603 6 : return nullptr;
2604 :
2605 39 : if (poWDS->m_apoOverviews.empty())
2606 12 : poWDS->m_apoOverviews.resize(nOvrCount);
2607 39 : if (!poWDS->m_apoOverviews[iOverview])
2608 19 : poWDS->m_apoOverviews[iOverview] =
2609 19 : poWDS->CreateImplicitOverview(iOverview);
2610 39 : if (!poWDS->m_apoOverviews[iOverview])
2611 0 : return nullptr;
2612 39 : return poWDS->m_apoOverviews[iOverview]->GetRasterBand(nBand);
2613 : }
2614 :
2615 : /*! @endcond */
|
