1 : /******************************************************************************
2 : * $Id: gdalbuildvrt.cpp 25026 2012-10-01 21:23:05Z rouault $
3 : *
4 : * Project: GDAL Utilities
5 : * Purpose: Commandline application to build VRT datasets from raster products or content of SHP tile index
6 : * Author: Even Rouault, even.rouault at mines-paris.org
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2007, Even Rouault
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 "gdal_proxy.h"
31 : #include "cpl_string.h"
32 : #include "vrt/gdal_vrt.h"
33 : #include "vrt/vrtdataset.h"
34 :
35 : #ifdef OGR_ENABLED
36 : #include "ogr_api.h"
37 : #endif
38 : #include "ogr_srs_api.h"
39 :
40 : CPL_CVSID("$Id: gdalbuildvrt.cpp 25026 2012-10-01 21:23:05Z rouault $");
41 :
42 : #define GEOTRSFRM_TOPLEFT_X 0
43 : #define GEOTRSFRM_WE_RES 1
44 : #define GEOTRSFRM_ROTATION_PARAM1 2
45 : #define GEOTRSFRM_TOPLEFT_Y 3
46 : #define GEOTRSFRM_ROTATION_PARAM2 4
47 : #define GEOTRSFRM_NS_RES 5
48 :
49 : typedef enum
50 : {
51 : LOWEST_RESOLUTION,
52 : HIGHEST_RESOLUTION,
53 : AVERAGE_RESOLUTION,
54 : USER_RESOLUTION
55 : } ResolutionStrategy;
56 :
57 : typedef struct
58 : {
59 : int isFileOK;
60 : int nRasterXSize;
61 : int nRasterYSize;
62 : double adfGeoTransform[6];
63 : int nBlockXSize;
64 : int nBlockYSize;
65 : GDALDataType firstBandType;
66 : int* panHasNoData;
67 : double* padfNoDataValues;
68 : int bHasDatasetMask;
69 : int nMaskBlockXSize;
70 : int nMaskBlockYSize;
71 : } DatasetProperty;
72 :
73 : typedef struct
74 : {
75 : GDALColorInterp colorInterpretation;
76 : GDALDataType dataType;
77 : GDALColorTableH colorTable;
78 : int bHasNoData;
79 : double noDataValue;
80 : } BandProperty;
81 :
82 : /************************************************************************/
83 : /* Usage() */
84 : /************************************************************************/
85 :
86 1 : static void Usage()
87 :
88 : {
89 : fprintf(stdout, "%s",
90 : "Usage: gdalbuildvrt [-tileindex field_name] [-resolution {highest|lowest|average|user}]\n"
91 : " [-tr xres yres] [-tap] [-separate] [-allow_projection_difference] [-q]\n"
92 : " [-te xmin ymin xmax ymax] [-addalpha] [-hidenodata] \n"
93 : " [-srcnodata \"value [value...]\"] [-vrtnodata \"value [value...]\"] \n"
94 : " [-a_srs srs_def]\n"
95 : " [-input_file_list my_liste.txt] [-overwrite] output.vrt [gdalfile]*\n"
96 : "\n"
97 : "eg.\n"
98 : " % gdalbuildvrt doq_index.vrt doq/*.tif\n"
99 : " % gdalbuildvrt -input_file_list my_liste.txt doq_index.vrt\n"
100 : "\n"
101 : "NOTES:\n"
102 : " o With -separate, each files goes into a separate band in the VRT band. Otherwise,\n"
103 : " the files are considered as tiles of a larger mosaic.\n"
104 : " o The default tile index field is 'location' unless otherwise specified by -tileindex.\n"
105 : " o In case the resolution of all input files is not the same, the -resolution flag.\n"
106 : " enable the user to control the way the output resolution is computed. average is the default.\n"
107 : " o Input files may be any valid GDAL dataset or a GDAL raster tile index.\n"
108 : " o For a GDAL raster tile index, all entries will be added to the VRT.\n"
109 : " o If one GDAL dataset is made of several subdatasets and has 0 raster bands, its\n"
110 : " datasets will be added to the VRT rather than the dataset itself.\n"
111 : " o By default, only datasets of same projection and band characteristics may be added to the VRT.\n"
112 1 : );
113 1 : exit( 1 );
114 : }
115 :
116 :
117 : /************************************************************************/
118 : /* GetSrcDstWin() */
119 : /************************************************************************/
120 :
121 41 : int GetSrcDstWin(DatasetProperty* psDP,
122 : double we_res, double ns_res,
123 : double minX, double minY, double maxX, double maxY,
124 : int* pnSrcXOff, int* pnSrcYOff, int* pnSrcXSize, int* pnSrcYSize,
125 : int* pnDstXOff, int* pnDstYOff, int* pnDstXSize, int* pnDstYSize)
126 : {
127 : /* Check that the destination bounding box intersects the source bounding box */
128 82 : if ( psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_X] +
129 : psDP->nRasterXSize *
130 41 : psDP->adfGeoTransform[GEOTRSFRM_WE_RES] < minX )
131 0 : return FALSE;
132 41 : if ( psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_X] > maxX )
133 0 : return FALSE;
134 82 : if ( psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_Y] +
135 : psDP->nRasterYSize *
136 41 : psDP->adfGeoTransform[GEOTRSFRM_NS_RES] > maxY )
137 0 : return FALSE;
138 41 : if ( psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_Y] < minY )
139 0 : return FALSE;
140 :
141 41 : *pnSrcXSize = psDP->nRasterXSize;
142 41 : *pnSrcYSize = psDP->nRasterYSize;
143 41 : if ( psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_X] < minX )
144 : {
145 1 : *pnSrcXOff = (int)((minX - psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_X]) /
146 1 : psDP->adfGeoTransform[GEOTRSFRM_WE_RES] + 0.5);
147 1 : *pnDstXOff = 0;
148 : }
149 : else
150 : {
151 40 : *pnSrcXOff = 0;
152 : *pnDstXOff = (int)
153 40 : (0.5 + (psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_X] - minX) / we_res);
154 : }
155 41 : if ( maxY < psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_Y])
156 : {
157 1 : *pnSrcYOff = (int)((psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_Y] - maxY) /
158 1 : -psDP->adfGeoTransform[GEOTRSFRM_NS_RES] + 0.5);
159 1 : *pnDstYOff = 0;
160 : }
161 : else
162 : {
163 40 : *pnSrcYOff = 0;
164 : *pnDstYOff = (int)
165 40 : (0.5 + (maxY - psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_Y]) / -ns_res);
166 : }
167 : *pnDstXSize = (int)
168 : (0.5 + psDP->nRasterXSize *
169 41 : psDP->adfGeoTransform[GEOTRSFRM_WE_RES] / we_res);
170 : *pnDstYSize = (int)
171 : (0.5 + psDP->nRasterYSize *
172 41 : psDP->adfGeoTransform[GEOTRSFRM_NS_RES] / ns_res);
173 :
174 41 : return TRUE;
175 : }
176 :
177 : /************************************************************************/
178 : /* VRTBuilder */
179 : /************************************************************************/
180 :
181 : class VRTBuilder
182 : {
183 : /* Input parameters */
184 : char *pszOutputFilename;
185 : int nInputFiles;
186 : char **ppszInputFilenames;
187 : ResolutionStrategy resolutionStrategy;
188 : double we_res;
189 : double ns_res;
190 : int bTargetAlignedPixels;
191 : double minX;
192 : double minY;
193 : double maxX;
194 : double maxY;
195 : int bSeparate;
196 : int bAllowProjectionDifference;
197 : int bAddAlpha;
198 : int bHideNoData;
199 : char *pszSrcNoData;
200 : char *pszVRTNoData;
201 : char *pszOutputSRS;
202 :
203 : /* Internal variables */
204 : char *pszProjectionRef;
205 : int nBands;
206 : BandProperty *pasBandProperties;
207 : int bFirst;
208 : int bHasGeoTransform;
209 : int nRasterXSize;
210 : int nRasterYSize;
211 : DatasetProperty *pasDatasetProperties;
212 : int bUserExtent;
213 : int bAllowSrcNoData;
214 : double *padfSrcNoData;
215 : int nSrcNoDataCount;
216 : int bAllowVRTNoData;
217 : double *padfVRTNoData;
218 : int nVRTNoDataCount;
219 : int bHasRunBuild;
220 : int bHasDatasetMask;
221 :
222 : int AnalyseRaster(GDALDatasetH hDS, const char* dsFileName,
223 : DatasetProperty* psDatasetProperties);
224 :
225 : void CreateVRTSeparate(VRTDatasetH hVRTDS);
226 : void CreateVRTNonSeparate(VRTDatasetH hVRTDS);
227 :
228 : public:
229 : VRTBuilder(const char* pszOutputFilename,
230 : int nInputFiles, const char* const * ppszInputFilenames,
231 : ResolutionStrategy resolutionStrategy,
232 : double we_res, double ns_res,
233 : int bTargetAlignedPixels,
234 : double minX, double minY, double maxX, double maxY,
235 : int bSeparate, int bAllowProjectionDifference,
236 : int bAddAlpha, int bHideNoData,
237 : const char* pszSrcNoData, const char* pszVRTNoData,
238 : const char* pszOutputSRS);
239 :
240 : ~VRTBuilder();
241 :
242 : int Build(GDALProgressFunc pfnProgress, void * pProgressData);
243 : };
244 :
245 :
246 : /************************************************************************/
247 : /* VRTBuilder() */
248 : /************************************************************************/
249 :
250 16 : VRTBuilder::VRTBuilder(const char* pszOutputFilename,
251 : int nInputFiles, const char* const * ppszInputFilenames,
252 : ResolutionStrategy resolutionStrategy,
253 : double we_res, double ns_res,
254 : int bTargetAlignedPixels,
255 : double minX, double minY, double maxX, double maxY,
256 : int bSeparate, int bAllowProjectionDifference,
257 : int bAddAlpha, int bHideNoData,
258 : const char* pszSrcNoData, const char* pszVRTNoData,
259 : const char* pszOutputSRS)
260 : {
261 16 : this->pszOutputFilename = CPLStrdup(pszOutputFilename);
262 16 : this->nInputFiles = nInputFiles;
263 :
264 16 : this->ppszInputFilenames = (char**) CPLMalloc(nInputFiles * sizeof(char*));
265 : int i;
266 61 : for(i=0;i<nInputFiles;i++)
267 : {
268 45 : this->ppszInputFilenames[i] = CPLStrdup(ppszInputFilenames[i]);
269 : }
270 :
271 16 : this->resolutionStrategy = resolutionStrategy;
272 16 : this->we_res = we_res;
273 16 : this->ns_res = ns_res;
274 16 : this->bTargetAlignedPixels = bTargetAlignedPixels;
275 16 : this->minX = minX;
276 16 : this->minY = minY;
277 16 : this->maxX = maxX;
278 16 : this->maxY = maxY;
279 16 : this->bSeparate = bSeparate;
280 16 : this->bAllowProjectionDifference = bAllowProjectionDifference;
281 16 : this->bAddAlpha = bAddAlpha;
282 16 : this->bHideNoData = bHideNoData;
283 16 : this->pszSrcNoData = (pszSrcNoData) ? CPLStrdup(pszSrcNoData) : NULL;
284 16 : this->pszVRTNoData = (pszVRTNoData) ? CPLStrdup(pszVRTNoData) : NULL;
285 16 : this->pszOutputSRS = (pszOutputSRS) ? CPLStrdup(pszOutputSRS) : NULL;
286 :
287 16 : bUserExtent = FALSE;
288 16 : pszProjectionRef = NULL;
289 16 : nBands = 0;
290 16 : pasBandProperties = NULL;
291 16 : bFirst = TRUE;
292 16 : bHasGeoTransform = FALSE;
293 16 : nRasterXSize = 0;
294 16 : nRasterYSize = 0;
295 16 : pasDatasetProperties = NULL;
296 16 : bAllowSrcNoData = TRUE;
297 16 : padfSrcNoData = NULL;
298 16 : nSrcNoDataCount = 0;
299 16 : bAllowVRTNoData = TRUE;
300 16 : padfVRTNoData = NULL;
301 16 : nVRTNoDataCount = 0;
302 16 : bHasRunBuild = FALSE;
303 16 : bHasDatasetMask = FALSE;
304 16 : }
305 :
306 : /************************************************************************/
307 : /* ~VRTBuilder() */
308 : /************************************************************************/
309 :
310 16 : VRTBuilder::~VRTBuilder()
311 : {
312 16 : CPLFree(pszOutputFilename);
313 16 : CPLFree(pszSrcNoData);
314 16 : CPLFree(pszVRTNoData);
315 :
316 : int i;
317 61 : for(i=0;i<nInputFiles;i++)
318 : {
319 45 : CPLFree(ppszInputFilenames[i]);
320 : }
321 16 : CPLFree(ppszInputFilenames);
322 :
323 16 : if (pasDatasetProperties != NULL)
324 : {
325 61 : for(i=0;i<nInputFiles;i++)
326 : {
327 45 : CPLFree(pasDatasetProperties[i].padfNoDataValues);
328 45 : CPLFree(pasDatasetProperties[i].panHasNoData);
329 : }
330 : }
331 16 : CPLFree(pasDatasetProperties);
332 :
333 16 : if (!bSeparate && pasBandProperties != NULL)
334 : {
335 : int j;
336 32 : for(j=0;j<nBands;j++)
337 : {
338 18 : GDALDestroyColorTable(pasBandProperties[j].colorTable);
339 : }
340 : }
341 16 : CPLFree(pasBandProperties);
342 :
343 16 : CPLFree(pszProjectionRef);
344 16 : CPLFree(padfSrcNoData);
345 16 : CPLFree(padfVRTNoData);
346 16 : CPLFree(pszOutputSRS);
347 16 : }
348 :
349 : /************************************************************************/
350 : /* ProjAreEqual() */
351 : /************************************************************************/
352 :
353 29 : static int ProjAreEqual(const char* pszWKT1, const char* pszWKT2)
354 : {
355 : int bRet;
356 : OGRSpatialReferenceH hSRS1, hSRS2;
357 :
358 29 : if (EQUAL(pszWKT1, pszWKT2))
359 28 : return TRUE;
360 :
361 1 : hSRS1 = OSRNewSpatialReference(pszWKT1);
362 1 : hSRS2 = OSRNewSpatialReference(pszWKT2);
363 1 : bRet = hSRS1 != NULL && hSRS2 != NULL && OSRIsSame(hSRS1,hSRS2);
364 1 : if (hSRS1)
365 1 : OSRDestroySpatialReference(hSRS1);
366 1 : if (hSRS2)
367 1 : OSRDestroySpatialReference(hSRS2);
368 1 : return bRet;
369 : }
370 :
371 : /************************************************************************/
372 : /* AnalyseRaster() */
373 : /************************************************************************/
374 :
375 45 : int VRTBuilder::AnalyseRaster( GDALDatasetH hDS, const char* dsFileName,
376 : DatasetProperty* psDatasetProperties)
377 : {
378 45 : char** papszMetadata = GDALGetMetadata( hDS, "SUBDATASETS" );
379 45 : if( CSLCount(papszMetadata) > 0 && GDALGetRasterCount(hDS) == 0 )
380 : {
381 : pasDatasetProperties =
382 : (DatasetProperty*) CPLRealloc(pasDatasetProperties,
383 0 : (nInputFiles+CSLCount(papszMetadata))*sizeof(DatasetProperty));
384 :
385 : ppszInputFilenames = (char**)CPLRealloc(ppszInputFilenames,
386 0 : sizeof(char*) * (nInputFiles+CSLCount(papszMetadata)));
387 0 : int count = 1;
388 : char subdatasetNameKey[256];
389 0 : sprintf(subdatasetNameKey, "SUBDATASET_%d_NAME", count);
390 0 : while(*papszMetadata != NULL)
391 : {
392 0 : if (EQUALN(*papszMetadata, subdatasetNameKey, strlen(subdatasetNameKey)))
393 : {
394 0 : memset(&pasDatasetProperties[nInputFiles], 0, sizeof(DatasetProperty));
395 0 : ppszInputFilenames[nInputFiles++] =
396 0 : CPLStrdup(*papszMetadata+strlen(subdatasetNameKey)+1);
397 0 : count++;
398 0 : sprintf(subdatasetNameKey, "SUBDATASET_%d_NAME", count);
399 : }
400 0 : papszMetadata++;
401 : }
402 0 : return FALSE;
403 : }
404 :
405 45 : const char* proj = GDALGetProjectionRef(hDS);
406 45 : double* padfGeoTransform = psDatasetProperties->adfGeoTransform;
407 45 : int bGotGeoTransform = GDALGetGeoTransform(hDS, padfGeoTransform) == CE_None;
408 45 : if (bSeparate)
409 : {
410 6 : if (bFirst)
411 : {
412 2 : bHasGeoTransform = bGotGeoTransform;
413 2 : if (!bHasGeoTransform)
414 : {
415 1 : if (bUserExtent)
416 : {
417 : CPLError(CE_Warning, CPLE_NotSupported,
418 0 : "User extent ignored by gdalbuildvrt -separate with ungeoreferenced images.");
419 : }
420 1 : if (resolutionStrategy == USER_RESOLUTION)
421 : {
422 : CPLError(CE_Warning, CPLE_NotSupported,
423 0 : "User resolution ignored by gdalbuildvrt -separate with ungeoreferenced images.");
424 : }
425 : }
426 : }
427 4 : else if (bHasGeoTransform != bGotGeoTransform)
428 : {
429 : CPLError(CE_Warning, CPLE_NotSupported,
430 : "gdalbuildvrt -separate cannot stack ungeoreferenced and georeferenced images. Skipping %s",
431 0 : dsFileName);
432 0 : return FALSE;
433 : }
434 4 : else if (!bHasGeoTransform &&
435 : (nRasterXSize != GDALGetRasterXSize(hDS) ||
436 : nRasterYSize != GDALGetRasterYSize(hDS)))
437 : {
438 : CPLError(CE_Warning, CPLE_NotSupported,
439 : "gdalbuildvrt -separate cannot stack ungeoreferenced images that have not the same dimensions. Skipping %s",
440 0 : dsFileName);
441 0 : return FALSE;
442 : }
443 : }
444 : else
445 : {
446 39 : if (!bGotGeoTransform)
447 : {
448 : CPLError(CE_Warning, CPLE_NotSupported,
449 : "gdalbuildvrt does not support ungeoreferenced image. Skipping %s",
450 0 : dsFileName);
451 0 : return FALSE;
452 : }
453 39 : bHasGeoTransform = TRUE;
454 : }
455 :
456 45 : if (bGotGeoTransform)
457 : {
458 86 : if (padfGeoTransform[GEOTRSFRM_ROTATION_PARAM1] != 0 ||
459 43 : padfGeoTransform[GEOTRSFRM_ROTATION_PARAM2] != 0)
460 : {
461 : CPLError(CE_Warning, CPLE_NotSupported,
462 : "gdalbuildvrt does not support rotated geo transforms. Skipping %s",
463 0 : dsFileName);
464 0 : return FALSE;
465 : }
466 43 : if (padfGeoTransform[GEOTRSFRM_NS_RES] >= 0)
467 : {
468 : CPLError(CE_Warning, CPLE_NotSupported,
469 : "gdalbuildvrt does not support positive NS resolution. Skipping %s",
470 0 : dsFileName);
471 0 : return FALSE;
472 : }
473 : }
474 :
475 45 : psDatasetProperties->nRasterXSize = GDALGetRasterXSize(hDS);
476 45 : psDatasetProperties->nRasterYSize = GDALGetRasterYSize(hDS);
477 45 : if (bFirst && bSeparate && !bGotGeoTransform)
478 : {
479 1 : nRasterXSize = GDALGetRasterXSize(hDS);
480 1 : nRasterYSize = GDALGetRasterYSize(hDS);
481 : }
482 :
483 45 : double ds_minX = padfGeoTransform[GEOTRSFRM_TOPLEFT_X];
484 45 : double ds_maxY = padfGeoTransform[GEOTRSFRM_TOPLEFT_Y];
485 : double ds_maxX = ds_minX +
486 : GDALGetRasterXSize(hDS) *
487 45 : padfGeoTransform[GEOTRSFRM_WE_RES];
488 : double ds_minY = ds_maxY +
489 : GDALGetRasterYSize(hDS) *
490 45 : padfGeoTransform[GEOTRSFRM_NS_RES];
491 :
492 : GDALGetBlockSize(GDALGetRasterBand( hDS, 1 ),
493 : &psDatasetProperties->nBlockXSize,
494 45 : &psDatasetProperties->nBlockYSize);
495 :
496 45 : int _nBands = GDALGetRasterCount(hDS);
497 45 : if (_nBands == 0)
498 : {
499 : CPLError(CE_Warning, CPLE_AppDefined,
500 0 : "Skipping %s as it has no bands", dsFileName);
501 0 : return FALSE;
502 : }
503 45 : else if (_nBands > 1 && bSeparate)
504 : {
505 : CPLError(CE_Warning, CPLE_AppDefined, "%s has %d bands. Only the first one will "
506 : "be taken into account in the -separate case",
507 0 : dsFileName, _nBands);
508 0 : _nBands = 1;
509 : }
510 :
511 : /* For the -separate case */
512 45 : psDatasetProperties->firstBandType = GDALGetRasterDataType(GDALGetRasterBand(hDS, 1));
513 :
514 45 : psDatasetProperties->padfNoDataValues = (double*)CPLCalloc(sizeof(double), _nBands);
515 45 : psDatasetProperties->panHasNoData = (int*)CPLCalloc(sizeof(int), _nBands);
516 :
517 45 : psDatasetProperties->bHasDatasetMask = GDALGetMaskFlags(GDALGetRasterBand(hDS, 1)) == GMF_PER_DATASET;
518 45 : if (psDatasetProperties->bHasDatasetMask)
519 1 : bHasDatasetMask = TRUE;
520 : GDALGetBlockSize(GDALGetMaskBand(GDALGetRasterBand( hDS, 1 )),
521 : &psDatasetProperties->nMaskBlockXSize,
522 45 : &psDatasetProperties->nMaskBlockYSize);
523 :
524 : int j;
525 97 : for(j=0;j<_nBands;j++)
526 : {
527 52 : if (nSrcNoDataCount > 0)
528 : {
529 2 : psDatasetProperties->panHasNoData[j] = TRUE;
530 2 : if (j < nSrcNoDataCount)
531 2 : psDatasetProperties->padfNoDataValues[j] = padfSrcNoData[j];
532 : else
533 0 : psDatasetProperties->padfNoDataValues[j] = padfSrcNoData[nSrcNoDataCount - 1];
534 : }
535 : else
536 : {
537 50 : psDatasetProperties->padfNoDataValues[j] =
538 : GDALGetRasterNoDataValue(GDALGetRasterBand(hDS, j+1),
539 50 : &psDatasetProperties->panHasNoData[j]);
540 : }
541 : }
542 :
543 45 : if (bFirst)
544 : {
545 16 : if (proj)
546 16 : pszProjectionRef = CPLStrdup(proj);
547 16 : if (!bUserExtent)
548 : {
549 14 : minX = ds_minX;
550 14 : minY = ds_minY;
551 14 : maxX = ds_maxX;
552 14 : maxY = ds_maxY;
553 : }
554 16 : nBands = _nBands;
555 :
556 16 : if (!bSeparate)
557 : {
558 14 : pasBandProperties = (BandProperty*)CPLMalloc(nBands*sizeof(BandProperty));
559 32 : for(j=0;j<nBands;j++)
560 : {
561 18 : GDALRasterBandH hRasterBand = GDALGetRasterBand( hDS, j+1 );
562 18 : pasBandProperties[j].colorInterpretation =
563 18 : GDALGetRasterColorInterpretation(hRasterBand);
564 18 : pasBandProperties[j].dataType = GDALGetRasterDataType(hRasterBand);
565 18 : if (pasBandProperties[j].colorInterpretation == GCI_PaletteIndex)
566 : {
567 1 : pasBandProperties[j].colorTable =
568 1 : GDALGetRasterColorTable( hRasterBand );
569 1 : if (pasBandProperties[j].colorTable)
570 : {
571 1 : pasBandProperties[j].colorTable =
572 1 : GDALCloneColorTable(pasBandProperties[j].colorTable);
573 : }
574 : }
575 : else
576 17 : pasBandProperties[j].colorTable = 0;
577 :
578 18 : if (nVRTNoDataCount > 0)
579 : {
580 1 : pasBandProperties[j].bHasNoData = TRUE;
581 1 : if (j < nVRTNoDataCount)
582 1 : pasBandProperties[j].noDataValue = padfVRTNoData[j];
583 : else
584 0 : pasBandProperties[j].noDataValue = padfVRTNoData[nVRTNoDataCount - 1];
585 : }
586 : else
587 : {
588 17 : pasBandProperties[j].noDataValue =
589 17 : GDALGetRasterNoDataValue(hRasterBand, &pasBandProperties[j].bHasNoData);
590 : }
591 : }
592 : }
593 : }
594 : else
595 : {
596 29 : if ((proj != NULL && pszProjectionRef == NULL) ||
597 : (proj == NULL && pszProjectionRef != NULL) ||
598 : (proj != NULL && pszProjectionRef != NULL && ProjAreEqual(proj, pszProjectionRef) == FALSE))
599 : {
600 1 : if (!bAllowProjectionDifference)
601 : {
602 : CPLError(CE_Warning, CPLE_NotSupported,
603 : "gdalbuildvrt does not support heterogenous projection. Skipping %s",
604 1 : dsFileName);
605 1 : return FALSE;
606 : }
607 : }
608 28 : if (!bSeparate)
609 : {
610 24 : if (nBands != _nBands)
611 : {
612 : CPLError(CE_Warning, CPLE_NotSupported,
613 : "gdalbuildvrt does not support heterogenous band numbers. Skipping %s",
614 1 : dsFileName);
615 1 : return FALSE;
616 : }
617 48 : for(j=0;j<nBands;j++)
618 : {
619 25 : GDALRasterBandH hRasterBand = GDALGetRasterBand( hDS, j+1 );
620 50 : if (pasBandProperties[j].colorInterpretation != GDALGetRasterColorInterpretation(hRasterBand) ||
621 25 : pasBandProperties[j].dataType != GDALGetRasterDataType(hRasterBand))
622 : {
623 : CPLError(CE_Warning, CPLE_NotSupported,
624 : "gdalbuildvrt does not support heterogenous band characteristics. Skipping %s",
625 0 : dsFileName);
626 0 : return FALSE;
627 : }
628 25 : if (pasBandProperties[j].colorTable)
629 : {
630 1 : GDALColorTableH colorTable = GDALGetRasterColorTable( hRasterBand );
631 1 : int nRefColorEntryCount = GDALGetColorEntryCount(pasBandProperties[j].colorTable);
632 : int i;
633 1 : if (colorTable == NULL ||
634 : GDALGetColorEntryCount(colorTable) != nRefColorEntryCount)
635 : {
636 : CPLError(CE_Warning, CPLE_NotSupported,
637 : "gdalbuildvrt does not support rasters with different color tables (different number of color table entries). Skipping %s",
638 0 : dsFileName);
639 0 : return FALSE;
640 : }
641 :
642 : /* Check that the palette are the same too */
643 : /* We just warn and still process the file. It is not a technical no-go, but the user */
644 : /* should check that the end result is OK for him. */
645 3 : for(i=0;i<nRefColorEntryCount;i++)
646 : {
647 2 : const GDALColorEntry* psEntry = GDALGetColorEntry(colorTable, i);
648 2 : const GDALColorEntry* psEntryRef = GDALGetColorEntry(pasBandProperties[j].colorTable, i);
649 2 : if (psEntry->c1 != psEntryRef->c1 || psEntry->c2 != psEntryRef->c2 ||
650 : psEntry->c3 != psEntryRef->c3 || psEntry->c4 != psEntryRef->c4)
651 : {
652 : static int bFirstWarningPCT = TRUE;
653 0 : if (bFirstWarningPCT)
654 : CPLError(CE_Warning, CPLE_NotSupported,
655 : "%s has different values than the first raster for some entries in the color table.\n"
656 : "The end result might produce weird colors.\n"
657 : "You're advised to preprocess your rasters with other tools, such as pct2rgb.py or gdal_translate -expand RGB\n"
658 0 : "to operate gdalbuildvrt on RGB rasters instead", dsFileName);
659 : else
660 : CPLError(CE_Warning, CPLE_NotSupported,
661 : "%s has different values than the first raster for some entries in the color table.",
662 0 : dsFileName);
663 0 : bFirstWarningPCT = FALSE;
664 0 : break;
665 : }
666 : }
667 : }
668 : }
669 :
670 : }
671 27 : if (!bUserExtent)
672 : {
673 24 : if (ds_minX < minX) minX = ds_minX;
674 24 : if (ds_minY < minY) minY = ds_minY;
675 24 : if (ds_maxX > maxX) maxX = ds_maxX;
676 24 : if (ds_maxY > maxY) maxY = ds_maxY;
677 : }
678 : }
679 :
680 43 : if (resolutionStrategy == AVERAGE_RESOLUTION)
681 : {
682 37 : we_res += padfGeoTransform[GEOTRSFRM_WE_RES];
683 37 : ns_res += padfGeoTransform[GEOTRSFRM_NS_RES];
684 : }
685 6 : else if (resolutionStrategy != USER_RESOLUTION)
686 : {
687 0 : if (bFirst)
688 : {
689 0 : we_res = padfGeoTransform[GEOTRSFRM_WE_RES];
690 0 : ns_res = padfGeoTransform[GEOTRSFRM_NS_RES];
691 : }
692 0 : else if (resolutionStrategy == HIGHEST_RESOLUTION)
693 : {
694 0 : we_res = MIN(we_res, padfGeoTransform[GEOTRSFRM_WE_RES]);
695 : /* Yes : as ns_res is negative, the highest resolution is the max value */
696 0 : ns_res = MAX(ns_res, padfGeoTransform[GEOTRSFRM_NS_RES]);
697 : }
698 : else
699 : {
700 0 : we_res = MAX(we_res, padfGeoTransform[GEOTRSFRM_WE_RES]);
701 : /* Yes : as ns_res is negative, the lowest resolution is the min value */
702 0 : ns_res = MIN(ns_res, padfGeoTransform[GEOTRSFRM_NS_RES]);
703 : }
704 : }
705 :
706 43 : return TRUE;
707 : }
708 :
709 : /************************************************************************/
710 : /* CreateVRTSeparate() */
711 : /************************************************************************/
712 :
713 2 : void VRTBuilder::CreateVRTSeparate(VRTDatasetH hVRTDS)
714 : {
715 : int i;
716 2 : int iBand = 1;
717 8 : for(i=0;i<nInputFiles;i++)
718 : {
719 6 : DatasetProperty* psDatasetProperties = &pasDatasetProperties[i];
720 :
721 6 : if (psDatasetProperties->isFileOK == FALSE)
722 0 : continue;
723 :
724 : int nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
725 : nDstXOff, nDstYOff, nDstXSize, nDstYSize;
726 6 : if (bHasGeoTransform)
727 : {
728 4 : if ( ! GetSrcDstWin(psDatasetProperties,
729 : we_res, ns_res, minX, minY, maxX, maxY,
730 : &nSrcXOff, &nSrcYOff, &nSrcXSize, &nSrcYSize,
731 : &nDstXOff, &nDstYOff, &nDstXSize, &nDstYSize) )
732 0 : continue;
733 : }
734 : else
735 : {
736 2 : nSrcXOff = nSrcYOff = nDstXOff = nDstYOff = 0;
737 2 : nSrcXSize = nDstXSize = nRasterXSize;
738 2 : nSrcYSize = nDstYSize = nRasterYSize;
739 : }
740 :
741 6 : const char* dsFileName = ppszInputFilenames[i];
742 :
743 6 : GDALAddBand(hVRTDS, psDatasetProperties->firstBandType, NULL);
744 :
745 : GDALProxyPoolDatasetH hProxyDS =
746 : GDALProxyPoolDatasetCreate(dsFileName,
747 : psDatasetProperties->nRasterXSize,
748 : psDatasetProperties->nRasterYSize,
749 : GA_ReadOnly, TRUE, pszProjectionRef,
750 6 : psDatasetProperties->adfGeoTransform);
751 : GDALProxyPoolDatasetAddSrcBandDescription(hProxyDS,
752 : psDatasetProperties->firstBandType,
753 : psDatasetProperties->nBlockXSize,
754 6 : psDatasetProperties->nBlockYSize);
755 :
756 : VRTSourcedRasterBandH hVRTBand =
757 6 : (VRTSourcedRasterBandH)GDALGetRasterBand(hVRTDS, iBand);
758 :
759 6 : if (bHideNoData)
760 0 : GDALSetMetadataItem(hVRTBand,"HideNoDataValue","1",NULL);
761 :
762 6 : if (bAllowSrcNoData && psDatasetProperties->panHasNoData[0])
763 : {
764 0 : GDALSetRasterNoDataValue(hVRTBand, psDatasetProperties->padfNoDataValues[0]);
765 : VRTAddComplexSource(hVRTBand, GDALGetRasterBand((GDALDatasetH)hProxyDS, 1),
766 : nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
767 : nDstXOff, nDstYOff, nDstXSize, nDstYSize,
768 0 : 0, 1, psDatasetProperties->padfNoDataValues[0]);
769 : }
770 : else
771 : /* Place the raster band at the right position in the VRT */
772 : VRTAddSimpleSource(hVRTBand, GDALGetRasterBand((GDALDatasetH)hProxyDS, 1),
773 : nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
774 : nDstXOff, nDstYOff, nDstXSize, nDstYSize,
775 6 : "near", VRT_NODATA_UNSET);
776 :
777 6 : GDALDereferenceDataset(hProxyDS);
778 :
779 6 : iBand ++;
780 : }
781 2 : }
782 :
783 : /************************************************************************/
784 : /* CreateVRTNonSeparate() */
785 : /************************************************************************/
786 :
787 14 : void VRTBuilder::CreateVRTNonSeparate(VRTDatasetH hVRTDS)
788 : {
789 : int i, j;
790 :
791 32 : for(j=0;j<nBands;j++)
792 : {
793 : GDALRasterBandH hBand;
794 18 : GDALAddBand(hVRTDS, pasBandProperties[j].dataType, NULL);
795 18 : hBand = GDALGetRasterBand(hVRTDS, j+1);
796 18 : GDALSetRasterColorInterpretation(hBand, pasBandProperties[j].colorInterpretation);
797 18 : if (pasBandProperties[j].colorInterpretation == GCI_PaletteIndex)
798 : {
799 1 : GDALSetRasterColorTable(hBand, pasBandProperties[j].colorTable);
800 : }
801 18 : if (bAllowVRTNoData && pasBandProperties[j].bHasNoData)
802 4 : GDALSetRasterNoDataValue(hBand, pasBandProperties[j].noDataValue);
803 18 : if ( bHideNoData )
804 0 : GDALSetMetadataItem(hBand,"HideNoDataValue","1",NULL);
805 : }
806 :
807 14 : VRTSourcedRasterBand* hMaskVRTBand = NULL;
808 14 : if (bAddAlpha)
809 : {
810 : GDALRasterBandH hBand;
811 0 : GDALAddBand(hVRTDS, GDT_Byte, NULL);
812 0 : hBand = GDALGetRasterBand(hVRTDS, nBands + 1);
813 0 : GDALSetRasterColorInterpretation(hBand, GCI_AlphaBand);
814 : }
815 14 : else if (bHasDatasetMask)
816 : {
817 1 : GDALCreateDatasetMaskBand(hVRTDS, GMF_PER_DATASET);
818 1 : hMaskVRTBand = (VRTSourcedRasterBand*)GDALGetMaskBand(GDALGetRasterBand(hVRTDS, 1));
819 : }
820 :
821 :
822 53 : for(i=0;i<nInputFiles;i++)
823 : {
824 39 : DatasetProperty* psDatasetProperties = &pasDatasetProperties[i];
825 :
826 39 : if (psDatasetProperties->isFileOK == FALSE)
827 2 : continue;
828 :
829 : int nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
830 : nDstXOff, nDstYOff, nDstXSize, nDstYSize;
831 37 : if ( ! GetSrcDstWin(psDatasetProperties,
832 : we_res, ns_res, minX, minY, maxX, maxY,
833 : &nSrcXOff, &nSrcYOff, &nSrcXSize, &nSrcYSize,
834 : &nDstXOff, &nDstYOff, &nDstXSize, &nDstYSize) )
835 0 : continue;
836 :
837 37 : const char* dsFileName = ppszInputFilenames[i];
838 :
839 : GDALProxyPoolDatasetH hProxyDS =
840 : GDALProxyPoolDatasetCreate(dsFileName,
841 : psDatasetProperties->nRasterXSize,
842 : psDatasetProperties->nRasterYSize,
843 : GA_ReadOnly, TRUE, pszProjectionRef,
844 37 : psDatasetProperties->adfGeoTransform);
845 :
846 80 : for(j=0;j<nBands;j++)
847 : {
848 : GDALProxyPoolDatasetAddSrcBandDescription(hProxyDS,
849 43 : pasBandProperties[j].dataType,
850 : psDatasetProperties->nBlockXSize,
851 86 : psDatasetProperties->nBlockYSize);
852 : }
853 37 : if (bHasDatasetMask && !bAddAlpha)
854 : {
855 : ((GDALProxyPoolRasterBand*)((GDALProxyPoolDataset*)hProxyDS)->GetRasterBand(1))->
856 : AddSrcMaskBandDescription (GDT_Byte,
857 : psDatasetProperties->nMaskBlockXSize,
858 1 : psDatasetProperties->nMaskBlockYSize);
859 : }
860 :
861 80 : for(j=0;j<nBands;j++)
862 : {
863 : VRTSourcedRasterBandH hVRTBand =
864 43 : (VRTSourcedRasterBandH)GDALGetRasterBand(hVRTDS, j + 1);
865 :
866 : /* Place the raster band at the right position in the VRT */
867 51 : if (bAllowSrcNoData && psDatasetProperties->panHasNoData[j])
868 : VRTAddComplexSource(hVRTBand, GDALGetRasterBand((GDALDatasetH)hProxyDS, j + 1),
869 : nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
870 : nDstXOff, nDstYOff, nDstXSize, nDstYSize,
871 8 : 0, 1, psDatasetProperties->padfNoDataValues[j]);
872 : else
873 : VRTAddSimpleSource(hVRTBand, GDALGetRasterBand((GDALDatasetH)hProxyDS, j + 1),
874 : nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
875 : nDstXOff, nDstYOff, nDstXSize, nDstYSize,
876 35 : "near", VRT_NODATA_UNSET);
877 : }
878 :
879 37 : if (bAddAlpha)
880 : {
881 : VRTSourcedRasterBandH hVRTBand =
882 0 : (VRTSourcedRasterBandH)GDALGetRasterBand(hVRTDS, nBands + 1);
883 : /* Little trick : we use an offset of 255 and a scaling of 0, so that in areas covered */
884 : /* by the source, the value of the alpha band will be 255, otherwise it will be 0 */
885 : VRTAddComplexSource(hVRTBand, GDALGetRasterBand((GDALDatasetH)hProxyDS, 1),
886 : nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
887 : nDstXOff, nDstYOff, nDstXSize, nDstYSize,
888 0 : 255, 0, VRT_NODATA_UNSET);
889 : }
890 37 : else if (bHasDatasetMask)
891 : {
892 : hMaskVRTBand->AddMaskBandSource((GDALRasterBand*)GDALGetRasterBand((GDALDatasetH)hProxyDS, 1),
893 : nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
894 1 : nDstXOff, nDstYOff, nDstXSize, nDstYSize);
895 : }
896 :
897 37 : GDALDereferenceDataset(hProxyDS);
898 : }
899 14 : }
900 :
901 : /************************************************************************/
902 : /* Build() */
903 : /************************************************************************/
904 :
905 16 : int VRTBuilder::Build(GDALProgressFunc pfnProgress, void * pProgressData)
906 : {
907 : int i;
908 :
909 16 : if (bHasRunBuild)
910 0 : return CE_Failure;
911 16 : bHasRunBuild = TRUE;
912 :
913 16 : if( pfnProgress == NULL )
914 0 : pfnProgress = GDALDummyProgress;
915 :
916 16 : bUserExtent = (minX != 0 || minY != 0 || maxX != 0 || maxY != 0);
917 16 : if (bUserExtent)
918 : {
919 2 : if (minX >= maxX || minY >= maxY )
920 : {
921 0 : CPLError(CE_Failure, CPLE_IllegalArg, "Invalid user extent");
922 0 : return CE_Failure;
923 : }
924 : }
925 :
926 16 : if (resolutionStrategy == USER_RESOLUTION)
927 : {
928 3 : if (we_res <= 0 || ns_res <= 0)
929 : {
930 0 : CPLError(CE_Failure, CPLE_IllegalArg, "Invalid user resolution");
931 0 : return CE_Failure;
932 : }
933 :
934 : /* We work with negative north-south resolution in all the following code */
935 3 : ns_res = -ns_res;
936 : }
937 : else
938 : {
939 13 : we_res = ns_res = 0;
940 : }
941 :
942 : pasDatasetProperties =
943 16 : (DatasetProperty*) CPLCalloc(nInputFiles, sizeof(DatasetProperty));
944 :
945 16 : if (pszSrcNoData != NULL)
946 : {
947 1 : if (EQUAL(pszSrcNoData, "none"))
948 : {
949 0 : bAllowSrcNoData = FALSE;
950 : }
951 : else
952 : {
953 1 : char **papszTokens = CSLTokenizeString( pszSrcNoData );
954 1 : nSrcNoDataCount = CSLCount(papszTokens);
955 1 : padfSrcNoData = (double *) CPLMalloc(sizeof(double) * nSrcNoDataCount);
956 2 : for(i=0;i<nSrcNoDataCount;i++)
957 1 : padfSrcNoData[i] = CPLAtofM(papszTokens[i]);
958 1 : CSLDestroy(papszTokens);
959 : }
960 : }
961 :
962 16 : if (pszVRTNoData != NULL)
963 : {
964 1 : if (EQUAL(pszVRTNoData, "none"))
965 : {
966 0 : bAllowVRTNoData = FALSE;
967 : }
968 : else
969 : {
970 1 : char **papszTokens = CSLTokenizeString( pszVRTNoData );
971 1 : nVRTNoDataCount = CSLCount(papszTokens);
972 1 : padfVRTNoData = (double *) CPLMalloc(sizeof(double) * nVRTNoDataCount);
973 2 : for(i=0;i<nVRTNoDataCount;i++)
974 1 : padfVRTNoData[i] = CPLAtofM(papszTokens[i]);
975 1 : CSLDestroy(papszTokens);
976 : }
977 : }
978 :
979 16 : int nCountValid = 0;
980 61 : for(i=0;i<nInputFiles;i++)
981 : {
982 45 : const char* dsFileName = ppszInputFilenames[i];
983 :
984 45 : if (!pfnProgress( 1.0 * (i+1) / nInputFiles, NULL, pProgressData))
985 : {
986 0 : return CE_Failure;
987 : }
988 :
989 45 : GDALDatasetH hDS = GDALOpen(ppszInputFilenames[i], GA_ReadOnly );
990 45 : pasDatasetProperties[i].isFileOK = FALSE;
991 :
992 45 : if (hDS)
993 : {
994 45 : if (AnalyseRaster( hDS, dsFileName, &pasDatasetProperties[i] ))
995 : {
996 43 : pasDatasetProperties[i].isFileOK = TRUE;
997 43 : nCountValid ++;
998 43 : bFirst = FALSE;
999 : }
1000 45 : GDALClose(hDS);
1001 : }
1002 : else
1003 : {
1004 : CPLError(CE_Warning, CPLE_AppDefined,
1005 0 : "Can't open %s. Skipping it", dsFileName);
1006 : }
1007 : }
1008 :
1009 16 : if (nCountValid == 0)
1010 0 : return CE_None;
1011 :
1012 16 : if (bHasGeoTransform)
1013 : {
1014 15 : if (resolutionStrategy == AVERAGE_RESOLUTION)
1015 : {
1016 12 : we_res /= nCountValid;
1017 12 : ns_res /= nCountValid;
1018 : }
1019 :
1020 15 : if ( bTargetAlignedPixels )
1021 : {
1022 1 : minX = floor(minX / we_res) * we_res;
1023 1 : maxX = ceil(maxX / we_res) * we_res;
1024 1 : minY = floor(minY / -ns_res) * -ns_res;
1025 1 : maxY = ceil(maxY / -ns_res) * -ns_res;
1026 : }
1027 :
1028 15 : nRasterXSize = (int)(0.5 + (maxX - minX) / we_res);
1029 15 : nRasterYSize = (int)(0.5 + (maxY - minY) / -ns_res);
1030 : }
1031 :
1032 16 : if (nRasterXSize == 0 || nRasterYSize == 0)
1033 : {
1034 : CPLError(CE_Failure, CPLE_AppDefined,
1035 0 : "Computed VRT dimension is invalid. You've probably specified unappropriate resolution.");
1036 0 : return CE_Failure;
1037 : }
1038 :
1039 16 : VRTDatasetH hVRTDS = VRTCreate(nRasterXSize, nRasterYSize);
1040 16 : GDALSetDescription(hVRTDS, pszOutputFilename);
1041 :
1042 16 : if( pszOutputSRS )
1043 : {
1044 1 : GDALSetProjection(hVRTDS, pszOutputSRS);
1045 : }
1046 15 : else if (pszProjectionRef)
1047 : {
1048 15 : GDALSetProjection(hVRTDS, pszProjectionRef);
1049 : }
1050 :
1051 16 : if (bHasGeoTransform)
1052 : {
1053 : double adfGeoTransform[6];
1054 15 : adfGeoTransform[GEOTRSFRM_TOPLEFT_X] = minX;
1055 15 : adfGeoTransform[GEOTRSFRM_WE_RES] = we_res;
1056 15 : adfGeoTransform[GEOTRSFRM_ROTATION_PARAM1] = 0;
1057 15 : adfGeoTransform[GEOTRSFRM_TOPLEFT_Y] = maxY;
1058 15 : adfGeoTransform[GEOTRSFRM_ROTATION_PARAM2] = 0;
1059 15 : adfGeoTransform[GEOTRSFRM_NS_RES] = ns_res;
1060 15 : GDALSetGeoTransform(hVRTDS, adfGeoTransform);
1061 : }
1062 :
1063 16 : if (bSeparate)
1064 : {
1065 2 : CreateVRTSeparate(hVRTDS);
1066 : }
1067 : else
1068 : {
1069 14 : CreateVRTNonSeparate(hVRTDS);
1070 : }
1071 :
1072 16 : GDALClose(hVRTDS);
1073 :
1074 16 : return CE_None;
1075 : }
1076 :
1077 : /************************************************************************/
1078 : /* add_file_to_list() */
1079 : /************************************************************************/
1080 :
1081 43 : static void add_file_to_list(const char* filename, const char* tile_index,
1082 : int* pnInputFiles, char*** pppszInputFilenames)
1083 : {
1084 :
1085 43 : int nInputFiles = *pnInputFiles;
1086 43 : char** ppszInputFilenames = *pppszInputFilenames;
1087 :
1088 43 : if (EQUAL(CPLGetExtension(filename), "SHP"))
1089 : {
1090 : #ifndef OGR_ENABLED
1091 : CPLError(CE_Failure, CPLE_AppDefined, "OGR support needed to read tileindex");
1092 : *pnInputFiles = 0;
1093 : *pppszInputFilenames = NULL;
1094 : #else
1095 : OGRDataSourceH hDS;
1096 : OGRLayerH hLayer;
1097 : OGRFeatureDefnH hFDefn;
1098 : int j, ti_field;
1099 :
1100 1 : OGRRegisterAll();
1101 :
1102 : /* Handle GDALTIndex Shapefile as a special case */
1103 1 : hDS = OGROpen( filename, FALSE, NULL );
1104 1 : if( hDS == NULL )
1105 : {
1106 : fprintf( stderr, "Unable to open shapefile `%s'.\n",
1107 0 : filename );
1108 0 : exit(2);
1109 : }
1110 :
1111 1 : hLayer = OGR_DS_GetLayer(hDS, 0);
1112 :
1113 1 : hFDefn = OGR_L_GetLayerDefn(hLayer);
1114 :
1115 1 : for( ti_field = 0; ti_field < OGR_FD_GetFieldCount(hFDefn); ti_field++ )
1116 : {
1117 1 : OGRFieldDefnH hFieldDefn = OGR_FD_GetFieldDefn( hFDefn, ti_field );
1118 1 : const char* pszName = OGR_Fld_GetNameRef(hFieldDefn);
1119 :
1120 1 : if (strcmp(pszName, "LOCATION") == 0 && strcmp("LOCATION", tile_index) != 0 )
1121 : {
1122 : fprintf( stderr, "This shapefile seems to be a tile index of "
1123 0 : "OGR features and not GDAL products.\n");
1124 : }
1125 1 : if( strcmp(pszName, tile_index) == 0 )
1126 1 : break;
1127 : }
1128 :
1129 1 : if( ti_field == OGR_FD_GetFieldCount(hFDefn) )
1130 : {
1131 : fprintf( stderr, "Unable to find field `%s' in DBF file `%s'.\n",
1132 0 : tile_index, filename );
1133 0 : return;
1134 : }
1135 :
1136 : /* Load in memory existing file names in SHP */
1137 1 : int nTileIndexFiles = OGR_L_GetFeatureCount(hLayer, TRUE);
1138 1 : if (nTileIndexFiles == 0)
1139 : {
1140 0 : fprintf( stderr, "Tile index %s is empty. Skipping it.\n", filename);
1141 0 : return;
1142 : }
1143 :
1144 : ppszInputFilenames = (char**)CPLRealloc(ppszInputFilenames,
1145 1 : sizeof(char*) * (nInputFiles+nTileIndexFiles));
1146 5 : for(j=0;j<nTileIndexFiles;j++)
1147 : {
1148 4 : OGRFeatureH hFeat = OGR_L_GetNextFeature(hLayer);
1149 8 : ppszInputFilenames[nInputFiles++] =
1150 4 : CPLStrdup(OGR_F_GetFieldAsString(hFeat, ti_field ));
1151 4 : OGR_F_Destroy(hFeat);
1152 : }
1153 :
1154 1 : OGR_DS_Destroy( hDS );
1155 : #endif
1156 : }
1157 : else
1158 : {
1159 : ppszInputFilenames = (char**)CPLRealloc(ppszInputFilenames,
1160 42 : sizeof(char*) * (nInputFiles+1));
1161 42 : ppszInputFilenames[nInputFiles++] = CPLStrdup(filename);
1162 : }
1163 :
1164 43 : *pnInputFiles = nInputFiles;
1165 43 : *pppszInputFilenames = ppszInputFilenames;
1166 : }
1167 :
1168 : /************************************************************************/
1169 : /* main() */
1170 : /************************************************************************/
1171 :
1172 18 : int main( int nArgc, char ** papszArgv )
1173 :
1174 : {
1175 18 : const char *tile_index = "location";
1176 18 : const char *resolution = NULL;
1177 18 : int nInputFiles = 0;
1178 18 : char ** ppszInputFilenames = NULL;
1179 18 : const char * pszOutputFilename = NULL;
1180 : int i, iArg;
1181 18 : int bSeparate = FALSE;
1182 18 : int bAllowProjectionDifference = FALSE;
1183 18 : int bQuiet = FALSE;
1184 18 : GDALProgressFunc pfnProgress = NULL;
1185 18 : double we_res = 0, ns_res = 0;
1186 18 : int bTargetAlignedPixels = FALSE;
1187 18 : double xmin = 0, ymin = 0, xmax = 0, ymax = 0;
1188 18 : int bAddAlpha = FALSE;
1189 18 : int bForceOverwrite = FALSE;
1190 18 : int bHideNoData = FALSE;
1191 18 : const char* pszSrcNoData = NULL;
1192 18 : const char* pszVRTNoData = NULL;
1193 18 : char* pszOutputSRS = NULL;
1194 :
1195 : /* Check strict compilation and runtime library version as we use C++ API */
1196 18 : if (! GDAL_CHECK_VERSION(papszArgv[0]))
1197 0 : exit(1);
1198 :
1199 18 : GDALAllRegister();
1200 :
1201 18 : nArgc = GDALGeneralCmdLineProcessor( nArgc, &papszArgv, 0 );
1202 18 : if( nArgc < 1 )
1203 0 : exit( -nArgc );
1204 :
1205 : /* -------------------------------------------------------------------- */
1206 : /* Parse commandline. */
1207 : /* -------------------------------------------------------------------- */
1208 86 : for( iArg = 1; iArg < nArgc; iArg++ )
1209 : {
1210 69 : if( EQUAL(papszArgv[iArg], "--utility_version") )
1211 : {
1212 : printf("%s was compiled against GDAL %s and is running against GDAL %s\n",
1213 1 : papszArgv[0], GDAL_RELEASE_NAME, GDALVersionInfo("RELEASE_NAME"));
1214 1 : return 0;
1215 : }
1216 68 : else if( EQUAL(papszArgv[iArg],"-tileindex") &&
1217 : iArg + 1 < nArgc)
1218 : {
1219 0 : tile_index = papszArgv[++iArg];
1220 : }
1221 68 : else if( EQUAL(papszArgv[iArg],"-resolution") &&
1222 : iArg + 1 < nArgc)
1223 : {
1224 0 : resolution = papszArgv[++iArg];
1225 : }
1226 69 : else if( EQUAL(papszArgv[iArg],"-input_file_list") &&
1227 : iArg + 1 < nArgc)
1228 : {
1229 1 : const char* input_file_list = papszArgv[++iArg];
1230 1 : FILE* f = VSIFOpen(input_file_list, "r");
1231 1 : if (f)
1232 : {
1233 4 : while(1)
1234 : {
1235 5 : const char* filename = CPLReadLine(f);
1236 5 : if (filename == NULL)
1237 : break;
1238 : add_file_to_list(filename, tile_index,
1239 4 : &nInputFiles, &ppszInputFilenames);
1240 : }
1241 1 : VSIFClose(f);
1242 : }
1243 : }
1244 67 : else if ( EQUAL(papszArgv[iArg],"-separate") )
1245 : {
1246 2 : bSeparate = TRUE;
1247 : }
1248 65 : else if ( EQUAL(papszArgv[iArg],"-allow_projection_difference") )
1249 : {
1250 0 : bAllowProjectionDifference = TRUE;
1251 : }
1252 : /* Alternate syntax for output file */
1253 65 : else if( EQUAL(papszArgv[iArg],"-o") &&
1254 : iArg + 1 < nArgc)
1255 : {
1256 0 : pszOutputFilename = papszArgv[++iArg];
1257 : }
1258 65 : else if ( EQUAL(papszArgv[iArg],"-q") || EQUAL(papszArgv[iArg],"-quiet") )
1259 : {
1260 0 : bQuiet = TRUE;
1261 : }
1262 68 : else if ( EQUAL(papszArgv[iArg],"-tr") && iArg + 2 < nArgc)
1263 : {
1264 3 : we_res = CPLAtofM(papszArgv[++iArg]);
1265 3 : ns_res = CPLAtofM(papszArgv[++iArg]);
1266 : }
1267 62 : else if( EQUAL(papszArgv[iArg],"-tap") )
1268 : {
1269 2 : bTargetAlignedPixels = TRUE;
1270 : }
1271 62 : else if ( EQUAL(papszArgv[iArg],"-te") && iArg + 4 < nArgc)
1272 : {
1273 2 : xmin = CPLAtofM(papszArgv[++iArg]);
1274 2 : ymin = CPLAtofM(papszArgv[++iArg]);
1275 2 : xmax = CPLAtofM(papszArgv[++iArg]);
1276 2 : ymax = CPLAtofM(papszArgv[++iArg]);
1277 : }
1278 58 : else if ( EQUAL(papszArgv[iArg],"-addalpha") )
1279 : {
1280 0 : bAddAlpha = TRUE;
1281 : }
1282 58 : else if ( EQUAL(papszArgv[iArg],"-hidenodata") )
1283 : {
1284 0 : bHideNoData = TRUE;
1285 : }
1286 58 : else if ( EQUAL(papszArgv[iArg],"-overwrite") )
1287 : {
1288 0 : bForceOverwrite = TRUE;
1289 : }
1290 59 : else if ( EQUAL(papszArgv[iArg],"-srcnodata") && iArg + 1 < nArgc)
1291 : {
1292 1 : pszSrcNoData = papszArgv[++iArg];
1293 : }
1294 57 : else if ( EQUAL(papszArgv[iArg],"-vrtnodata") && iArg + 1 < nArgc)
1295 : {
1296 0 : pszVRTNoData = papszArgv[++iArg];
1297 : }
1298 58 : else if( EQUAL(papszArgv[iArg],"-a_srs") && iArg + 1 < nArgc)
1299 : {
1300 1 : OGRSpatialReferenceH hOutputSRS = OSRNewSpatialReference(NULL);
1301 :
1302 1 : if( OSRSetFromUserInput( hOutputSRS, papszArgv[iArg+1] ) != OGRERR_NONE )
1303 : {
1304 : fprintf( stderr, "Failed to process SRS definition: %s\n",
1305 0 : papszArgv[iArg+1] );
1306 0 : GDALDestroyDriverManager();
1307 0 : exit( 1 );
1308 : }
1309 :
1310 1 : OSRExportToWkt( hOutputSRS, &pszOutputSRS );
1311 1 : OSRDestroySpatialReference( hOutputSRS );
1312 1 : iArg++;
1313 : }
1314 56 : else if ( papszArgv[iArg][0] == '-' )
1315 : {
1316 0 : printf("Unrecognized option : %s\n", papszArgv[iArg]);
1317 0 : exit(-1);
1318 : }
1319 56 : else if( pszOutputFilename == NULL )
1320 : {
1321 17 : pszOutputFilename = papszArgv[iArg];
1322 : }
1323 : else
1324 : {
1325 39 : add_file_to_list(papszArgv[iArg], tile_index,
1326 78 : &nInputFiles, &ppszInputFilenames);
1327 : }
1328 : }
1329 :
1330 17 : if( pszOutputFilename == NULL || nInputFiles == 0 )
1331 0 : Usage();
1332 :
1333 17 : if (!bQuiet)
1334 17 : pfnProgress = GDALTermProgress;
1335 :
1336 : /* Avoid overwriting a non VRT dataset if the user did not put the */
1337 : /* filenames in the right order */
1338 : VSIStatBuf sBuf;
1339 17 : if (!bForceOverwrite)
1340 : {
1341 17 : int bExists = (VSIStat(pszOutputFilename, &sBuf) == 0);
1342 17 : if (bExists)
1343 : {
1344 6 : GDALDriverH hDriver = GDALIdentifyDriver( pszOutputFilename, NULL );
1345 6 : if (hDriver && !EQUAL(GDALGetDriverShortName(hDriver), "VRT"))
1346 : {
1347 : fprintf(stderr,
1348 : "'%s' is an existing GDAL dataset managed by %s driver.\n"
1349 : "There is an high chance you did not put filenames in the right order.\n"
1350 : "If you want to overwrite %s, add -overwrite option to the command line.\n\n",
1351 0 : pszOutputFilename, GDALGetDriverShortName(hDriver), pszOutputFilename);
1352 0 : Usage();
1353 : }
1354 : }
1355 : }
1356 :
1357 17 : if (we_res != 0 && ns_res != 0 &&
1358 : resolution != NULL && !EQUAL(resolution, "user"))
1359 : {
1360 0 : fprintf(stderr, "-tr option is not compatible with -resolution %s\n", resolution);
1361 0 : Usage();
1362 : }
1363 :
1364 17 : if (bTargetAlignedPixels && we_res == 0 && ns_res == 0)
1365 : {
1366 1 : fprintf( stderr, "-tap option cannot be used without using -tr\n");
1367 1 : Usage();
1368 : }
1369 :
1370 16 : if (bAddAlpha && bSeparate)
1371 : {
1372 0 : fprintf(stderr, "-addalpha option is not compatible with -separate\n");
1373 0 : Usage();
1374 : }
1375 :
1376 16 : ResolutionStrategy eStrategy = AVERAGE_RESOLUTION;
1377 32 : if ( resolution == NULL || EQUAL(resolution, "user") )
1378 : {
1379 19 : if ( we_res != 0 || ns_res != 0)
1380 3 : eStrategy = USER_RESOLUTION;
1381 13 : else if ( resolution != NULL && EQUAL(resolution, "user") )
1382 : {
1383 0 : fprintf(stderr, "-tr option must be used with -resolution user\n");
1384 0 : Usage();
1385 : }
1386 : }
1387 0 : else if ( EQUAL(resolution, "average") )
1388 0 : eStrategy = AVERAGE_RESOLUTION;
1389 0 : else if ( EQUAL(resolution, "highest") )
1390 0 : eStrategy = HIGHEST_RESOLUTION;
1391 0 : else if ( EQUAL(resolution, "lowest") )
1392 0 : eStrategy = LOWEST_RESOLUTION;
1393 : else
1394 : {
1395 0 : fprintf(stderr, "invalid value (%s) for -resolution\n", resolution);
1396 0 : Usage();
1397 : }
1398 :
1399 : /* If -srcnodata is specified, use it as the -vrtnodata if the latter is not */
1400 : /* specified */
1401 16 : if (pszSrcNoData != NULL && pszVRTNoData == NULL)
1402 1 : pszVRTNoData = pszSrcNoData;
1403 :
1404 : VRTBuilder oBuilder(pszOutputFilename, nInputFiles, ppszInputFilenames,
1405 : eStrategy, we_res, ns_res, bTargetAlignedPixels, xmin, ymin, xmax, ymax,
1406 : bSeparate, bAllowProjectionDifference, bAddAlpha, bHideNoData,
1407 16 : pszSrcNoData, pszVRTNoData, pszOutputSRS);
1408 :
1409 16 : oBuilder.Build(pfnProgress, NULL);
1410 :
1411 61 : for(i=0;i<nInputFiles;i++)
1412 : {
1413 45 : CPLFree(ppszInputFilenames[i]);
1414 : }
1415 16 : CPLFree(ppszInputFilenames);
1416 16 : CPLFree(pszOutputSRS);
1417 :
1418 16 : CSLDestroy( papszArgv );
1419 16 : GDALDumpOpenDatasets( stderr );
1420 16 : GDALDestroyDriverManager();
1421 : #ifdef OGR_ENABLED
1422 16 : OGRCleanupAll();
1423 : #endif
1424 :
1425 16 : return 0;
1426 : }
|
