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