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