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