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