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