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