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