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