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