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