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