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