Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL Utilities
4 : * Purpose: Command line application to build VRT datasets from raster products
5 : * or content of SHP tile index
6 : * Author: Even Rouault, <even dot rouault at spatialys dot com>
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2007-2016, Even Rouault <even dot rouault at spatialys dot com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "ogr_api.h"
15 : #include "ogr_srs_api.h"
16 :
17 : #include "cpl_port.h"
18 : #include "gdal_utils.h"
19 : #include "gdal_utils_priv.h"
20 : #include "gdalargumentparser.h"
21 :
22 : #include <cassert>
23 : #include <cmath>
24 : #include <cstdio>
25 : #include <cstdlib>
26 : #include <cstring>
27 :
28 : #include <algorithm>
29 : #include <memory>
30 : #include <optional>
31 : #include <set>
32 : #include <string>
33 : #include <vector>
34 :
35 : #include "commonutils.h"
36 : #include "cpl_conv.h"
37 : #include "cpl_error.h"
38 : #include "cpl_float.h"
39 : #include "cpl_progress.h"
40 : #include "cpl_string.h"
41 : #include "cpl_vsi.h"
42 : #include "cpl_vsi_virtual.h"
43 : #include "gdal.h"
44 : #include "gdal_vrt.h"
45 : #include "gdal_priv.h"
46 : #include "gdal_proxy.h"
47 : #include "ogr_api.h"
48 : #include "ogr_core.h"
49 : #include "ogr_srs_api.h"
50 : #include "ogr_spatialref.h"
51 : #include "ogrsf_frmts.h"
52 : #include "vrtdataset.h"
53 :
54 : #define GEOTRSFRM_TOPLEFT_X 0
55 : #define GEOTRSFRM_WE_RES 1
56 : #define GEOTRSFRM_ROTATION_PARAM1 2
57 : #define GEOTRSFRM_TOPLEFT_Y 3
58 : #define GEOTRSFRM_ROTATION_PARAM2 4
59 : #define GEOTRSFRM_NS_RES 5
60 :
61 : namespace gdal::GDALBuildVRT
62 : {
63 : typedef enum
64 : {
65 : LOWEST_RESOLUTION,
66 : HIGHEST_RESOLUTION,
67 : AVERAGE_RESOLUTION,
68 : SAME_RESOLUTION,
69 : USER_RESOLUTION,
70 : COMMON_RESOLUTION,
71 : } ResolutionStrategy;
72 :
73 : struct DatasetProperty
74 : {
75 : int isFileOK = FALSE;
76 : int nRasterXSize = 0;
77 : int nRasterYSize = 0;
78 : GDALGeoTransform gt{};
79 : int nBlockXSize = 0;
80 : int nBlockYSize = 0;
81 : std::vector<GDALDataType> aeBandType{};
82 : std::vector<bool> abHasNoData{};
83 : std::vector<double> adfNoDataValues{};
84 : std::vector<bool> abHasOffset{};
85 : std::vector<double> adfOffset{};
86 : std::vector<bool> abHasScale{};
87 : std::vector<bool> abHasMaskBand{};
88 : std::vector<double> adfScale{};
89 : int bHasDatasetMask = 0;
90 : bool bLastBandIsAlpha = false;
91 : int nMaskBlockXSize = 0;
92 : int nMaskBlockYSize = 0;
93 : std::vector<int> anOverviewFactors{};
94 : std::vector<std::string> aosDescriptions{};
95 : std::map<int, std::map<std::string, std::string>> mapBandMetadata{};
96 : };
97 :
98 : struct BandProperty
99 : {
100 : GDALColorInterp colorInterpretation = GCI_Undefined;
101 : GDALDataType dataType = GDT_Unknown;
102 : std::unique_ptr<GDALColorTable> colorTable{};
103 : bool bHasNoData = false;
104 : double noDataValue = 0;
105 : bool bHasOffset = false;
106 : double dfOffset = 0;
107 : bool bHasScale = false;
108 : double dfScale = 0;
109 : std::string osDescription = "";
110 : std::map<std::string, std::string> mapBandMetadata{};
111 : };
112 : } // namespace gdal::GDALBuildVRT
113 :
114 : using namespace gdal::GDALBuildVRT;
115 :
116 : /************************************************************************/
117 : /* GetSrcDstWin() */
118 : /************************************************************************/
119 :
120 2387 : static int GetSrcDstWin(DatasetProperty *psDP, double we_res, double ns_res,
121 : double minX, double minY, double maxX, double maxY,
122 : int nTargetXSize, int nTargetYSize, double *pdfSrcXOff,
123 : double *pdfSrcYOff, double *pdfSrcXSize,
124 : double *pdfSrcYSize, double *pdfDstXOff,
125 : double *pdfDstYOff, double *pdfDstXSize,
126 : double *pdfDstYSize)
127 : {
128 2387 : if (we_res == 0 || ns_res == 0)
129 : {
130 : // should not happen. to please Coverity
131 0 : return FALSE;
132 : }
133 :
134 : /* Check that the destination bounding box intersects the source bounding
135 : * box */
136 2387 : if (psDP->gt[GEOTRSFRM_TOPLEFT_X] +
137 2387 : psDP->nRasterXSize * psDP->gt[GEOTRSFRM_WE_RES] <=
138 : minX)
139 0 : return FALSE;
140 2387 : if (psDP->gt[GEOTRSFRM_TOPLEFT_X] >= maxX)
141 1 : return FALSE;
142 2386 : if (psDP->gt[GEOTRSFRM_TOPLEFT_Y] +
143 2386 : psDP->nRasterYSize * psDP->gt[GEOTRSFRM_NS_RES] >=
144 : maxY)
145 0 : return FALSE;
146 2386 : if (psDP->gt[GEOTRSFRM_TOPLEFT_Y] <= minY)
147 0 : return FALSE;
148 :
149 2386 : if (psDP->gt[GEOTRSFRM_TOPLEFT_X] < minX)
150 : {
151 4 : *pdfSrcXOff =
152 4 : (minX - psDP->gt[GEOTRSFRM_TOPLEFT_X]) / psDP->gt[GEOTRSFRM_WE_RES];
153 4 : *pdfDstXOff = 0.0;
154 : }
155 : else
156 : {
157 2382 : *pdfSrcXOff = 0.0;
158 2382 : *pdfDstXOff = ((psDP->gt[GEOTRSFRM_TOPLEFT_X] - minX) / we_res);
159 : }
160 2386 : if (maxY < psDP->gt[GEOTRSFRM_TOPLEFT_Y])
161 : {
162 7 : *pdfSrcYOff = (psDP->gt[GEOTRSFRM_TOPLEFT_Y] - maxY) /
163 7 : -psDP->gt[GEOTRSFRM_NS_RES];
164 7 : *pdfDstYOff = 0.0;
165 : }
166 : else
167 : {
168 2379 : *pdfSrcYOff = 0.0;
169 2379 : *pdfDstYOff = ((maxY - psDP->gt[GEOTRSFRM_TOPLEFT_Y]) / -ns_res);
170 : }
171 :
172 2386 : *pdfSrcXSize = psDP->nRasterXSize;
173 2386 : *pdfSrcYSize = psDP->nRasterYSize;
174 2386 : if (*pdfSrcXOff > 0)
175 4 : *pdfSrcXSize -= *pdfSrcXOff;
176 2386 : if (*pdfSrcYOff > 0)
177 7 : *pdfSrcYSize -= *pdfSrcYOff;
178 :
179 2386 : const double dfSrcToDstXSize = psDP->gt[GEOTRSFRM_WE_RES] / we_res;
180 2386 : *pdfDstXSize = *pdfSrcXSize * dfSrcToDstXSize;
181 2386 : const double dfSrcToDstYSize = psDP->gt[GEOTRSFRM_NS_RES] / ns_res;
182 2386 : *pdfDstYSize = *pdfSrcYSize * dfSrcToDstYSize;
183 :
184 2386 : if (*pdfDstXOff + *pdfDstXSize > nTargetXSize)
185 : {
186 8 : *pdfDstXSize = nTargetXSize - *pdfDstXOff;
187 8 : *pdfSrcXSize = *pdfDstXSize / dfSrcToDstXSize;
188 : }
189 :
190 2386 : if (*pdfDstYOff + *pdfDstYSize > nTargetYSize)
191 : {
192 7 : *pdfDstYSize = nTargetYSize - *pdfDstYOff;
193 7 : *pdfSrcYSize = *pdfDstYSize / dfSrcToDstYSize;
194 : }
195 :
196 4772 : return *pdfSrcXSize > 0 && *pdfDstXSize > 0 && *pdfSrcYSize > 0 &&
197 4772 : *pdfDstYSize > 0;
198 : }
199 :
200 : /************************************************************************/
201 : /* VRTBuilder */
202 : /************************************************************************/
203 :
204 : class VRTBuilder
205 : {
206 : /* Input parameters */
207 : bool bStrict = false;
208 : char *pszOutputFilename = nullptr;
209 : int nInputFiles = 0;
210 : char **ppszInputFilenames = nullptr;
211 : int nSrcDSCount = 0;
212 : GDALDatasetH *pahSrcDS = nullptr;
213 : int nTotalBands = 0;
214 : bool bLastBandIsAlpha = false;
215 : bool bExplicitBandList = false;
216 : int nMaxSelectedBandNo = 0;
217 : int nSelectedBands = 0;
218 : int *panSelectedBandList = nullptr;
219 : ResolutionStrategy resolutionStrategy = AVERAGE_RESOLUTION;
220 : int nCountValid = 0;
221 : double we_res = 0;
222 : double ns_res = 0;
223 : int bTargetAlignedPixels = 0;
224 : double minX = 0;
225 : double minY = 0;
226 : double maxX = 0;
227 : double maxY = 0;
228 : int bSeparate = 0;
229 : int bAllowProjectionDifference = 0;
230 : int bAddAlpha = 0;
231 : int bHideNoData = 0;
232 : int nSubdataset = 0;
233 : char *pszSrcNoData = nullptr;
234 : char *pszVRTNoData = nullptr;
235 : char *pszOutputSRS = nullptr;
236 : char *pszResampling = nullptr;
237 : char **papszOpenOptions = nullptr;
238 : bool bUseSrcMaskBand = true;
239 : bool bNoDataFromMask = false;
240 : double dfMaskValueThreshold = 0;
241 : const CPLStringList aosCreateOptions;
242 : std::string osPixelFunction{};
243 : const CPLStringList aosPixelFunctionArgs;
244 : const bool bWriteAbsolutePath;
245 :
246 : /* Internal variables */
247 : char *pszProjectionRef = nullptr;
248 : std::vector<BandProperty> asBandProperties{};
249 : int bFirst = TRUE;
250 : int bHasGeoTransform = 0;
251 : int nRasterXSize = 0;
252 : int nRasterYSize = 0;
253 : std::vector<DatasetProperty> asDatasetProperties{};
254 : int bUserExtent = 0;
255 : int bAllowSrcNoData = TRUE;
256 : double *padfSrcNoData = nullptr;
257 : int nSrcNoDataCount = 0;
258 : int bAllowVRTNoData = TRUE;
259 : double *padfVRTNoData = nullptr;
260 : int nVRTNoDataCount = 0;
261 : int bHasRunBuild = 0;
262 : int bHasDatasetMask = 0;
263 :
264 : std::string AnalyseRaster(GDALDatasetH hDS,
265 : DatasetProperty *psDatasetProperties);
266 :
267 : void CreateVRTSeparate(VRTDataset *poVTDS);
268 : void CreateVRTNonSeparate(VRTDataset *poVRTDS);
269 :
270 : CPL_DISALLOW_COPY_ASSIGN(VRTBuilder)
271 :
272 : public:
273 : VRTBuilder(bool bStrictIn, const char *pszOutputFilename, int nInputFiles,
274 : const char *const *ppszInputFilenames, GDALDatasetH *pahSrcDSIn,
275 : const int *panSelectedBandListIn, int nBandCount,
276 : ResolutionStrategy resolutionStrategy, double we_res,
277 : double ns_res, int bTargetAlignedPixels, double minX,
278 : double minY, double maxX, double maxY, int bSeparate,
279 : int bAllowProjectionDifference, int bAddAlpha, int bHideNoData,
280 : int nSubdataset, const char *pszSrcNoData,
281 : const char *pszVRTNoData, bool bUseSrcMaskBand,
282 : bool bNoDataFromMask, double dfMaskValueThreshold,
283 : const char *pszOutputSRS, const char *pszResampling,
284 : const char *pszPixelFunctionName,
285 : const CPLStringList &aosPixelFunctionArgs,
286 : const char *const *papszOpenOptionsIn,
287 : const CPLStringList &aosCreateOptionsIn,
288 : bool bWriteAbsolutePathIn);
289 :
290 : ~VRTBuilder();
291 :
292 : std::unique_ptr<GDALDataset> Build(GDALProgressFunc pfnProgress,
293 : void *pProgressData);
294 :
295 : std::string m_osProgramName{};
296 : };
297 :
298 : /************************************************************************/
299 : /* VRTBuilder() */
300 : /************************************************************************/
301 :
302 252 : VRTBuilder::VRTBuilder(
303 : bool bStrictIn, const char *pszOutputFilenameIn, int nInputFilesIn,
304 : const char *const *ppszInputFilenamesIn, GDALDatasetH *pahSrcDSIn,
305 : const int *panSelectedBandListIn, int nBandCount,
306 : ResolutionStrategy resolutionStrategyIn, double we_resIn, double ns_resIn,
307 : int bTargetAlignedPixelsIn, double minXIn, double minYIn, double maxXIn,
308 : double maxYIn, int bSeparateIn, int bAllowProjectionDifferenceIn,
309 : int bAddAlphaIn, int bHideNoDataIn, int nSubdatasetIn,
310 : const char *pszSrcNoDataIn, const char *pszVRTNoDataIn,
311 : bool bUseSrcMaskBandIn, bool bNoDataFromMaskIn,
312 : double dfMaskValueThresholdIn, const char *pszOutputSRSIn,
313 : const char *pszResamplingIn, const char *pszPixelFunctionIn,
314 : const CPLStringList &aosPixelFunctionArgsIn,
315 : const char *const *papszOpenOptionsIn,
316 252 : const CPLStringList &aosCreateOptionsIn, bool bWriteAbsolutePathIn)
317 : : bStrict(bStrictIn), aosCreateOptions(aosCreateOptionsIn),
318 : aosPixelFunctionArgs(aosPixelFunctionArgsIn),
319 252 : bWriteAbsolutePath(bWriteAbsolutePathIn)
320 : {
321 252 : pszOutputFilename = CPLStrdup(pszOutputFilenameIn);
322 252 : nInputFiles = nInputFilesIn;
323 252 : papszOpenOptions = CSLDuplicate(const_cast<char **>(papszOpenOptionsIn));
324 :
325 252 : if (pszPixelFunctionIn != nullptr)
326 : {
327 3 : osPixelFunction = pszPixelFunctionIn;
328 : }
329 :
330 252 : if (ppszInputFilenamesIn)
331 : {
332 158 : ppszInputFilenames =
333 158 : static_cast<char **>(CPLMalloc(nInputFiles * sizeof(char *)));
334 1417 : for (int i = 0; i < nInputFiles; i++)
335 : {
336 1259 : ppszInputFilenames[i] = CPLStrdup(ppszInputFilenamesIn[i]);
337 : }
338 : }
339 94 : else if (pahSrcDSIn)
340 : {
341 94 : nSrcDSCount = nInputFiles;
342 94 : pahSrcDS = static_cast<GDALDatasetH *>(
343 94 : CPLMalloc(nInputFiles * sizeof(GDALDatasetH)));
344 94 : memcpy(pahSrcDS, pahSrcDSIn, nInputFiles * sizeof(GDALDatasetH));
345 94 : ppszInputFilenames =
346 94 : static_cast<char **>(CPLMalloc(nInputFiles * sizeof(char *)));
347 1241 : for (int i = 0; i < nInputFiles; i++)
348 : {
349 2294 : ppszInputFilenames[i] =
350 1147 : CPLStrdup(GDALGetDescription(pahSrcDSIn[i]));
351 : }
352 : }
353 :
354 252 : bExplicitBandList = nBandCount != 0;
355 252 : nSelectedBands = nBandCount;
356 252 : if (nBandCount)
357 : {
358 19 : panSelectedBandList =
359 19 : static_cast<int *>(CPLMalloc(nSelectedBands * sizeof(int)));
360 19 : memcpy(panSelectedBandList, panSelectedBandListIn,
361 19 : nSelectedBands * sizeof(int));
362 : }
363 :
364 252 : resolutionStrategy = resolutionStrategyIn;
365 252 : we_res = we_resIn;
366 252 : ns_res = ns_resIn;
367 252 : bTargetAlignedPixels = bTargetAlignedPixelsIn;
368 252 : minX = minXIn;
369 252 : minY = minYIn;
370 252 : maxX = maxXIn;
371 252 : maxY = maxYIn;
372 252 : bSeparate = bSeparateIn;
373 252 : bAllowProjectionDifference = bAllowProjectionDifferenceIn;
374 252 : bAddAlpha = bAddAlphaIn;
375 252 : bHideNoData = bHideNoDataIn;
376 252 : nSubdataset = nSubdatasetIn;
377 252 : pszSrcNoData = (pszSrcNoDataIn) ? CPLStrdup(pszSrcNoDataIn) : nullptr;
378 252 : pszVRTNoData = (pszVRTNoDataIn) ? CPLStrdup(pszVRTNoDataIn) : nullptr;
379 252 : pszOutputSRS = (pszOutputSRSIn) ? CPLStrdup(pszOutputSRSIn) : nullptr;
380 252 : pszResampling = (pszResamplingIn) ? CPLStrdup(pszResamplingIn) : nullptr;
381 252 : bUseSrcMaskBand = bUseSrcMaskBandIn;
382 252 : bNoDataFromMask = bNoDataFromMaskIn;
383 252 : dfMaskValueThreshold = dfMaskValueThresholdIn;
384 252 : }
385 :
386 : /************************************************************************/
387 : /* ~VRTBuilder() */
388 : /************************************************************************/
389 :
390 252 : VRTBuilder::~VRTBuilder()
391 : {
392 252 : CPLFree(pszOutputFilename);
393 252 : CPLFree(pszSrcNoData);
394 252 : CPLFree(pszVRTNoData);
395 252 : CPLFree(panSelectedBandList);
396 :
397 252 : if (ppszInputFilenames)
398 : {
399 2658 : for (int i = 0; i < nInputFiles; i++)
400 : {
401 2406 : CPLFree(ppszInputFilenames[i]);
402 : }
403 : }
404 252 : CPLFree(ppszInputFilenames);
405 252 : CPLFree(pahSrcDS);
406 :
407 252 : CPLFree(pszProjectionRef);
408 252 : CPLFree(padfSrcNoData);
409 252 : CPLFree(padfVRTNoData);
410 252 : CPLFree(pszOutputSRS);
411 252 : CPLFree(pszResampling);
412 252 : CSLDestroy(papszOpenOptions);
413 252 : }
414 :
415 : /************************************************************************/
416 : /* ProjAreEqual() */
417 : /************************************************************************/
418 :
419 2152 : static int ProjAreEqual(const char *pszWKT1, const char *pszWKT2)
420 : {
421 2152 : if (EQUAL(pszWKT1, pszWKT2))
422 2150 : return TRUE;
423 :
424 2 : OGRSpatialReferenceH hSRS1 = OSRNewSpatialReference(pszWKT1);
425 2 : OGRSpatialReferenceH hSRS2 = OSRNewSpatialReference(pszWKT2);
426 2 : int bRet = hSRS1 != nullptr && hSRS2 != nullptr && OSRIsSame(hSRS1, hSRS2);
427 2 : if (hSRS1)
428 2 : OSRDestroySpatialReference(hSRS1);
429 2 : if (hSRS2)
430 2 : OSRDestroySpatialReference(hSRS2);
431 2 : return bRet;
432 : }
433 :
434 : /************************************************************************/
435 : /* GetProjectionName() */
436 : /************************************************************************/
437 :
438 4 : static CPLString GetProjectionName(const char *pszProjection)
439 : {
440 4 : if (!pszProjection)
441 0 : return "(null)";
442 :
443 8 : OGRSpatialReference oSRS;
444 4 : oSRS.SetFromUserInput(pszProjection);
445 :
446 4 : const char *pszName = oSRS.GetName();
447 4 : return pszName ? pszName : "(null)";
448 : }
449 :
450 : /************************************************************************/
451 : /* checkNoDataValues() */
452 : /************************************************************************/
453 :
454 2393 : static void checkNoDataValues(const std::vector<BandProperty> &asProperties)
455 : {
456 6098 : for (const auto &oProps : asProperties)
457 : {
458 3763 : if (oProps.bHasNoData && GDALDataTypeIsInteger(oProps.dataType) &&
459 58 : !GDALIsValueExactAs(oProps.noDataValue, oProps.dataType))
460 : {
461 2 : CPLError(CE_Warning, CPLE_NotSupported,
462 : "Band data type of %s cannot represent the specified "
463 : "NoData value of %g",
464 2 : GDALGetDataTypeName(oProps.dataType), oProps.noDataValue);
465 : }
466 : }
467 2393 : }
468 :
469 : /************************************************************************/
470 : /* AnalyseRaster() */
471 : /************************************************************************/
472 :
473 2404 : std::string VRTBuilder::AnalyseRaster(GDALDatasetH hDS,
474 : DatasetProperty *psDatasetProperties)
475 : {
476 2404 : GDALDataset *poDS = GDALDataset::FromHandle(hDS);
477 2404 : const char *dsFileName = poDS->GetDescription();
478 2404 : CSLConstList papszMetadata = poDS->GetMetadata(GDAL_MDD_SUBDATASETS);
479 2404 : if (CSLCount(papszMetadata) > 0 && poDS->GetRasterCount() == 0)
480 : {
481 0 : ppszInputFilenames = static_cast<char **>(CPLRealloc(
482 0 : ppszInputFilenames,
483 0 : sizeof(char *) * (nInputFiles + CSLCount(papszMetadata))));
484 0 : if (nSubdataset < 0)
485 : {
486 0 : int count = 1;
487 : char subdatasetNameKey[80];
488 0 : snprintf(subdatasetNameKey, sizeof(subdatasetNameKey),
489 : "SUBDATASET_%d_NAME", count);
490 0 : while (*papszMetadata != nullptr)
491 : {
492 0 : if (EQUALN(*papszMetadata, subdatasetNameKey,
493 : strlen(subdatasetNameKey)))
494 : {
495 0 : asDatasetProperties.resize(nInputFiles + 1);
496 0 : ppszInputFilenames[nInputFiles] = CPLStrdup(
497 0 : *papszMetadata + strlen(subdatasetNameKey) + 1);
498 0 : nInputFiles++;
499 0 : count++;
500 0 : snprintf(subdatasetNameKey, sizeof(subdatasetNameKey),
501 : "SUBDATASET_%d_NAME", count);
502 : }
503 0 : papszMetadata++;
504 : }
505 : }
506 : else
507 : {
508 : char subdatasetNameKey[80];
509 : const char *pszSubdatasetName;
510 :
511 0 : snprintf(subdatasetNameKey, sizeof(subdatasetNameKey),
512 : "SUBDATASET_%d_NAME", nSubdataset);
513 : pszSubdatasetName =
514 0 : CSLFetchNameValue(papszMetadata, subdatasetNameKey);
515 0 : if (pszSubdatasetName)
516 : {
517 0 : asDatasetProperties.resize(nInputFiles + 1);
518 0 : ppszInputFilenames[nInputFiles] = CPLStrdup(pszSubdatasetName);
519 0 : nInputFiles++;
520 : }
521 : }
522 0 : return "SILENTLY_IGNORE";
523 : }
524 :
525 2404 : const char *proj = poDS->GetProjectionRef();
526 2404 : auto > = psDatasetProperties->gt;
527 2404 : int bGotGeoTransform = poDS->GetGeoTransform(gt) == CE_None;
528 2404 : if (bSeparate)
529 : {
530 43 : std::string osProgramName(m_osProgramName);
531 43 : if (osProgramName == "gdalbuildvrt")
532 36 : osProgramName += " -separate";
533 :
534 43 : if (bFirst)
535 : {
536 24 : bHasGeoTransform = bGotGeoTransform;
537 24 : if (!bHasGeoTransform)
538 : {
539 1 : if (bUserExtent)
540 : {
541 0 : CPLError(CE_Warning, CPLE_NotSupported, "%s",
542 0 : ("User extent ignored by " + osProgramName +
543 : "with ungeoreferenced images.")
544 : .c_str());
545 : }
546 1 : if (resolutionStrategy == USER_RESOLUTION)
547 : {
548 0 : CPLError(CE_Warning, CPLE_NotSupported, "%s",
549 0 : ("User resolution ignored by " + osProgramName +
550 : " with ungeoreferenced images.")
551 : .c_str());
552 : }
553 : }
554 : }
555 19 : else if (bHasGeoTransform != bGotGeoTransform)
556 : {
557 : return osProgramName + " cannot stack ungeoreferenced and "
558 0 : "georeferenced images.";
559 : }
560 20 : else if (!bHasGeoTransform && (nRasterXSize != poDS->GetRasterXSize() ||
561 1 : nRasterYSize != poDS->GetRasterYSize()))
562 : {
563 : return osProgramName + " cannot stack ungeoreferenced images "
564 0 : "that have not the same dimensions.";
565 : }
566 : }
567 : else
568 : {
569 2361 : if (!bGotGeoTransform)
570 : {
571 0 : return m_osProgramName + " does not support ungeoreferenced image.";
572 : }
573 2361 : bHasGeoTransform = TRUE;
574 : }
575 :
576 2404 : if (bGotGeoTransform)
577 : {
578 4804 : if (gt[GEOTRSFRM_ROTATION_PARAM1] != 0 ||
579 2402 : gt[GEOTRSFRM_ROTATION_PARAM2] != 0)
580 : {
581 0 : return m_osProgramName +
582 0 : " does not support rotated geo transforms.";
583 : }
584 2402 : if (gt[GEOTRSFRM_NS_RES] >= 0)
585 : {
586 0 : return m_osProgramName +
587 0 : " does not support positive NS resolution.";
588 : }
589 : }
590 :
591 2404 : psDatasetProperties->nRasterXSize = poDS->GetRasterXSize();
592 2404 : psDatasetProperties->nRasterYSize = poDS->GetRasterYSize();
593 2404 : if (bFirst && bSeparate && !bGotGeoTransform)
594 : {
595 1 : nRasterXSize = poDS->GetRasterXSize();
596 1 : nRasterYSize = poDS->GetRasterYSize();
597 : }
598 :
599 2404 : double ds_minX = gt[GEOTRSFRM_TOPLEFT_X];
600 2404 : double ds_maxY = gt[GEOTRSFRM_TOPLEFT_Y];
601 2404 : double ds_maxX = ds_minX + GDALGetRasterXSize(hDS) * gt[GEOTRSFRM_WE_RES];
602 2404 : double ds_minY = ds_maxY + GDALGetRasterYSize(hDS) * gt[GEOTRSFRM_NS_RES];
603 :
604 2404 : int _nBands = GDALGetRasterCount(hDS);
605 2404 : if (_nBands == 0)
606 : {
607 0 : return "Dataset has no bands";
608 : }
609 2413 : if (bNoDataFromMask &&
610 9 : poDS->GetRasterBand(_nBands)->GetColorInterpretation() == GCI_AlphaBand)
611 3 : _nBands--;
612 :
613 2404 : GDALRasterBand *poFirstBand = poDS->GetRasterBand(1);
614 2404 : poFirstBand->GetBlockSize(&psDatasetProperties->nBlockXSize,
615 : &psDatasetProperties->nBlockYSize);
616 :
617 : /* For the -separate case */
618 2404 : psDatasetProperties->aeBandType.resize(_nBands);
619 2404 : psDatasetProperties->aosDescriptions.resize(_nBands);
620 :
621 2404 : psDatasetProperties->adfNoDataValues.resize(_nBands);
622 2404 : psDatasetProperties->abHasNoData.resize(_nBands);
623 :
624 2404 : psDatasetProperties->adfOffset.resize(_nBands);
625 2404 : psDatasetProperties->abHasOffset.resize(_nBands);
626 :
627 2404 : psDatasetProperties->adfScale.resize(_nBands);
628 2404 : psDatasetProperties->abHasScale.resize(_nBands);
629 :
630 2404 : psDatasetProperties->abHasMaskBand.resize(_nBands);
631 :
632 2404 : psDatasetProperties->bHasDatasetMask =
633 2404 : poFirstBand->GetMaskFlags() == GMF_PER_DATASET;
634 2404 : if (psDatasetProperties->bHasDatasetMask)
635 17 : bHasDatasetMask = TRUE;
636 2404 : poFirstBand->GetMaskBand()->GetBlockSize(
637 : &psDatasetProperties->nMaskBlockXSize,
638 : &psDatasetProperties->nMaskBlockYSize);
639 :
640 2404 : psDatasetProperties->bLastBandIsAlpha = false;
641 2404 : if (poDS->GetRasterBand(_nBands)->GetColorInterpretation() == GCI_AlphaBand)
642 14 : psDatasetProperties->bLastBandIsAlpha = true;
643 :
644 : // Collect overview factors. We only handle power-of-two situations for now
645 2404 : const int nOverviews = poFirstBand->GetOverviewCount();
646 2404 : int nExpectedOvFactor = 2;
647 2416 : for (int j = 0; j < nOverviews; j++)
648 : {
649 23 : GDALRasterBand *poOverview = poFirstBand->GetOverview(j);
650 23 : if (!poOverview)
651 0 : continue;
652 23 : if (poOverview->GetXSize() < 128 && poOverview->GetYSize() < 128)
653 : {
654 11 : break;
655 : }
656 :
657 12 : const int nOvFactor = GDALComputeOvFactor(
658 : poOverview->GetXSize(), poFirstBand->GetXSize(),
659 12 : poOverview->GetYSize(), poFirstBand->GetYSize());
660 :
661 12 : if (nOvFactor != nExpectedOvFactor)
662 0 : break;
663 :
664 12 : psDatasetProperties->anOverviewFactors.push_back(nOvFactor);
665 12 : nExpectedOvFactor *= 2;
666 : }
667 :
668 6003 : for (int j = 0; j < _nBands; j++)
669 : {
670 3599 : GDALRasterBand *poBand = poDS->GetRasterBand(j + 1);
671 :
672 3599 : psDatasetProperties->aeBandType[j] = poBand->GetRasterDataType();
673 :
674 : // Only used by separate mode
675 3599 : if (bSeparate)
676 : {
677 53 : psDatasetProperties->aosDescriptions[j] = poBand->GetDescription();
678 : // Add metadata items
679 53 : CSLConstList papszMD(poBand->GetMetadata());
680 4 : for (const auto &[pszKey, pszValue] :
681 55 : cpl::IterateNameValue(papszMD))
682 : {
683 2 : psDatasetProperties->mapBandMetadata[j][pszKey] = pszValue;
684 : }
685 : }
686 :
687 3599 : if (!bSeparate && nSrcNoDataCount > 0)
688 : {
689 4 : psDatasetProperties->abHasNoData[j] = true;
690 4 : if (j < nSrcNoDataCount)
691 4 : psDatasetProperties->adfNoDataValues[j] = padfSrcNoData[j];
692 : else
693 0 : psDatasetProperties->adfNoDataValues[j] =
694 0 : padfSrcNoData[nSrcNoDataCount - 1];
695 : }
696 : else
697 : {
698 3595 : int bHasNoData = false;
699 7190 : psDatasetProperties->adfNoDataValues[j] =
700 3595 : poBand->GetNoDataValue(&bHasNoData);
701 3595 : psDatasetProperties->abHasNoData[j] = bHasNoData != 0;
702 : }
703 :
704 3599 : int bHasOffset = false;
705 3599 : psDatasetProperties->adfOffset[j] = poBand->GetOffset(&bHasOffset);
706 7198 : psDatasetProperties->abHasOffset[j] =
707 3599 : bHasOffset != 0 && psDatasetProperties->adfOffset[j] != 0.0;
708 :
709 3599 : int bHasScale = false;
710 3599 : psDatasetProperties->adfScale[j] = poBand->GetScale(&bHasScale);
711 7198 : psDatasetProperties->abHasScale[j] =
712 3599 : bHasScale != 0 && psDatasetProperties->adfScale[j] != 1.0;
713 :
714 3599 : const int nMaskFlags = poBand->GetMaskFlags();
715 3599 : psDatasetProperties->abHasMaskBand[j] =
716 7138 : (nMaskFlags != GMF_ALL_VALID && nMaskFlags != GMF_NODATA) ||
717 7138 : poBand->GetColorInterpretation() == GCI_AlphaBand;
718 : }
719 :
720 2404 : if (bSeparate)
721 : {
722 49 : for (int j = 0; j < nSelectedBands; j++)
723 : {
724 7 : if (panSelectedBandList[j] > _nBands)
725 : {
726 : return CPLSPrintf("%s has %d bands, but %d is requested",
727 1 : dsFileName, _nBands, panSelectedBandList[j]);
728 : }
729 : }
730 : }
731 :
732 2403 : if (bFirst)
733 : {
734 251 : nTotalBands = _nBands;
735 251 : if (bAddAlpha && psDatasetProperties->bLastBandIsAlpha)
736 : {
737 4 : bLastBandIsAlpha = true;
738 4 : nTotalBands--;
739 : }
740 :
741 251 : if (proj)
742 251 : pszProjectionRef = CPLStrdup(proj);
743 251 : if (!bUserExtent)
744 : {
745 237 : minX = ds_minX;
746 237 : minY = ds_minY;
747 237 : maxX = ds_maxX;
748 237 : maxY = ds_maxY;
749 : }
750 :
751 251 : if (!bSeparate)
752 : {
753 : // if not provided an explicit band list, take the one of the first
754 : // dataset
755 228 : if (nSelectedBands == 0)
756 : {
757 211 : nSelectedBands = nTotalBands;
758 211 : CPLFree(panSelectedBandList);
759 211 : panSelectedBandList =
760 211 : static_cast<int *>(CPLMalloc(nSelectedBands * sizeof(int)));
761 511 : for (int j = 0; j < nSelectedBands; j++)
762 : {
763 300 : panSelectedBandList[j] = j + 1;
764 : }
765 : }
766 757 : for (int j = 0; j < nSelectedBands; j++)
767 : {
768 529 : nMaxSelectedBandNo =
769 529 : std::max(nMaxSelectedBandNo, panSelectedBandList[j]);
770 : }
771 :
772 228 : asBandProperties.resize(nSelectedBands);
773 756 : for (int j = 0; j < nSelectedBands; j++)
774 : {
775 529 : const int nSelBand = panSelectedBandList[j];
776 529 : if (nSelBand <= 0 || nSelBand > nTotalBands)
777 : {
778 1 : return CPLSPrintf("Invalid band number: %d", nSelBand);
779 : }
780 528 : GDALRasterBand *poBand = poDS->GetRasterBand(nSelBand);
781 1056 : asBandProperties[j].colorInterpretation =
782 528 : poBand->GetColorInterpretation();
783 528 : asBandProperties[j].dataType = poBand->GetRasterDataType();
784 528 : asBandProperties[j].osDescription = poBand->GetDescription();
785 : // Add metadata items
786 528 : const CSLConstList aosMD(poBand->GetMetadata());
787 30 : for (const auto &[pszKey, pszValue] :
788 543 : cpl::IterateNameValue(aosMD))
789 : {
790 15 : asBandProperties[j].mapBandMetadata[pszKey] = pszValue;
791 : }
792 :
793 528 : if (asBandProperties[j].colorInterpretation == GCI_PaletteIndex)
794 : {
795 5 : auto colorTable = poBand->GetColorTable();
796 5 : if (colorTable)
797 : {
798 5 : asBandProperties[j].colorTable.reset(
799 : colorTable->Clone());
800 : }
801 : }
802 : else
803 523 : asBandProperties[j].colorTable = nullptr;
804 :
805 528 : if (nVRTNoDataCount > 0)
806 : {
807 26 : asBandProperties[j].bHasNoData = true;
808 26 : if (j < nVRTNoDataCount)
809 20 : asBandProperties[j].noDataValue = padfVRTNoData[j];
810 : else
811 6 : asBandProperties[j].noDataValue =
812 6 : padfVRTNoData[nVRTNoDataCount - 1];
813 : }
814 : else
815 : {
816 502 : int bHasNoData = false;
817 1004 : asBandProperties[j].noDataValue =
818 502 : poBand->GetNoDataValue(&bHasNoData);
819 502 : asBandProperties[j].bHasNoData = bHasNoData != 0;
820 : }
821 :
822 528 : int bHasOffset = false;
823 528 : asBandProperties[j].dfOffset = poBand->GetOffset(&bHasOffset);
824 528 : asBandProperties[j].bHasOffset =
825 528 : bHasOffset != 0 && asBandProperties[j].dfOffset != 0.0;
826 :
827 528 : int bHasScale = false;
828 528 : asBandProperties[j].dfScale = poBand->GetScale(&bHasScale);
829 528 : asBandProperties[j].bHasScale =
830 528 : bHasScale != 0 && asBandProperties[j].dfScale != 1.0;
831 : }
832 : }
833 : }
834 : else
835 : {
836 2152 : if ((proj != nullptr && pszProjectionRef == nullptr) ||
837 6456 : (proj == nullptr && pszProjectionRef != nullptr) ||
838 2152 : (proj != nullptr && pszProjectionRef != nullptr &&
839 2152 : ProjAreEqual(proj, pszProjectionRef) == FALSE))
840 : {
841 2 : if (!bAllowProjectionDifference)
842 : {
843 4 : CPLString osExpected = GetProjectionName(pszProjectionRef);
844 4 : CPLString osGot = GetProjectionName(proj);
845 2 : return m_osProgramName +
846 : CPLSPrintf(" does not support heterogeneous "
847 : "projection: expected \"%s\", got \"%s\".",
848 2 : osExpected.c_str(), osGot.c_str());
849 : }
850 : }
851 2150 : if (!bSeparate)
852 : {
853 2131 : if (!bExplicitBandList && _nBands != nTotalBands)
854 : {
855 9 : if (bAddAlpha && _nBands == nTotalBands + 1 &&
856 4 : psDatasetProperties->bLastBandIsAlpha)
857 : {
858 4 : bLastBandIsAlpha = true;
859 : }
860 : else
861 : {
862 5 : return m_osProgramName +
863 : CPLSPrintf(" does not support heterogeneous band "
864 : "numbers: expected %d, got %d.",
865 5 : nTotalBands, _nBands);
866 : }
867 : }
868 2122 : else if (bExplicitBandList && _nBands < nMaxSelectedBandNo)
869 : {
870 0 : return m_osProgramName +
871 : CPLSPrintf(" does not support heterogeneous band "
872 : "numbers: expected at least %d, got %d.",
873 0 : nMaxSelectedBandNo, _nBands);
874 : }
875 :
876 5305 : for (int j = 0; j < nSelectedBands; j++)
877 : {
878 3179 : const int nSelBand = panSelectedBandList[j];
879 3179 : CPLAssert(nSelBand >= 1 && nSelBand <= _nBands);
880 3179 : GDALRasterBand *poBand = poDS->GetRasterBand(nSelBand);
881 : // In normal mode we only preserve description if identical across
882 3179 : if (asBandProperties[j].osDescription !=
883 3179 : poBand->GetDescription())
884 : {
885 1 : asBandProperties[j].osDescription = "";
886 : }
887 : // same for metadata
888 3179 : const CPLStringList aosMD(poBand->GetMetadata());
889 3179 : std::vector<std::string> keysToErase;
890 22 : for (const auto &[pszKey, pszValue] :
891 3201 : cpl::IterateNameValue(aosMD))
892 : {
893 : const auto &existingValue =
894 11 : asBandProperties[j].mapBandMetadata[pszKey];
895 19 : if (existingValue.empty() ||
896 8 : !EQUAL(existingValue.c_str(), pszValue))
897 : {
898 8 : keysToErase.push_back(pszKey);
899 : }
900 : }
901 : // Also expand keysToErase to those that are not in the current band
902 3192 : for (const auto &pair : asBandProperties[j].mapBandMetadata)
903 : {
904 13 : if (aosMD.FetchNameValue(pair.first.c_str()) == nullptr)
905 : {
906 2 : keysToErase.push_back(pair.first);
907 : }
908 : }
909 :
910 3189 : for (const auto &key : keysToErase)
911 : {
912 10 : asBandProperties[j].mapBandMetadata.erase(key);
913 : }
914 :
915 3179 : if (asBandProperties[j].colorInterpretation !=
916 3179 : poBand->GetColorInterpretation())
917 : {
918 0 : return m_osProgramName +
919 : CPLSPrintf(
920 : " does not support heterogeneous "
921 : "band color interpretation: expected %s, got "
922 : "%s.",
923 : GDALGetColorInterpretationName(
924 0 : asBandProperties[j].colorInterpretation),
925 : GDALGetColorInterpretationName(
926 0 : poBand->GetColorInterpretation()));
927 : }
928 3179 : if (asBandProperties[j].dataType != poBand->GetRasterDataType())
929 : {
930 0 : return m_osProgramName +
931 : CPLSPrintf(" does not support heterogeneous "
932 : "band data type: expected %s, got %s.",
933 : GDALGetDataTypeName(
934 0 : asBandProperties[j].dataType),
935 : GDALGetDataTypeName(
936 0 : poBand->GetRasterDataType()));
937 : }
938 3179 : if (asBandProperties[j].colorTable)
939 : {
940 4 : const GDALColorTable *colorTable = poBand->GetColorTable();
941 : int nRefColorEntryCount =
942 4 : asBandProperties[j].colorTable->GetColorEntryCount();
943 8 : if (colorTable == nullptr ||
944 4 : (colorTable->GetColorEntryCount() !=
945 2 : nRefColorEntryCount &&
946 : // For NITF CADRG tiles that may have 256 or 257 colors,
947 : // with the 257th color when present being transparent
948 2 : !(std::min(colorTable->GetColorEntryCount(),
949 2 : nRefColorEntryCount) +
950 : 1 ==
951 8 : std::max(colorTable->GetColorEntryCount(),
952 2 : nRefColorEntryCount))))
953 :
954 : {
955 0 : return m_osProgramName +
956 : " does not support rasters with "
957 : "different color tables (different number of "
958 0 : "color table entries)";
959 : }
960 :
961 : /* Check that the palette are the same too */
962 : /* We just warn and still process the file. It is not a
963 : * technical no-go, but the user */
964 : /* should check that the end result is OK for him. */
965 694 : for (int i = 0;
966 694 : i < std::min(nRefColorEntryCount,
967 694 : colorTable->GetColorEntryCount());
968 : i++)
969 : {
970 : const GDALColorEntry *psEntry =
971 690 : colorTable->GetColorEntry(i);
972 : const GDALColorEntry *psEntryRef =
973 690 : asBandProperties[j].colorTable->GetColorEntry(i);
974 690 : if (psEntry->c1 != psEntryRef->c1 ||
975 690 : psEntry->c2 != psEntryRef->c2 ||
976 690 : psEntry->c3 != psEntryRef->c3 ||
977 690 : psEntry->c4 != psEntryRef->c4)
978 : {
979 : static int bFirstWarningPCT = TRUE;
980 0 : if (bFirstWarningPCT)
981 0 : CPLError(
982 : CE_Warning, CPLE_NotSupported,
983 : "%s has different values than the first "
984 : "raster for some entries in the color "
985 : "table.\n"
986 : "The end result might produce weird "
987 : "colors.\n"
988 : "You're advised to pre-process your "
989 : "rasters with other tools, such as "
990 : "pct2rgb.py or gdal_translate -expand RGB\n"
991 : "to operate %s on RGB rasters "
992 : "instead",
993 : dsFileName, m_osProgramName.c_str());
994 : else
995 0 : CPLError(CE_Warning, CPLE_NotSupported,
996 : "%s has different values than the "
997 : "first raster for some entries in the "
998 : "color table.",
999 : dsFileName);
1000 0 : bFirstWarningPCT = FALSE;
1001 0 : break;
1002 : }
1003 : }
1004 :
1005 : // For NITF CADRG tiles that may have 256 or 257 colors,
1006 : // with the 257th color when present being transparent
1007 8 : if (nRefColorEntryCount + 1 ==
1008 4 : colorTable->GetColorEntryCount() &&
1009 1 : colorTable->GetColorEntry(nRefColorEntryCount)->c1 ==
1010 1 : 0 &&
1011 1 : colorTable->GetColorEntry(nRefColorEntryCount)->c2 ==
1012 1 : 0 &&
1013 1 : colorTable->GetColorEntry(nRefColorEntryCount)->c3 ==
1014 5 : 0 &&
1015 1 : colorTable->GetColorEntry(nRefColorEntryCount)->c4 == 0)
1016 : {
1017 1 : int bHasNoData = false;
1018 : const double dfNoData =
1019 1 : poBand->GetNoDataValue(&bHasNoData);
1020 2 : if (bHasNoData && !asBandProperties[j].bHasNoData &&
1021 1 : dfNoData == nRefColorEntryCount)
1022 : {
1023 1 : asBandProperties[j].bHasNoData = true;
1024 1 : asBandProperties[j].noDataValue = dfNoData;
1025 : }
1026 :
1027 1 : asBandProperties[j].colorTable.reset(
1028 : colorTable->Clone());
1029 : }
1030 : }
1031 :
1032 3179 : if (psDatasetProperties->abHasOffset[j] !=
1033 6358 : asBandProperties[j].bHasOffset ||
1034 3179 : (asBandProperties[j].bHasOffset &&
1035 0 : psDatasetProperties->adfOffset[j] !=
1036 0 : asBandProperties[j].dfOffset))
1037 : {
1038 0 : return m_osProgramName +
1039 : CPLSPrintf(
1040 : " does not support heterogeneous "
1041 : "band offset: expected (%d,%f), got (%d,%f).",
1042 0 : static_cast<int>(asBandProperties[j].bHasOffset),
1043 0 : asBandProperties[j].dfOffset,
1044 : static_cast<int>(
1045 0 : psDatasetProperties->abHasOffset[j]),
1046 0 : psDatasetProperties->adfOffset[j]);
1047 : }
1048 :
1049 3179 : if (psDatasetProperties->abHasScale[j] !=
1050 6358 : asBandProperties[j].bHasScale ||
1051 3179 : (asBandProperties[j].bHasScale &&
1052 0 : psDatasetProperties->adfScale[j] !=
1053 0 : asBandProperties[j].dfScale))
1054 : {
1055 0 : return m_osProgramName +
1056 : CPLSPrintf(
1057 : " does not support heterogeneous "
1058 : "band scale: expected (%d,%f), got (%d,%f).",
1059 0 : static_cast<int>(asBandProperties[j].bHasScale),
1060 0 : asBandProperties[j].dfScale,
1061 : static_cast<int>(
1062 0 : psDatasetProperties->abHasScale[j]),
1063 0 : psDatasetProperties->adfScale[j]);
1064 : }
1065 : }
1066 : }
1067 2145 : if (!bUserExtent)
1068 : {
1069 2141 : if (ds_minX < minX)
1070 11 : minX = ds_minX;
1071 2141 : if (ds_minY < minY)
1072 31 : minY = ds_minY;
1073 2141 : if (ds_maxX > maxX)
1074 699 : maxX = ds_maxX;
1075 2141 : if (ds_maxY > maxY)
1076 29 : maxY = ds_maxY;
1077 : }
1078 : }
1079 :
1080 2395 : if (resolutionStrategy == AVERAGE_RESOLUTION)
1081 : {
1082 2314 : ++nCountValid;
1083 : {
1084 2314 : const double dfDelta = gt[GEOTRSFRM_WE_RES] - we_res;
1085 2314 : we_res += dfDelta / nCountValid;
1086 : }
1087 : {
1088 2314 : const double dfDelta = gt[GEOTRSFRM_NS_RES] - ns_res;
1089 2314 : ns_res += dfDelta / nCountValid;
1090 : }
1091 : }
1092 81 : else if (resolutionStrategy == SAME_RESOLUTION)
1093 : {
1094 44 : if (bFirst)
1095 : {
1096 38 : we_res = gt[GEOTRSFRM_WE_RES];
1097 38 : ns_res = gt[GEOTRSFRM_NS_RES];
1098 : }
1099 11 : else if (we_res != gt[GEOTRSFRM_WE_RES] ||
1100 5 : ns_res != gt[GEOTRSFRM_NS_RES])
1101 : {
1102 : return CPLSPrintf(
1103 : "Dataset %s has resolution %.17g x %.17g, whereas "
1104 : "previous sources have resolution %.17g x %.17g. To mosaic "
1105 : "these data sources, a different resolution strategy should be "
1106 : "specified.",
1107 2 : dsFileName, gt[GEOTRSFRM_WE_RES], gt[GEOTRSFRM_NS_RES], we_res,
1108 1 : ns_res);
1109 : }
1110 : }
1111 37 : else if (resolutionStrategy != USER_RESOLUTION)
1112 : {
1113 24 : if (bFirst)
1114 : {
1115 12 : we_res = gt[GEOTRSFRM_WE_RES];
1116 12 : ns_res = gt[GEOTRSFRM_NS_RES];
1117 : }
1118 12 : else if (resolutionStrategy == HIGHEST_RESOLUTION)
1119 : {
1120 1 : we_res = std::min(we_res, gt[GEOTRSFRM_WE_RES]);
1121 : // ns_res is negative, the highest resolution is the max value.
1122 1 : ns_res = std::max(ns_res, gt[GEOTRSFRM_NS_RES]);
1123 : }
1124 11 : else if (resolutionStrategy == COMMON_RESOLUTION)
1125 : {
1126 10 : we_res = CPLGreatestCommonDivisor(we_res, gt[GEOTRSFRM_WE_RES]);
1127 10 : if (we_res == 0)
1128 : {
1129 1 : return "Failed to get common resolution";
1130 : }
1131 9 : ns_res = CPLGreatestCommonDivisor(ns_res, gt[GEOTRSFRM_NS_RES]);
1132 9 : if (ns_res == 0)
1133 : {
1134 0 : return "Failed to get common resolution";
1135 : }
1136 : }
1137 : else
1138 : {
1139 1 : we_res = std::max(we_res, gt[GEOTRSFRM_WE_RES]);
1140 : // ns_res is negative, the lowest resolution is the min value.
1141 1 : ns_res = std::min(ns_res, gt[GEOTRSFRM_NS_RES]);
1142 : }
1143 : }
1144 :
1145 2393 : checkNoDataValues(asBandProperties);
1146 :
1147 2393 : return "";
1148 : }
1149 :
1150 : /************************************************************************/
1151 : /* WriteAbsolutePath() */
1152 : /************************************************************************/
1153 :
1154 7 : static void WriteAbsolutePath(VRTSimpleSource *poSource, const char *dsFileName)
1155 : {
1156 7 : if (dsFileName[0])
1157 : {
1158 7 : if (CPLIsFilenameRelative(dsFileName))
1159 : {
1160 : VSIStatBufL sStat;
1161 2 : if (VSIStatL(dsFileName, &sStat) == 0)
1162 : {
1163 2 : if (char *pszCurDir = CPLGetCurrentDir())
1164 : {
1165 2 : poSource->SetSourceDatasetName(
1166 4 : CPLFormFilenameSafe(pszCurDir, dsFileName, nullptr)
1167 : .c_str(),
1168 : false);
1169 2 : CPLFree(pszCurDir);
1170 : }
1171 : }
1172 : }
1173 : else
1174 : {
1175 5 : poSource->SetSourceDatasetName(dsFileName, false);
1176 : }
1177 : }
1178 7 : }
1179 :
1180 : /************************************************************************/
1181 : /* IsTransientSrcDataset() */
1182 : /************************************************************************/
1183 :
1184 1137 : static bool IsTransientSrcDataset(const char *dsFileName, GDALDatasetH hDS)
1185 : {
1186 1137 : auto hDriver = GDALGetDatasetDriver(hDS);
1187 1137 : return !hDriver || dsFileName[0] == '\0' || // could be a unnamed VRT file
1188 : // Inner pipeline
1189 56 : (dsFileName[0] == '[' &&
1190 2274 : dsFileName[strlen(dsFileName) - 1] == ']') ||
1191 1192 : EQUAL(GDALGetDescription(hDriver), "MEM");
1192 : }
1193 :
1194 : /************************************************************************/
1195 : /* CreateVRTSeparate() */
1196 : /************************************************************************/
1197 :
1198 23 : void VRTBuilder::CreateVRTSeparate(VRTDataset *poVRTDS)
1199 : {
1200 23 : int iBand = 1;
1201 65 : for (int i = 0; ppszInputFilenames != nullptr && i < nInputFiles; i++)
1202 : {
1203 42 : DatasetProperty *psDatasetProperties = &asDatasetProperties[i];
1204 :
1205 42 : if (psDatasetProperties->isFileOK == FALSE)
1206 0 : continue;
1207 :
1208 42 : const char *dsFileName = ppszInputFilenames[i];
1209 :
1210 : double dfSrcXOff, dfSrcYOff, dfSrcXSize, dfSrcYSize, dfDstXOff,
1211 : dfDstYOff, dfDstXSize, dfDstYSize;
1212 42 : if (bHasGeoTransform)
1213 : {
1214 40 : if (!GetSrcDstWin(psDatasetProperties, we_res, ns_res, minX, minY,
1215 : maxX, maxY, nRasterXSize, nRasterYSize,
1216 : &dfSrcXOff, &dfSrcYOff, &dfSrcXSize, &dfSrcYSize,
1217 : &dfDstXOff, &dfDstYOff, &dfDstXSize, &dfDstYSize))
1218 : {
1219 0 : CPLDebug("BuildVRT",
1220 : "Skipping %s as not intersecting area of interest",
1221 : dsFileName);
1222 0 : continue;
1223 : }
1224 : }
1225 : else
1226 : {
1227 2 : dfSrcXOff = dfSrcYOff = dfDstXOff = dfDstYOff = 0;
1228 2 : dfSrcXSize = dfDstXSize = nRasterXSize;
1229 2 : dfSrcYSize = dfDstYSize = nRasterYSize;
1230 : }
1231 :
1232 : GDALDatasetH hSourceDS;
1233 42 : bool bDropRef = false;
1234 64 : if (nSrcDSCount == nInputFiles &&
1235 22 : IsTransientSrcDataset(dsFileName, pahSrcDS[i]))
1236 : {
1237 22 : hSourceDS = pahSrcDS[i];
1238 : }
1239 : else
1240 : {
1241 20 : bDropRef = true;
1242 40 : GDALProxyPoolDatasetH hProxyDS = GDALProxyPoolDatasetCreate(
1243 : dsFileName, psDatasetProperties->nRasterXSize,
1244 : psDatasetProperties->nRasterYSize, GA_ReadOnly, TRUE,
1245 20 : pszProjectionRef, psDatasetProperties->gt.data());
1246 20 : hSourceDS = static_cast<GDALDatasetH>(hProxyDS);
1247 : cpl::down_cast<GDALProxyPoolDataset *>(
1248 : GDALDataset::FromHandle(hProxyDS))
1249 20 : ->SetOpenOptions(papszOpenOptions);
1250 :
1251 43 : for (int jBand = 0;
1252 43 : jBand <
1253 43 : static_cast<int>(psDatasetProperties->aeBandType.size());
1254 : ++jBand)
1255 : {
1256 23 : GDALProxyPoolDatasetAddSrcBandDescription(
1257 23 : hProxyDS, psDatasetProperties->aeBandType[jBand],
1258 : psDatasetProperties->nBlockXSize,
1259 : psDatasetProperties->nBlockYSize);
1260 : }
1261 : }
1262 :
1263 : const int nBandsToIter =
1264 42 : nSelectedBands > 0
1265 82 : ? nSelectedBands
1266 40 : : static_cast<int>(psDatasetProperties->aeBandType.size());
1267 92 : for (int iBandToIter = 0; iBandToIter < nBandsToIter; ++iBandToIter)
1268 : {
1269 : // 0-based
1270 100 : const int nSrcBandIdx = nSelectedBands > 0
1271 50 : ? panSelectedBandList[iBandToIter] - 1
1272 : : iBandToIter;
1273 50 : assert(nSrcBandIdx >= 0);
1274 50 : const auto eBandType = psDatasetProperties->aeBandType[nSrcBandIdx];
1275 50 : poVRTDS->AddBand(eBandType, nullptr);
1276 :
1277 : VRTSourcedRasterBand *poVRTBand =
1278 : static_cast<VRTSourcedRasterBand *>(
1279 50 : poVRTDS->GetRasterBand(iBand));
1280 :
1281 50 : poVRTBand->SetDescription(
1282 50 : psDatasetProperties->aosDescriptions[nSrcBandIdx].c_str());
1283 50 : if (!psDatasetProperties->mapBandMetadata[nSrcBandIdx].empty())
1284 : {
1285 4 : for (const auto &[key, value] :
1286 6 : psDatasetProperties->mapBandMetadata[nSrcBandIdx])
1287 : {
1288 2 : poVRTBand->SetMetadataItem(key.c_str(), value.c_str());
1289 : }
1290 : }
1291 :
1292 50 : if (bHideNoData)
1293 0 : poVRTBand->SetMetadataItem("HideNoDataValue", "1", nullptr);
1294 :
1295 50 : if (bAllowVRTNoData)
1296 : {
1297 48 : std::optional<double> noData;
1298 48 : if (nVRTNoDataCount > 0)
1299 : {
1300 5 : if (iBand - 1 < nVRTNoDataCount)
1301 5 : noData = padfVRTNoData[iBand - 1];
1302 : else
1303 0 : noData = padfVRTNoData[nVRTNoDataCount - 1];
1304 : }
1305 43 : else if (psDatasetProperties->abHasNoData[nSrcBandIdx])
1306 : {
1307 2 : noData = psDatasetProperties->adfNoDataValues[nSrcBandIdx];
1308 : }
1309 48 : if (noData.has_value())
1310 : {
1311 14 : if (GDALDataTypeIsInteger(eBandType) &&
1312 7 : !GDALIsValueExactAs(*noData, eBandType))
1313 : {
1314 1 : CPLError(CE_Warning, CPLE_NotSupported,
1315 : "Band data type of %s cannot represent the "
1316 : "specified "
1317 : "NoData value of %g",
1318 1 : GDALGetDataTypeName(eBandType), *noData);
1319 : }
1320 7 : poVRTBand->SetNoDataValue(*noData);
1321 : }
1322 : }
1323 :
1324 : VRTSimpleSource *poSimpleSource;
1325 98 : if (bAllowSrcNoData &&
1326 91 : (nSrcNoDataCount > 0 ||
1327 93 : psDatasetProperties->abHasNoData[nSrcBandIdx]))
1328 : {
1329 7 : auto poComplexSource = new VRTComplexSource();
1330 7 : poSimpleSource = poComplexSource;
1331 7 : if (nSrcNoDataCount > 0)
1332 : {
1333 5 : if (iBand - 1 < nSrcNoDataCount)
1334 5 : poComplexSource->SetNoDataValue(
1335 5 : padfSrcNoData[iBand - 1]);
1336 : else
1337 0 : poComplexSource->SetNoDataValue(
1338 0 : padfSrcNoData[nSrcNoDataCount - 1]);
1339 : }
1340 : else /* if (psDatasetProperties->abHasNoData[nSrcBandIdx]) */
1341 : {
1342 2 : poComplexSource->SetNoDataValue(
1343 2 : psDatasetProperties->adfNoDataValues[nSrcBandIdx]);
1344 : }
1345 : }
1346 86 : else if (bUseSrcMaskBand &&
1347 86 : psDatasetProperties->abHasMaskBand[nSrcBandIdx])
1348 : {
1349 1 : auto poSource = new VRTComplexSource();
1350 1 : poSource->SetUseMaskBand(true);
1351 1 : poSimpleSource = poSource;
1352 : }
1353 : else
1354 42 : poSimpleSource = new VRTSimpleSource();
1355 :
1356 50 : if (pszResampling)
1357 0 : poSimpleSource->SetResampling(pszResampling);
1358 100 : poVRTBand->ConfigureSource(
1359 : poSimpleSource,
1360 : static_cast<GDALRasterBand *>(
1361 50 : GDALGetRasterBand(hSourceDS, nSrcBandIdx + 1)),
1362 : FALSE, dfSrcXOff, dfSrcYOff, dfSrcXSize, dfSrcYSize, dfDstXOff,
1363 : dfDstYOff, dfDstXSize, dfDstYSize);
1364 :
1365 50 : if (bWriteAbsolutePath)
1366 3 : WriteAbsolutePath(poSimpleSource, dsFileName);
1367 :
1368 50 : if (psDatasetProperties->abHasOffset[nSrcBandIdx])
1369 0 : poVRTBand->SetOffset(
1370 0 : psDatasetProperties->adfOffset[nSrcBandIdx]);
1371 :
1372 50 : if (psDatasetProperties->abHasScale[nSrcBandIdx])
1373 0 : poVRTBand->SetScale(psDatasetProperties->adfScale[nSrcBandIdx]);
1374 :
1375 50 : poVRTBand->AddSource(poSimpleSource);
1376 :
1377 50 : iBand++;
1378 : }
1379 :
1380 42 : if (bDropRef)
1381 : {
1382 20 : GDALDereferenceDataset(hSourceDS);
1383 : }
1384 : }
1385 23 : }
1386 :
1387 : /************************************************************************/
1388 : /* CreateVRTNonSeparate() */
1389 : /************************************************************************/
1390 :
1391 223 : void VRTBuilder::CreateVRTNonSeparate(VRTDataset *poVRTDS)
1392 : {
1393 446 : CPLStringList aosOptions;
1394 :
1395 223 : if (!osPixelFunction.empty())
1396 : {
1397 3 : aosOptions.AddNameValue("subclass", "VRTDerivedRasterBand");
1398 3 : aosOptions.AddNameValue("PixelFunctionType", osPixelFunction.c_str());
1399 3 : aosOptions.AddNameValue("SkipNonContributingSources", "1");
1400 6 : CPLString osName;
1401 2 : for (const auto &[pszKey, pszValue] :
1402 5 : cpl::IterateNameValue(aosPixelFunctionArgs))
1403 : {
1404 1 : osName.Printf("_PIXELFN_ARG_%s", pszKey);
1405 1 : aosOptions.AddNameValue(osName.c_str(), pszValue);
1406 : }
1407 : }
1408 :
1409 747 : for (int j = 0; j < nSelectedBands; j++)
1410 : {
1411 524 : const char *pszSourceTransferType = "Float64";
1412 1047 : if (osPixelFunction == "mean" || osPixelFunction == "min" ||
1413 523 : osPixelFunction == "max")
1414 : {
1415 : pszSourceTransferType =
1416 1 : GDALGetDataTypeName(asBandProperties[j].dataType);
1417 : }
1418 524 : aosOptions.AddNameValue("SourceTransferType", pszSourceTransferType);
1419 :
1420 524 : poVRTDS->AddBand(asBandProperties[j].dataType, aosOptions.List());
1421 524 : GDALRasterBand *poBand = poVRTDS->GetRasterBand(j + 1);
1422 524 : poBand->SetColorInterpretation(asBandProperties[j].colorInterpretation);
1423 524 : poBand->SetDescription(asBandProperties[j].osDescription.c_str());
1424 524 : if (!asBandProperties[j].mapBandMetadata.empty())
1425 : {
1426 12 : for (const auto &[key, value] : asBandProperties[j].mapBandMetadata)
1427 : {
1428 8 : poBand->SetMetadataItem(key.c_str(), value.c_str());
1429 : }
1430 : }
1431 524 : if (asBandProperties[j].colorInterpretation == GCI_PaletteIndex)
1432 : {
1433 5 : poBand->SetColorTable(asBandProperties[j].colorTable.get());
1434 : }
1435 524 : if (bAllowVRTNoData && asBandProperties[j].bHasNoData)
1436 42 : poBand->SetNoDataValue(asBandProperties[j].noDataValue);
1437 524 : if (bHideNoData)
1438 2 : poBand->SetMetadataItem("HideNoDataValue", "1");
1439 :
1440 524 : if (asBandProperties[j].bHasOffset)
1441 0 : poBand->SetOffset(asBandProperties[j].dfOffset);
1442 :
1443 524 : if (asBandProperties[j].bHasScale)
1444 0 : poBand->SetScale(asBandProperties[j].dfScale);
1445 : }
1446 :
1447 223 : VRTSourcedRasterBand *poMaskVRTBand = nullptr;
1448 223 : if (bAddAlpha)
1449 : {
1450 11 : poVRTDS->AddBand(GDT_UInt8);
1451 11 : GDALRasterBand *poBand = poVRTDS->GetRasterBand(nSelectedBands + 1);
1452 11 : poBand->SetColorInterpretation(GCI_AlphaBand);
1453 : }
1454 212 : else if (bHasDatasetMask)
1455 : {
1456 12 : poVRTDS->CreateMaskBand(GMF_PER_DATASET);
1457 : poMaskVRTBand = static_cast<VRTSourcedRasterBand *>(
1458 12 : poVRTDS->GetRasterBand(1)->GetMaskBand());
1459 : }
1460 :
1461 223 : bool bCanCollectOverviewFactors = true;
1462 446 : std::set<int> anOverviewFactorsSet;
1463 446 : std::vector<int> anIdxValidDatasets;
1464 :
1465 2577 : for (int i = 0; ppszInputFilenames != nullptr && i < nInputFiles; i++)
1466 : {
1467 2354 : DatasetProperty *psDatasetProperties = &asDatasetProperties[i];
1468 :
1469 2354 : if (psDatasetProperties->isFileOK == FALSE)
1470 8 : continue;
1471 :
1472 2347 : const char *dsFileName = ppszInputFilenames[i];
1473 :
1474 : double dfSrcXOff;
1475 : double dfSrcYOff;
1476 : double dfSrcXSize;
1477 : double dfSrcYSize;
1478 : double dfDstXOff;
1479 : double dfDstYOff;
1480 : double dfDstXSize;
1481 : double dfDstYSize;
1482 2347 : if (!GetSrcDstWin(psDatasetProperties, we_res, ns_res, minX, minY, maxX,
1483 : maxY, nRasterXSize, nRasterYSize, &dfSrcXOff,
1484 : &dfSrcYOff, &dfSrcXSize, &dfSrcYSize, &dfDstXOff,
1485 : &dfDstYOff, &dfDstXSize, &dfDstYSize))
1486 : {
1487 1 : CPLDebug("BuildVRT",
1488 : "Skipping %s as not intersecting area of interest",
1489 : dsFileName);
1490 1 : continue;
1491 : }
1492 :
1493 2346 : anIdxValidDatasets.push_back(i);
1494 :
1495 2346 : if (bCanCollectOverviewFactors)
1496 : {
1497 2331 : if (std::abs(psDatasetProperties->gt.xscale - we_res) >
1498 4643 : 1e-8 * std::abs(we_res) ||
1499 2312 : std::abs(psDatasetProperties->gt.yscale - ns_res) >
1500 2312 : 1e-8 * std::abs(ns_res))
1501 : {
1502 19 : bCanCollectOverviewFactors = false;
1503 19 : anOverviewFactorsSet.clear();
1504 : }
1505 : }
1506 2346 : if (bCanCollectOverviewFactors)
1507 : {
1508 2321 : for (int nOvFactor : psDatasetProperties->anOverviewFactors)
1509 9 : anOverviewFactorsSet.insert(nOvFactor);
1510 : }
1511 :
1512 : GDALDatasetH hSourceDS;
1513 2346 : bool bDropRef = false;
1514 :
1515 3461 : if (nSrcDSCount == nInputFiles &&
1516 1115 : IsTransientSrcDataset(dsFileName, pahSrcDS[i]))
1517 : {
1518 1096 : hSourceDS = pahSrcDS[i];
1519 : }
1520 : else
1521 : {
1522 1250 : bDropRef = true;
1523 2500 : GDALProxyPoolDatasetH hProxyDS = GDALProxyPoolDatasetCreate(
1524 : dsFileName, psDatasetProperties->nRasterXSize,
1525 : psDatasetProperties->nRasterYSize, GA_ReadOnly, TRUE,
1526 1250 : pszProjectionRef, psDatasetProperties->gt.data());
1527 : cpl::down_cast<GDALProxyPoolDataset *>(
1528 : GDALDataset::FromHandle(hProxyDS))
1529 1250 : ->SetOpenOptions(papszOpenOptions);
1530 :
1531 3614 : for (int j = 0;
1532 3614 : j < nMaxSelectedBandNo +
1533 42 : (bAddAlpha && psDatasetProperties->bLastBandIsAlpha
1534 3656 : ? 1
1535 : : 0);
1536 : j++)
1537 : {
1538 2364 : GDALProxyPoolDatasetAddSrcBandDescription(
1539 : hProxyDS,
1540 2364 : j < static_cast<int>(asBandProperties.size())
1541 2359 : ? asBandProperties[j].dataType
1542 : : GDT_UInt8,
1543 : psDatasetProperties->nBlockXSize,
1544 : psDatasetProperties->nBlockYSize);
1545 : }
1546 1250 : if (bHasDatasetMask && !bAddAlpha)
1547 : {
1548 : static_cast<GDALProxyPoolRasterBand *>(
1549 13 : cpl::down_cast<GDALProxyPoolDataset *>(
1550 : GDALDataset::FromHandle(hProxyDS))
1551 13 : ->GetRasterBand(1))
1552 13 : ->AddSrcMaskBandDescription(
1553 : GDT_UInt8, psDatasetProperties->nMaskBlockXSize,
1554 : psDatasetProperties->nMaskBlockYSize);
1555 : }
1556 :
1557 1250 : hSourceDS = static_cast<GDALDatasetH>(hProxyDS);
1558 : }
1559 :
1560 6054 : for (int j = 0;
1561 6054 : j <
1562 6054 : nSelectedBands +
1563 6054 : (bAddAlpha && psDatasetProperties->bLastBandIsAlpha ? 1 : 0);
1564 : j++)
1565 : {
1566 : VRTSourcedRasterBandH hVRTBand = static_cast<VRTSourcedRasterBandH>(
1567 3708 : poVRTDS->GetRasterBand(j + 1));
1568 3708 : const int nSelBand = j == nSelectedBands ? nSelectedBands + 1
1569 3700 : : panSelectedBandList[j];
1570 :
1571 : /* Place the raster band at the right position in the VRT */
1572 3708 : VRTSourcedRasterBand *poVRTBand =
1573 : static_cast<VRTSourcedRasterBand *>(hVRTBand);
1574 :
1575 : VRTSimpleSource *poSimpleSource;
1576 3708 : if (bNoDataFromMask)
1577 : {
1578 15 : auto poNoDataFromMaskSource = new VRTNoDataFromMaskSource();
1579 15 : poSimpleSource = poNoDataFromMaskSource;
1580 15 : poNoDataFromMaskSource->SetParameters(
1581 15 : (nVRTNoDataCount > 0)
1582 15 : ? ((j < nVRTNoDataCount)
1583 15 : ? padfVRTNoData[j]
1584 6 : : padfVRTNoData[nVRTNoDataCount - 1])
1585 : : 0,
1586 : dfMaskValueThreshold);
1587 : }
1588 7386 : else if (bAllowSrcNoData &&
1589 7386 : psDatasetProperties->abHasNoData[nSelBand - 1])
1590 : {
1591 30 : auto poComplexSource = new VRTComplexSource();
1592 30 : poSimpleSource = poComplexSource;
1593 30 : poComplexSource->SetNoDataValue(
1594 30 : psDatasetProperties->adfNoDataValues[nSelBand - 1]);
1595 : }
1596 7326 : else if (bUseSrcMaskBand &&
1597 7326 : psDatasetProperties->abHasMaskBand[nSelBand - 1])
1598 : {
1599 58 : auto poSource = new VRTComplexSource();
1600 58 : poSource->SetUseMaskBand(true);
1601 58 : poSimpleSource = poSource;
1602 : }
1603 : else
1604 3605 : poSimpleSource = new VRTSimpleSource();
1605 3708 : if (pszResampling)
1606 23 : poSimpleSource->SetResampling(pszResampling);
1607 3708 : auto poSrcBand = GDALRasterBand::FromHandle(
1608 : GDALGetRasterBand(hSourceDS, nSelBand));
1609 3708 : poVRTBand->ConfigureSource(poSimpleSource, poSrcBand, FALSE,
1610 : dfSrcXOff, dfSrcYOff, dfSrcXSize,
1611 : dfSrcYSize, dfDstXOff, dfDstYOff,
1612 : dfDstXSize, dfDstYSize);
1613 :
1614 3708 : if (bWriteAbsolutePath)
1615 3 : WriteAbsolutePath(poSimpleSource, dsFileName);
1616 :
1617 3708 : poVRTBand->AddSource(poSimpleSource);
1618 : }
1619 :
1620 2346 : if (bAddAlpha && !psDatasetProperties->bLastBandIsAlpha)
1621 : {
1622 : VRTSourcedRasterBand *poVRTBand =
1623 : static_cast<VRTSourcedRasterBand *>(
1624 12 : poVRTDS->GetRasterBand(nSelectedBands + 1));
1625 12 : if (psDatasetProperties->bHasDatasetMask && bUseSrcMaskBand)
1626 : {
1627 1 : auto poComplexSource = new VRTComplexSource();
1628 1 : poComplexSource->SetUseMaskBand(true);
1629 1 : poVRTBand->ConfigureSource(
1630 : poComplexSource,
1631 : GDALRasterBand::FromHandle(GDALGetRasterBand(hSourceDS, 1)),
1632 : TRUE, dfSrcXOff, dfSrcYOff, dfSrcXSize, dfSrcYSize,
1633 : dfDstXOff, dfDstYOff, dfDstXSize, dfDstYSize);
1634 :
1635 1 : if (bWriteAbsolutePath)
1636 0 : WriteAbsolutePath(poComplexSource, dsFileName);
1637 :
1638 1 : poVRTBand->AddSource(poComplexSource);
1639 : }
1640 : else
1641 : {
1642 : /* Little trick : we use an offset of 255 and a scaling of 0, so
1643 : * that in areas covered */
1644 : /* by the source, the value of the alpha band will be 255, otherwise
1645 : * it will be 0 */
1646 11 : poVRTBand->AddComplexSource(
1647 : GDALRasterBand::FromHandle(GDALGetRasterBand(hSourceDS, 1)),
1648 : dfSrcXOff, dfSrcYOff, dfSrcXSize, dfSrcYSize, dfDstXOff,
1649 : dfDstYOff, dfDstXSize, dfDstYSize, 255, 0,
1650 : VRT_NODATA_UNSET);
1651 12 : }
1652 : }
1653 2334 : else if (bHasDatasetMask)
1654 : {
1655 : VRTSimpleSource *poSource;
1656 15 : if (bUseSrcMaskBand)
1657 : {
1658 15 : auto poComplexSource = new VRTComplexSource();
1659 15 : poComplexSource->SetUseMaskBand(true);
1660 15 : poSource = poComplexSource;
1661 : }
1662 : else
1663 : {
1664 0 : poSource = new VRTSimpleSource();
1665 : }
1666 15 : if (pszResampling)
1667 4 : poSource->SetResampling(pszResampling);
1668 15 : assert(poMaskVRTBand);
1669 15 : poMaskVRTBand->ConfigureSource(
1670 : poSource,
1671 : GDALRasterBand::FromHandle(GDALGetRasterBand(hSourceDS, 1)),
1672 : TRUE, dfSrcXOff, dfSrcYOff, dfSrcXSize, dfSrcYSize, dfDstXOff,
1673 : dfDstYOff, dfDstXSize, dfDstYSize);
1674 :
1675 15 : if (bWriteAbsolutePath)
1676 1 : WriteAbsolutePath(poSource, dsFileName);
1677 :
1678 15 : poMaskVRTBand->AddSource(poSource);
1679 : }
1680 :
1681 2346 : if (bDropRef)
1682 : {
1683 1250 : GDALDereferenceDataset(hSourceDS);
1684 : }
1685 : }
1686 :
1687 2569 : for (int i : anIdxValidDatasets)
1688 : {
1689 2346 : const DatasetProperty *psDatasetProperties = &asDatasetProperties[i];
1690 2346 : for (auto oIter = anOverviewFactorsSet.begin();
1691 2355 : oIter != anOverviewFactorsSet.end();)
1692 : {
1693 9 : const int nGlobalOvrFactor = *oIter;
1694 9 : auto oIterNext = oIter;
1695 9 : ++oIterNext;
1696 :
1697 9 : if (psDatasetProperties->nRasterXSize / nGlobalOvrFactor < 128 &&
1698 0 : psDatasetProperties->nRasterYSize / nGlobalOvrFactor < 128)
1699 : {
1700 0 : break;
1701 : }
1702 9 : if (std::find(psDatasetProperties->anOverviewFactors.begin(),
1703 : psDatasetProperties->anOverviewFactors.end(),
1704 9 : nGlobalOvrFactor) ==
1705 18 : psDatasetProperties->anOverviewFactors.end())
1706 : {
1707 0 : anOverviewFactorsSet.erase(oIter);
1708 : }
1709 :
1710 9 : oIter = oIterNext;
1711 : }
1712 : }
1713 226 : if (!anOverviewFactorsSet.empty() &&
1714 3 : CPLTestBool(CPLGetConfigOption("VRT_VIRTUAL_OVERVIEWS", "YES")))
1715 : {
1716 6 : std::vector<int> anOverviewFactors;
1717 3 : anOverviewFactors.insert(anOverviewFactors.end(),
1718 : anOverviewFactorsSet.begin(),
1719 6 : anOverviewFactorsSet.end());
1720 3 : const char *const apszOptions[] = {"VRT_VIRTUAL_OVERVIEWS=YES",
1721 : nullptr};
1722 3 : poVRTDS->BuildOverviews(pszResampling ? pszResampling : "nearest",
1723 3 : static_cast<int>(anOverviewFactors.size()),
1724 3 : &anOverviewFactors[0], 0, nullptr, nullptr,
1725 : nullptr, apszOptions);
1726 : }
1727 223 : }
1728 :
1729 : /************************************************************************/
1730 : /* Build() */
1731 : /************************************************************************/
1732 :
1733 252 : std::unique_ptr<GDALDataset> VRTBuilder::Build(GDALProgressFunc pfnProgress,
1734 : void *pProgressData)
1735 : {
1736 252 : if (bHasRunBuild)
1737 0 : return nullptr;
1738 252 : bHasRunBuild = TRUE;
1739 :
1740 252 : if (pfnProgress == nullptr)
1741 0 : pfnProgress = GDALDummyProgress;
1742 :
1743 252 : bUserExtent = (minX != 0 || minY != 0 || maxX != 0 || maxY != 0);
1744 252 : if (bUserExtent)
1745 : {
1746 14 : if (minX >= maxX || minY >= maxY)
1747 : {
1748 0 : CPLError(CE_Failure, CPLE_IllegalArg, "Invalid user extent");
1749 0 : return nullptr;
1750 : }
1751 : }
1752 :
1753 252 : if (resolutionStrategy == USER_RESOLUTION)
1754 : {
1755 8 : if (we_res <= 0 || ns_res <= 0)
1756 : {
1757 0 : CPLError(CE_Failure, CPLE_IllegalArg, "Invalid user resolution");
1758 0 : return nullptr;
1759 : }
1760 :
1761 : /* We work with negative north-south resolution in all the following
1762 : * code */
1763 8 : ns_res = -ns_res;
1764 : }
1765 : else
1766 : {
1767 244 : we_res = ns_res = 0;
1768 : }
1769 :
1770 252 : asDatasetProperties.resize(nInputFiles);
1771 :
1772 252 : if (pszSrcNoData != nullptr)
1773 : {
1774 7 : if (EQUAL(pszSrcNoData, "none"))
1775 : {
1776 1 : bAllowSrcNoData = FALSE;
1777 : }
1778 : else
1779 : {
1780 6 : char **papszTokens = CSLTokenizeString(pszSrcNoData);
1781 6 : nSrcNoDataCount = CSLCount(papszTokens);
1782 6 : padfSrcNoData = static_cast<double *>(
1783 6 : CPLMalloc(sizeof(double) * nSrcNoDataCount));
1784 14 : for (int i = 0; i < nSrcNoDataCount; i++)
1785 : {
1786 8 : if (!ArgIsNumeric(papszTokens[i]) &&
1787 0 : !EQUAL(papszTokens[i], "nan") &&
1788 8 : !EQUAL(papszTokens[i], "-inf") &&
1789 0 : !EQUAL(papszTokens[i], "inf"))
1790 : {
1791 0 : CPLError(CE_Failure, CPLE_IllegalArg,
1792 : "Invalid -srcnodata value");
1793 0 : CSLDestroy(papszTokens);
1794 0 : return nullptr;
1795 : }
1796 8 : padfSrcNoData[i] = CPLAtofM(papszTokens[i]);
1797 : }
1798 6 : CSLDestroy(papszTokens);
1799 : }
1800 : }
1801 :
1802 252 : if (pszVRTNoData != nullptr)
1803 : {
1804 24 : if (EQUAL(pszVRTNoData, "none"))
1805 : {
1806 1 : bAllowVRTNoData = FALSE;
1807 : }
1808 : else
1809 : {
1810 23 : char **papszTokens = CSLTokenizeString(pszVRTNoData);
1811 23 : nVRTNoDataCount = CSLCount(papszTokens);
1812 23 : padfVRTNoData = static_cast<double *>(
1813 23 : CPLMalloc(sizeof(double) * nVRTNoDataCount));
1814 48 : for (int i = 0; i < nVRTNoDataCount; i++)
1815 : {
1816 25 : if (!ArgIsNumeric(papszTokens[i]) &&
1817 1 : !EQUAL(papszTokens[i], "nan") &&
1818 26 : !EQUAL(papszTokens[i], "-inf") &&
1819 0 : !EQUAL(papszTokens[i], "inf"))
1820 : {
1821 0 : CPLError(CE_Failure, CPLE_IllegalArg,
1822 : "Invalid -vrtnodata value");
1823 0 : CSLDestroy(papszTokens);
1824 0 : return nullptr;
1825 : }
1826 25 : padfVRTNoData[i] = CPLAtofM(papszTokens[i]);
1827 : }
1828 23 : CSLDestroy(papszTokens);
1829 : }
1830 : }
1831 :
1832 252 : bool bFoundValid = false;
1833 2654 : for (int i = 0; ppszInputFilenames != nullptr && i < nInputFiles; i++)
1834 : {
1835 2406 : const char *dsFileName = ppszInputFilenames[i];
1836 :
1837 2406 : if (!pfnProgress(1.0 * (i + 1) / nInputFiles, nullptr, pProgressData))
1838 : {
1839 0 : return nullptr;
1840 : }
1841 :
1842 2406 : GDALDatasetH hDS = (pahSrcDS)
1843 2406 : ? pahSrcDS[i]
1844 1259 : : GDALOpenEx(dsFileName, GDAL_OF_RASTER, nullptr,
1845 1259 : papszOpenOptions, nullptr);
1846 2406 : asDatasetProperties[i].isFileOK = FALSE;
1847 :
1848 2406 : if (hDS)
1849 : {
1850 2404 : const auto osErrorMsg = AnalyseRaster(hDS, &asDatasetProperties[i]);
1851 2404 : if (osErrorMsg.empty())
1852 : {
1853 2393 : asDatasetProperties[i].isFileOK = TRUE;
1854 2393 : bFoundValid = true;
1855 2393 : bFirst = FALSE;
1856 : }
1857 2404 : if (pahSrcDS == nullptr)
1858 1257 : GDALClose(hDS);
1859 2404 : if (!osErrorMsg.empty() && osErrorMsg != "SILENTLY_IGNORE")
1860 : {
1861 11 : if (bStrict)
1862 : {
1863 3 : CPLError(CE_Failure, CPLE_AppDefined, "%s",
1864 : osErrorMsg.c_str());
1865 3 : return nullptr;
1866 : }
1867 : else
1868 : {
1869 8 : CPLError(CE_Warning, CPLE_AppDefined, "%s Skipping %s",
1870 : osErrorMsg.c_str(), dsFileName);
1871 : }
1872 : }
1873 : }
1874 : else
1875 : {
1876 2 : if (bStrict)
1877 : {
1878 1 : CPLError(CE_Failure, CPLE_AppDefined, "Can't open %s.",
1879 : dsFileName);
1880 1 : return nullptr;
1881 : }
1882 : else
1883 : {
1884 1 : CPLError(CE_Warning, CPLE_AppDefined,
1885 : "Can't open %s. Skipping it", dsFileName);
1886 : }
1887 : }
1888 : }
1889 :
1890 248 : if (!bFoundValid)
1891 2 : return nullptr;
1892 :
1893 246 : if (bHasGeoTransform)
1894 : {
1895 245 : if (bTargetAlignedPixels)
1896 : {
1897 2 : minX = floor(minX / we_res) * we_res;
1898 2 : maxX = ceil(maxX / we_res) * we_res;
1899 2 : minY = floor(minY / -ns_res) * -ns_res;
1900 2 : maxY = ceil(maxY / -ns_res) * -ns_res;
1901 : }
1902 :
1903 245 : nRasterXSize = static_cast<int>(0.5 + (maxX - minX) / we_res);
1904 245 : nRasterYSize = static_cast<int>(0.5 + (maxY - minY) / -ns_res);
1905 : }
1906 :
1907 246 : if (nRasterXSize == 0 || nRasterYSize == 0)
1908 : {
1909 0 : CPLError(CE_Failure, CPLE_AppDefined,
1910 : "Computed VRT dimension is invalid. You've probably "
1911 : "specified inappropriate resolution.");
1912 0 : return nullptr;
1913 : }
1914 :
1915 246 : auto poDS = VRTDataset::CreateVRTDataset(pszOutputFilename, nRasterXSize,
1916 : nRasterYSize, 0, GDT_Unknown,
1917 492 : aosCreateOptions.List());
1918 246 : if (!poDS)
1919 : {
1920 0 : return nullptr;
1921 : }
1922 :
1923 246 : if (pszOutputSRS)
1924 : {
1925 1 : poDS->SetProjection(pszOutputSRS);
1926 : }
1927 245 : else if (pszProjectionRef)
1928 : {
1929 245 : poDS->SetProjection(pszProjectionRef);
1930 : }
1931 :
1932 246 : if (bHasGeoTransform)
1933 : {
1934 245 : GDALGeoTransform gt;
1935 245 : gt[GEOTRSFRM_TOPLEFT_X] = minX;
1936 245 : gt[GEOTRSFRM_WE_RES] = we_res;
1937 245 : gt[GEOTRSFRM_ROTATION_PARAM1] = 0;
1938 245 : gt[GEOTRSFRM_TOPLEFT_Y] = maxY;
1939 245 : gt[GEOTRSFRM_ROTATION_PARAM2] = 0;
1940 245 : gt[GEOTRSFRM_NS_RES] = ns_res;
1941 245 : poDS->SetGeoTransform(gt);
1942 : }
1943 :
1944 246 : if (bSeparate)
1945 : {
1946 23 : CreateVRTSeparate(poDS.get());
1947 : }
1948 : else
1949 : {
1950 223 : CreateVRTNonSeparate(poDS.get());
1951 : }
1952 :
1953 246 : return poDS;
1954 : }
1955 :
1956 : /************************************************************************/
1957 : /* add_file_to_list() */
1958 : /************************************************************************/
1959 :
1960 45 : static bool add_file_to_list(const char *filename, const char *tile_index,
1961 : CPLStringList &aosList)
1962 : {
1963 :
1964 45 : if (EQUAL(CPLGetExtensionSafe(filename).c_str(), "SHP"))
1965 : {
1966 : /* Handle gdaltindex Shapefile as a special case */
1967 1 : auto poDS = std::unique_ptr<GDALDataset>(GDALDataset::Open(filename));
1968 1 : if (poDS == nullptr)
1969 : {
1970 0 : CPLError(CE_Failure, CPLE_AppDefined,
1971 : "Unable to open shapefile `%s'.", filename);
1972 0 : return false;
1973 : }
1974 :
1975 1 : auto poLayer = poDS->GetLayer(0);
1976 1 : const auto poFDefn = poLayer->GetLayerDefn();
1977 :
1978 2 : if (poFDefn->GetFieldIndex("LOCATION") >= 0 &&
1979 1 : strcmp("LOCATION", tile_index) != 0)
1980 : {
1981 1 : CPLError(CE_Failure, CPLE_AppDefined,
1982 : "This shapefile seems to be a tile index of "
1983 : "OGR features and not GDAL products.");
1984 : }
1985 1 : const int ti_field = poFDefn->GetFieldIndex(tile_index);
1986 1 : if (ti_field < 0)
1987 : {
1988 0 : CPLError(CE_Failure, CPLE_AppDefined,
1989 : "Unable to find field `%s' in DBF file `%s'.", tile_index,
1990 : filename);
1991 0 : return false;
1992 : }
1993 :
1994 : /* Load in memory existing file names in SHP */
1995 1 : const auto nTileIndexFiles = poLayer->GetFeatureCount(TRUE);
1996 1 : if (nTileIndexFiles == 0)
1997 : {
1998 0 : CPLError(CE_Warning, CPLE_AppDefined,
1999 : "Tile index %s is empty. Skipping it.", filename);
2000 0 : return true;
2001 : }
2002 1 : if (nTileIndexFiles > 100 * 1024 * 1024)
2003 : {
2004 0 : CPLError(CE_Failure, CPLE_AppDefined,
2005 : "Too large feature count in tile index");
2006 0 : return false;
2007 : }
2008 :
2009 5 : for (auto &&poFeature : poLayer)
2010 : {
2011 4 : aosList.AddString(poFeature->GetFieldAsString(ti_field));
2012 : }
2013 : }
2014 : else
2015 : {
2016 44 : aosList.AddString(filename);
2017 : }
2018 :
2019 45 : return true;
2020 : }
2021 :
2022 : /************************************************************************/
2023 : /* GDALBuildVRTOptions */
2024 : /************************************************************************/
2025 :
2026 : /** Options for use with GDALBuildVRT(). GDALBuildVRTOptions* must be allocated
2027 : * and freed with GDALBuildVRTOptionsNew() and GDALBuildVRTOptionsFree()
2028 : * respectively.
2029 : */
2030 : struct GDALBuildVRTOptions
2031 : {
2032 : std::string osProgramName = "gdalbuildvrt";
2033 : std::string osTileIndex = "location";
2034 : bool bStrict = false;
2035 : std::string osResolution{};
2036 : bool bSeparate = false;
2037 : bool bAllowProjectionDifference = false;
2038 : double we_res = 0;
2039 : double ns_res = 0;
2040 : bool bTargetAlignedPixels = false;
2041 : double xmin = 0;
2042 : double ymin = 0;
2043 : double xmax = 0;
2044 : double ymax = 0;
2045 : bool bAddAlpha = false;
2046 : bool bHideNoData = false;
2047 : int nSubdataset = -1;
2048 : std::string osSrcNoData{};
2049 : std::string osVRTNoData{};
2050 : std::string osOutputSRS{};
2051 : std::vector<int> anSelectedBandList{};
2052 : std::string osResampling{};
2053 : CPLStringList aosOpenOptions{};
2054 : CPLStringList aosCreateOptions{};
2055 : bool bUseSrcMaskBand = true;
2056 : bool bNoDataFromMask = false;
2057 : double dfMaskValueThreshold = 0;
2058 : bool bWriteAbsolutePath = false;
2059 : std::string osPixelFunction{};
2060 : CPLStringList aosPixelFunctionArgs{};
2061 :
2062 : /*! allow or suppress progress monitor and other non-error output */
2063 : bool bQuiet = true;
2064 :
2065 : /*! the progress function to use */
2066 : GDALProgressFunc pfnProgress = GDALDummyProgress;
2067 :
2068 : /*! pointer to the progress data variable */
2069 : void *pProgressData = nullptr;
2070 : };
2071 :
2072 : /************************************************************************/
2073 : /* GDALBuildVRT() */
2074 : /************************************************************************/
2075 :
2076 : /* clang-format off */
2077 : /**
2078 : * Build a VRT from a list of datasets.
2079 : *
2080 : * This is the equivalent of the
2081 : * <a href="/programs/gdalbuildvrt.html">gdalbuildvrt</a> utility.
2082 : *
2083 : * GDALBuildVRTOptions* must be allocated and freed with
2084 : * GDALBuildVRTOptionsNew() and GDALBuildVRTOptionsFree() respectively. pahSrcDS
2085 : * and papszSrcDSNames cannot be used at the same time.
2086 : *
2087 : * @param pszDest the destination dataset path.
2088 : * @param nSrcCount the number of input datasets.
2089 : * @param pahSrcDS the list of input datasets (or NULL, exclusive with
2090 : * papszSrcDSNames). For practical purposes, the type
2091 : * of this argument should be considered as "const GDALDatasetH* const*", that
2092 : * is neither the array nor its values are mutated by this function.
2093 : * @param papszSrcDSNames the list of input dataset names (or NULL, exclusive
2094 : * with pahSrcDS)
2095 : * @param psOptionsIn the options struct returned by GDALBuildVRTOptionsNew() or
2096 : * NULL.
2097 : * @param pbUsageError pointer to a integer output variable to store if any
2098 : * usage error has occurred.
2099 : * @return the output dataset (new dataset that must be closed using
2100 : * GDALClose()) or NULL in case of error. If using pahSrcDS, the returned VRT
2101 : * dataset has a reference to each pahSrcDS[] element. Hence pahSrcDS[] elements
2102 : * should be closed after the returned dataset if using GDALClose().
2103 : * A safer alternative is to use GDALReleaseDataset() instead of using
2104 : * GDALClose(), in which case you can close datasets in any order.
2105 :
2106 : *
2107 : * @since GDAL 2.1
2108 : */
2109 : /* clang-format on */
2110 :
2111 253 : GDALDatasetH GDALBuildVRT(const char *pszDest, int nSrcCount,
2112 : GDALDatasetH *pahSrcDS,
2113 : const char *const *papszSrcDSNames,
2114 : const GDALBuildVRTOptions *psOptionsIn,
2115 : int *pbUsageError)
2116 : {
2117 253 : if (pszDest == nullptr)
2118 0 : pszDest = "";
2119 :
2120 253 : if (nSrcCount == 0)
2121 : {
2122 0 : CPLError(CE_Failure, CPLE_AppDefined, "No input dataset specified.");
2123 :
2124 0 : if (pbUsageError)
2125 0 : *pbUsageError = TRUE;
2126 0 : return nullptr;
2127 : }
2128 :
2129 : // cppcheck-suppress unreadVariable
2130 : GDALBuildVRTOptions sOptions(psOptionsIn ? *psOptionsIn
2131 506 : : GDALBuildVRTOptions());
2132 :
2133 8 : if (sOptions.we_res != 0 && sOptions.ns_res != 0 &&
2134 261 : !sOptions.osResolution.empty() &&
2135 0 : !EQUAL(sOptions.osResolution.c_str(), "user"))
2136 : {
2137 0 : CPLError(CE_Failure, CPLE_NotSupported,
2138 : "-tr option is not compatible with -resolution %s",
2139 : sOptions.osResolution.c_str());
2140 0 : if (pbUsageError)
2141 0 : *pbUsageError = TRUE;
2142 0 : return nullptr;
2143 : }
2144 :
2145 253 : if (sOptions.bTargetAlignedPixels && sOptions.we_res == 0 &&
2146 1 : sOptions.ns_res == 0)
2147 : {
2148 1 : CPLError(CE_Failure, CPLE_NotSupported,
2149 : "-tap option cannot be used without using -tr");
2150 1 : if (pbUsageError)
2151 1 : *pbUsageError = TRUE;
2152 1 : return nullptr;
2153 : }
2154 :
2155 252 : if (sOptions.bAddAlpha && sOptions.bSeparate)
2156 : {
2157 0 : CPLError(CE_Failure, CPLE_NotSupported,
2158 : "-addalpha option is not compatible with -separate.");
2159 0 : if (pbUsageError)
2160 0 : *pbUsageError = TRUE;
2161 0 : return nullptr;
2162 : }
2163 :
2164 252 : ResolutionStrategy eStrategy = AVERAGE_RESOLUTION;
2165 303 : if (sOptions.osResolution.empty() ||
2166 51 : EQUAL(sOptions.osResolution.c_str(), "user"))
2167 : {
2168 201 : if (sOptions.we_res != 0 || sOptions.ns_res != 0)
2169 8 : eStrategy = USER_RESOLUTION;
2170 193 : else if (EQUAL(sOptions.osResolution.c_str(), "user"))
2171 : {
2172 0 : CPLError(CE_Failure, CPLE_NotSupported,
2173 : "-tr option must be used with -resolution user.");
2174 0 : if (pbUsageError)
2175 0 : *pbUsageError = TRUE;
2176 0 : return nullptr;
2177 : }
2178 : }
2179 51 : else if (EQUAL(sOptions.osResolution.c_str(), "average"))
2180 1 : eStrategy = AVERAGE_RESOLUTION;
2181 50 : else if (EQUAL(sOptions.osResolution.c_str(), "highest"))
2182 1 : eStrategy = HIGHEST_RESOLUTION;
2183 49 : else if (EQUAL(sOptions.osResolution.c_str(), "lowest"))
2184 1 : eStrategy = LOWEST_RESOLUTION;
2185 48 : else if (EQUAL(sOptions.osResolution.c_str(), "same"))
2186 38 : eStrategy = SAME_RESOLUTION;
2187 10 : else if (EQUAL(sOptions.osResolution.c_str(), "common"))
2188 10 : eStrategy = COMMON_RESOLUTION;
2189 :
2190 : /* If -srcnodata is specified, use it as the -vrtnodata if the latter is not
2191 : */
2192 : /* specified */
2193 252 : if (!sOptions.osSrcNoData.empty() && sOptions.osVRTNoData.empty())
2194 3 : sOptions.osVRTNoData = sOptions.osSrcNoData;
2195 :
2196 : VRTBuilder oBuilder(
2197 252 : sOptions.bStrict, pszDest, nSrcCount, papszSrcDSNames, pahSrcDS,
2198 252 : sOptions.anSelectedBandList.empty()
2199 : ? nullptr
2200 19 : : sOptions.anSelectedBandList.data(),
2201 252 : static_cast<int>(sOptions.anSelectedBandList.size()), eStrategy,
2202 252 : sOptions.we_res, sOptions.ns_res, sOptions.bTargetAlignedPixels,
2203 : sOptions.xmin, sOptions.ymin, sOptions.xmax, sOptions.ymax,
2204 252 : sOptions.bSeparate, sOptions.bAllowProjectionDifference,
2205 252 : sOptions.bAddAlpha, sOptions.bHideNoData, sOptions.nSubdataset,
2206 259 : sOptions.osSrcNoData.empty() ? nullptr : sOptions.osSrcNoData.c_str(),
2207 24 : sOptions.osVRTNoData.empty() ? nullptr : sOptions.osVRTNoData.c_str(),
2208 252 : sOptions.bUseSrcMaskBand, sOptions.bNoDataFromMask,
2209 : sOptions.dfMaskValueThreshold,
2210 253 : sOptions.osOutputSRS.empty() ? nullptr : sOptions.osOutputSRS.c_str(),
2211 267 : sOptions.osResampling.empty() ? nullptr : sOptions.osResampling.c_str(),
2212 252 : sOptions.osPixelFunction.empty() ? nullptr
2213 3 : : sOptions.osPixelFunction.c_str(),
2214 252 : sOptions.aosPixelFunctionArgs, sOptions.aosOpenOptions.List(),
2215 1287 : sOptions.aosCreateOptions, sOptions.bWriteAbsolutePath);
2216 252 : oBuilder.m_osProgramName = sOptions.osProgramName;
2217 :
2218 252 : return GDALDataset::ToHandle(
2219 504 : oBuilder.Build(sOptions.pfnProgress, sOptions.pProgressData).release());
2220 : }
2221 :
2222 : /************************************************************************/
2223 : /* SanitizeSRS */
2224 : /************************************************************************/
2225 :
2226 1 : static char *SanitizeSRS(const char *pszUserInput)
2227 :
2228 : {
2229 : OGRSpatialReferenceH hSRS;
2230 1 : char *pszResult = nullptr;
2231 :
2232 1 : CPLErrorReset();
2233 :
2234 1 : hSRS = OSRNewSpatialReference(nullptr);
2235 1 : if (OSRSetFromUserInput(hSRS, pszUserInput) == OGRERR_NONE)
2236 1 : OSRExportToWkt(hSRS, &pszResult);
2237 : else
2238 : {
2239 0 : CPLError(CE_Failure, CPLE_AppDefined, "Translating SRS failed:\n%s",
2240 : pszUserInput);
2241 : }
2242 :
2243 1 : OSRDestroySpatialReference(hSRS);
2244 :
2245 1 : return pszResult;
2246 : }
2247 :
2248 : /************************************************************************/
2249 : /* GDALBuildVRTOptionsGetParser() */
2250 : /************************************************************************/
2251 :
2252 : static std::unique_ptr<GDALArgumentParser>
2253 256 : GDALBuildVRTOptionsGetParser(GDALBuildVRTOptions *psOptions,
2254 : GDALBuildVRTOptionsForBinary *psOptionsForBinary)
2255 : {
2256 : auto argParser = std::make_unique<GDALArgumentParser>(
2257 256 : "gdalbuildvrt", /* bForBinary=*/psOptionsForBinary != nullptr);
2258 :
2259 256 : argParser->add_description(_("Builds a VRT from a list of datasets."));
2260 :
2261 256 : argParser->add_epilog(_(
2262 : "\n"
2263 : "e.g.\n"
2264 : " % gdalbuildvrt doq_index.vrt doq/*.tif\n"
2265 : " % gdalbuildvrt -input_file_list my_list.txt doq_index.vrt\n"
2266 : "\n"
2267 : "NOTES:\n"
2268 : " o With -separate, each files goes into a separate band in the VRT "
2269 : "band.\n"
2270 : " Otherwise, the files are considered as tiles of a larger mosaic.\n"
2271 : " o -b option selects a band to add into vrt. Multiple bands can be "
2272 : "listed.\n"
2273 : " By default all bands are queried.\n"
2274 : " o The default tile index field is 'location' unless otherwise "
2275 : "specified by\n"
2276 : " -tileindex.\n"
2277 : " o In case the resolution of all input files is not the same, the "
2278 : "-resolution\n"
2279 : " flag enable the user to control the way the output resolution is "
2280 : "computed.\n"
2281 : " Average is the default.\n"
2282 : " o Input files may be any valid GDAL dataset or a GDAL raster tile "
2283 : "index.\n"
2284 : " o For a GDAL raster tile index, all entries will be added to the "
2285 : "VRT.\n"
2286 : " o If one GDAL dataset is made of several subdatasets and has 0 "
2287 : "raster bands,\n"
2288 : " its datasets will be added to the VRT rather than the dataset "
2289 : "itself.\n"
2290 : " Single subdataset could be selected by its number using the -sd "
2291 : "option.\n"
2292 : " o By default, only datasets of same projection and band "
2293 : "characteristics\n"
2294 : " may be added to the VRT.\n"
2295 : "\n"
2296 : "For more details, consult "
2297 256 : "https://gdal.org/programs/gdalbuildvrt.html"));
2298 :
2299 : argParser->add_quiet_argument(
2300 256 : psOptionsForBinary ? &psOptionsForBinary->bQuiet : nullptr);
2301 :
2302 : {
2303 256 : auto &group = argParser->add_mutually_exclusive_group();
2304 :
2305 256 : group.add_argument("-strict")
2306 256 : .flag()
2307 256 : .store_into(psOptions->bStrict)
2308 256 : .help(_("Turn warnings as failures."));
2309 :
2310 256 : group.add_argument("-non_strict")
2311 256 : .flag()
2312 0 : .action([psOptions](const std::string &)
2313 256 : { psOptions->bStrict = false; })
2314 : .help(_("Skip source datasets that have issues with warnings, and "
2315 256 : "continue processing."));
2316 : }
2317 :
2318 256 : argParser->add_argument("-tile_index")
2319 512 : .metavar("<field_name>")
2320 256 : .store_into(psOptions->osTileIndex)
2321 : .help(_("Use the specified value as the tile index field, instead of "
2322 256 : "the default value which is 'location'."));
2323 :
2324 256 : argParser->add_argument("-resolution")
2325 512 : .metavar("user|average|common|highest|lowest|same")
2326 : .action(
2327 310 : [psOptions](const std::string &s)
2328 : {
2329 51 : psOptions->osResolution = s;
2330 51 : if (!EQUAL(psOptions->osResolution.c_str(), "user") &&
2331 51 : !EQUAL(psOptions->osResolution.c_str(), "average") &&
2332 50 : !EQUAL(psOptions->osResolution.c_str(), "highest") &&
2333 49 : !EQUAL(psOptions->osResolution.c_str(), "lowest") &&
2334 112 : !EQUAL(psOptions->osResolution.c_str(), "same") &&
2335 10 : !EQUAL(psOptions->osResolution.c_str(), "common"))
2336 : {
2337 : throw std::invalid_argument(
2338 : CPLSPrintf("Illegal resolution value (%s).",
2339 0 : psOptions->osResolution.c_str()));
2340 : }
2341 307 : })
2342 256 : .help(_("Control the way the output resolution is computed."));
2343 :
2344 256 : argParser->add_argument("-tr")
2345 512 : .metavar("<xres> <yres>")
2346 256 : .nargs(2)
2347 256 : .scan<'g', double>()
2348 256 : .help(_("Set target resolution."));
2349 :
2350 256 : if (psOptionsForBinary)
2351 : {
2352 20 : argParser->add_argument("-input_file_list")
2353 40 : .metavar("<filename>")
2354 : .action(
2355 5 : [psOptions, psOptionsForBinary](const std::string &s)
2356 : {
2357 1 : const char *input_file_list = s.c_str();
2358 : auto f = VSIVirtualHandleUniquePtr(
2359 2 : VSIFOpenL(input_file_list, "r"));
2360 1 : if (f)
2361 : {
2362 : while (1)
2363 : {
2364 5 : const char *filename = CPLReadLineL(f.get());
2365 5 : if (filename == nullptr)
2366 1 : break;
2367 4 : if (!add_file_to_list(
2368 : filename, psOptions->osTileIndex.c_str(),
2369 4 : psOptionsForBinary->aosSrcFiles))
2370 : {
2371 : throw std::invalid_argument(
2372 0 : std::string("Cannot add ")
2373 0 : .append(filename)
2374 0 : .append(" to input file list"));
2375 : }
2376 4 : }
2377 : }
2378 21 : })
2379 20 : .help(_("Text file with an input filename on each line"));
2380 : }
2381 :
2382 : {
2383 256 : auto &group = argParser->add_mutually_exclusive_group();
2384 :
2385 256 : group.add_argument("-separate")
2386 256 : .flag()
2387 256 : .store_into(psOptions->bSeparate)
2388 256 : .help(_("Place each input file into a separate band."));
2389 :
2390 256 : group.add_argument("-pixel-function")
2391 512 : .metavar("<function>")
2392 : .action(
2393 9 : [psOptions](const std::string &s)
2394 : {
2395 : auto *poPixFun =
2396 5 : VRTDerivedRasterBand::GetPixelFunction(s.c_str());
2397 5 : if (poPixFun == nullptr)
2398 : {
2399 : throw std::invalid_argument(
2400 1 : s + " is not a registered pixel function.");
2401 : }
2402 :
2403 4 : psOptions->osPixelFunction = s;
2404 260 : })
2405 :
2406 256 : .help("Function to calculate value from overlapping inputs");
2407 : }
2408 :
2409 256 : argParser->add_argument("-pixel-function-arg")
2410 512 : .metavar("<NAME>=<VALUE>")
2411 256 : .append()
2412 2 : .action([psOptions](const std::string &s)
2413 258 : { psOptions->aosPixelFunctionArgs.AddString(s); })
2414 256 : .help(_("Pixel function argument(s)"));
2415 :
2416 256 : argParser->add_argument("-allow_projection_difference")
2417 256 : .flag()
2418 256 : .store_into(psOptions->bAllowProjectionDifference)
2419 : .help(_("Accept source files not in the same projection (but without "
2420 256 : "reprojecting them!)."));
2421 :
2422 256 : argParser->add_argument("-sd")
2423 512 : .metavar("<n>")
2424 256 : .store_into(psOptions->nSubdataset)
2425 : .help(_("Use subdataset of specified index (starting at 1), instead of "
2426 256 : "the source dataset itself."));
2427 :
2428 256 : argParser->add_argument("-tap")
2429 256 : .flag()
2430 256 : .store_into(psOptions->bTargetAlignedPixels)
2431 : .help(_("Align the coordinates of the extent of the output file to the "
2432 256 : "values of the resolution."));
2433 :
2434 256 : argParser->add_argument("-te")
2435 512 : .metavar("<xmin> <ymin> <xmax> <ymax>")
2436 256 : .nargs(4)
2437 256 : .scan<'g', double>()
2438 256 : .help(_("Set georeferenced extents of output file to be created."));
2439 :
2440 256 : argParser->add_argument("-addalpha")
2441 256 : .flag()
2442 256 : .store_into(psOptions->bAddAlpha)
2443 : .help(_("Adds an alpha mask band to the VRT when the source raster "
2444 256 : "have none."));
2445 :
2446 256 : argParser->add_argument("-b")
2447 512 : .metavar("<band>")
2448 256 : .append()
2449 256 : .store_into(psOptions->anSelectedBandList)
2450 256 : .help(_("Specify input band(s) number."));
2451 :
2452 256 : argParser->add_argument("-hidenodata")
2453 256 : .flag()
2454 256 : .store_into(psOptions->bHideNoData)
2455 256 : .help(_("Makes the VRT band not report the NoData."));
2456 :
2457 256 : if (psOptionsForBinary)
2458 : {
2459 20 : argParser->add_argument("-overwrite")
2460 20 : .flag()
2461 20 : .store_into(psOptionsForBinary->bOverwrite)
2462 20 : .help(_("Overwrite the VRT if it already exists."));
2463 : }
2464 :
2465 256 : argParser->add_argument("-srcnodata")
2466 512 : .metavar("\"<value>[ <value>]...\"")
2467 256 : .store_into(psOptions->osSrcNoData)
2468 256 : .help(_("Set nodata values for input bands."));
2469 :
2470 256 : argParser->add_argument("-vrtnodata")
2471 512 : .metavar("\"<value>[ <value>]...\"")
2472 256 : .store_into(psOptions->osVRTNoData)
2473 256 : .help(_("Set nodata values at the VRT band level."));
2474 :
2475 256 : argParser->add_argument("-a_srs")
2476 512 : .metavar("<srs_def>")
2477 : .action(
2478 2 : [psOptions](const std::string &s)
2479 : {
2480 1 : char *pszSRS = SanitizeSRS(s.c_str());
2481 1 : if (pszSRS == nullptr)
2482 : {
2483 0 : throw std::invalid_argument("Invalid value for -a_srs");
2484 : }
2485 1 : psOptions->osOutputSRS = pszSRS;
2486 1 : CPLFree(pszSRS);
2487 257 : })
2488 256 : .help(_("Override the projection for the output file.."));
2489 :
2490 256 : argParser->add_argument("-r")
2491 512 : .metavar("nearest|bilinear|cubic|cubicspline|lanczos|average|mode")
2492 256 : .store_into(psOptions->osResampling)
2493 256 : .help(_("Resampling algorithm."));
2494 :
2495 256 : argParser->add_open_options_argument(&psOptions->aosOpenOptions);
2496 :
2497 256 : argParser->add_creation_options_argument(psOptions->aosCreateOptions);
2498 :
2499 256 : argParser->add_argument("-write_absolute_path")
2500 256 : .flag()
2501 256 : .store_into(psOptions->bWriteAbsolutePath)
2502 : .help(_("Write the absolute path of the raster files in the tile index "
2503 256 : "file."));
2504 :
2505 256 : argParser->add_argument("-ignore_srcmaskband")
2506 256 : .flag()
2507 0 : .action([psOptions](const std::string &)
2508 256 : { psOptions->bUseSrcMaskBand = false; })
2509 256 : .help(_("Cause mask band of sources will not be taken into account."));
2510 :
2511 256 : argParser->add_argument("-nodata_max_mask_threshold")
2512 512 : .metavar("<threshold>")
2513 256 : .scan<'g', double>()
2514 : .action(
2515 18 : [psOptions](const std::string &s)
2516 : {
2517 9 : psOptions->bNoDataFromMask = true;
2518 9 : psOptions->dfMaskValueThreshold = CPLAtofM(s.c_str());
2519 256 : })
2520 : .help(_("Replaces the value of the source with the value of -vrtnodata "
2521 : "when the value of the mask band of the source is less or "
2522 256 : "equal to the threshold."));
2523 :
2524 256 : argParser->add_argument("-program_name")
2525 256 : .store_into(psOptions->osProgramName)
2526 256 : .hidden();
2527 :
2528 256 : if (psOptionsForBinary)
2529 : {
2530 20 : if (psOptionsForBinary->osDstFilename.empty())
2531 : {
2532 : // We normally go here, unless undocumented -o switch is used
2533 20 : argParser->add_argument("vrt_dataset_name")
2534 40 : .metavar("<vrt_dataset_name>")
2535 20 : .store_into(psOptionsForBinary->osDstFilename)
2536 20 : .help(_("Output VRT."));
2537 : }
2538 :
2539 20 : argParser->add_argument("src_dataset_name")
2540 40 : .metavar("<src_dataset_name>")
2541 20 : .nargs(argparse::nargs_pattern::any)
2542 : .action(
2543 41 : [psOptions, psOptionsForBinary](const std::string &s)
2544 : {
2545 41 : if (!add_file_to_list(s.c_str(),
2546 : psOptions->osTileIndex.c_str(),
2547 41 : psOptionsForBinary->aosSrcFiles))
2548 : {
2549 : throw std::invalid_argument(
2550 0 : std::string("Cannot add ")
2551 0 : .append(s)
2552 0 : .append(" to input file list"));
2553 : }
2554 61 : })
2555 20 : .help(_("Input dataset(s)."));
2556 : }
2557 :
2558 256 : return argParser;
2559 : }
2560 :
2561 : /************************************************************************/
2562 : /* GDALBuildVRTGetParserUsage() */
2563 : /************************************************************************/
2564 :
2565 1 : std::string GDALBuildVRTGetParserUsage()
2566 : {
2567 : try
2568 : {
2569 2 : GDALBuildVRTOptions sOptions;
2570 2 : GDALBuildVRTOptionsForBinary sOptionsForBinary;
2571 : auto argParser =
2572 2 : GDALBuildVRTOptionsGetParser(&sOptions, &sOptionsForBinary);
2573 1 : return argParser->usage();
2574 : }
2575 0 : catch (const std::exception &err)
2576 : {
2577 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unexpected exception: %s",
2578 0 : err.what());
2579 0 : return std::string();
2580 : }
2581 : }
2582 :
2583 : /************************************************************************/
2584 : /* GDALBuildVRTOptionsNew() */
2585 : /************************************************************************/
2586 :
2587 : /**
2588 : * Allocates a GDALBuildVRTOptions struct.
2589 : *
2590 : * @param papszArgv NULL terminated list of options (potentially including
2591 : * filename and open options too), or NULL. The accepted options are the ones of
2592 : * the <a href="/programs/gdalbuildvrt.html">gdalbuildvrt</a> utility.
2593 : * @param psOptionsForBinary (output) may be NULL (and should generally be
2594 : * NULL), otherwise (gdalbuildvrt_bin.cpp use case) must be allocated with
2595 : * GDALBuildVRTOptionsForBinaryNew() prior to this function. Will be filled
2596 : * with potentially present filename, open options,...
2597 : * @return pointer to the allocated GDALBuildVRTOptions struct. Must be freed
2598 : * with GDALBuildVRTOptionsFree().
2599 : *
2600 : * @since GDAL 2.1
2601 : */
2602 :
2603 : GDALBuildVRTOptions *
2604 255 : GDALBuildVRTOptionsNew(char **papszArgv,
2605 : GDALBuildVRTOptionsForBinary *psOptionsForBinary)
2606 : {
2607 510 : auto psOptions = std::make_unique<GDALBuildVRTOptions>();
2608 :
2609 510 : CPLStringList aosArgv;
2610 255 : const int nArgc = CSLCount(papszArgv);
2611 1268 : for (int i = 0;
2612 1268 : i < nArgc && papszArgv != nullptr && papszArgv[i] != nullptr; i++)
2613 : {
2614 1013 : if (psOptionsForBinary && EQUAL(papszArgv[i], "-o") && i + 1 < nArgc &&
2615 0 : papszArgv[i + 1] != nullptr)
2616 : {
2617 : // Undocumented alternate way of specifying the destination file
2618 0 : psOptionsForBinary->osDstFilename = papszArgv[i + 1];
2619 0 : ++i;
2620 : }
2621 : // argparser will be confused if the value of a string argument
2622 : // starts with a negative sign.
2623 1013 : else if (EQUAL(papszArgv[i], "-srcnodata") && i + 1 < nArgc)
2624 : {
2625 7 : ++i;
2626 7 : psOptions->osSrcNoData = papszArgv[i];
2627 : }
2628 : // argparser will be confused if the value of a string argument
2629 : // starts with a negative sign.
2630 1006 : else if (EQUAL(papszArgv[i], "-vrtnodata") && i + 1 < nArgc)
2631 : {
2632 21 : ++i;
2633 21 : psOptions->osVRTNoData = papszArgv[i];
2634 : }
2635 :
2636 : else
2637 : {
2638 985 : aosArgv.AddString(papszArgv[i]);
2639 : }
2640 : }
2641 :
2642 : try
2643 : {
2644 : auto argParser =
2645 510 : GDALBuildVRTOptionsGetParser(psOptions.get(), psOptionsForBinary);
2646 :
2647 255 : argParser->parse_args_without_binary_name(aosArgv.List());
2648 :
2649 261 : if (auto adfTargetRes = argParser->present<std::vector<double>>("-tr"))
2650 : {
2651 8 : psOptions->we_res = (*adfTargetRes)[0];
2652 8 : psOptions->ns_res = (*adfTargetRes)[1];
2653 : }
2654 :
2655 267 : if (auto oTE = argParser->present<std::vector<double>>("-te"))
2656 : {
2657 14 : psOptions->xmin = (*oTE)[0];
2658 14 : psOptions->ymin = (*oTE)[1];
2659 14 : psOptions->xmax = (*oTE)[2];
2660 14 : psOptions->ymax = (*oTE)[3];
2661 : }
2662 :
2663 503 : if (psOptions->osPixelFunction.empty() &&
2664 250 : !psOptions->aosPixelFunctionArgs.empty())
2665 : {
2666 : throw std::runtime_error(
2667 1 : "Pixel function arguments provided without a pixel function");
2668 : }
2669 :
2670 252 : return psOptions.release();
2671 : }
2672 3 : catch (const std::exception &err)
2673 : {
2674 3 : CPLError(CE_Failure, CPLE_AppDefined, "%s", err.what());
2675 3 : return nullptr;
2676 : }
2677 : }
2678 :
2679 : /************************************************************************/
2680 : /* GDALBuildVRTOptionsFree() */
2681 : /************************************************************************/
2682 :
2683 : /**
2684 : * Frees the GDALBuildVRTOptions struct.
2685 : *
2686 : * @param psOptions the options struct for GDALBuildVRT().
2687 : *
2688 : * @since GDAL 2.1
2689 : */
2690 :
2691 251 : void GDALBuildVRTOptionsFree(GDALBuildVRTOptions *psOptions)
2692 : {
2693 251 : delete psOptions;
2694 251 : }
2695 :
2696 : /************************************************************************/
2697 : /* GDALBuildVRTOptionsSetProgress() */
2698 : /************************************************************************/
2699 :
2700 : /**
2701 : * Set a progress function.
2702 : *
2703 : * @param psOptions the options struct for GDALBuildVRT().
2704 : * @param pfnProgress the progress callback.
2705 : * @param pProgressData the user data for the progress callback.
2706 : *
2707 : * @since GDAL 2.1
2708 : */
2709 :
2710 20 : void GDALBuildVRTOptionsSetProgress(GDALBuildVRTOptions *psOptions,
2711 : GDALProgressFunc pfnProgress,
2712 : void *pProgressData)
2713 : {
2714 20 : psOptions->pfnProgress = pfnProgress ? pfnProgress : GDALDummyProgress;
2715 20 : psOptions->pProgressData = pProgressData;
2716 20 : if (pfnProgress == GDALTermProgress)
2717 19 : psOptions->bQuiet = false;
2718 20 : }
|
