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