1 : /******************************************************************************
2 : * $Id: gdalbuildvrt.cpp 18306 2009-12-15 18:57:11Z rouault $
3 : *
4 : * Project: GDAL Utilities
5 : * Purpose: Commandline application to build VRT datasets from raster products or content of SHP tile index
6 : * Author: Even Rouault, even.rouault at mines-paris.org
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2007, Even Rouault
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "gdal_proxy.h"
31 : #include "cpl_string.h"
32 : #include "vrt/gdal_vrt.h"
33 :
34 : #ifdef OGR_ENABLED
35 : #include "ogr_api.h"
36 : #endif
37 :
38 : CPL_CVSID("$Id: gdalbuildvrt.cpp 18306 2009-12-15 18:57:11Z rouault $");
39 :
40 : #define GEOTRSFRM_TOPLEFT_X 0
41 : #define GEOTRSFRM_WE_RES 1
42 : #define GEOTRSFRM_ROTATION_PARAM1 2
43 : #define GEOTRSFRM_TOPLEFT_Y 3
44 : #define GEOTRSFRM_ROTATION_PARAM2 4
45 : #define GEOTRSFRM_NS_RES 5
46 :
47 : typedef enum
48 : {
49 : LOWEST_RESOLUTION,
50 : HIGHEST_RESOLUTION,
51 : AVERAGE_RESOLUTION,
52 : USER_RESOLUTION
53 : } ResolutionStrategy;
54 :
55 : typedef struct
56 : {
57 : int isFileOK;
58 : int nRasterXSize;
59 : int nRasterYSize;
60 : double adfGeoTransform[6];
61 : int nBlockXSize;
62 : int nBlockYSize;
63 : GDALDataType firstBandType;
64 : int* panHasNoData;
65 : double* padfNoDataValues;
66 : } DatasetProperty;
67 :
68 : typedef struct
69 : {
70 : GDALColorInterp colorInterpretation;
71 : GDALDataType dataType;
72 : GDALColorTableH colorTable;
73 : int bHasNoData;
74 : double noDataValue;
75 : } BandProperty;
76 :
77 : /************************************************************************/
78 : /* Usage() */
79 : /************************************************************************/
80 :
81 0 : static void Usage()
82 :
83 : {
84 : fprintf(stdout, "%s",
85 : "Usage: gdalbuildvrt [-tileindex field_name] [-resolution {highest|lowest|average|user}]\n"
86 : " [-tr xres yres] [-separate] [-allow_projection_difference] [-q]\n"
87 : " [-te xmin ymin xmax ymax] [-addalpha] \n"
88 : " [-srcnodata \"value [value...]\"] [-vrtnodata \"value [value...]\"] \n"
89 : " [-input_file_list my_liste.txt] output.vrt [gdalfile]*\n"
90 : "\n"
91 : "eg.\n"
92 : " % gdalbuildvrt doq_index.vrt doq/*.tif\n"
93 : " % gdalbuildvrt -input_file_list my_liste.txt doq_index.vrt\n"
94 : "\n"
95 : "NOTES:\n"
96 : " o With -separate, each files goes into a separate band in the VRT band. Otherwise,\n"
97 : " the files are considered as tiles of a larger mosaic.\n"
98 : " o The default tile index field is 'location' unless otherwise specified by -tileindex.\n"
99 : " o In case the resolution of all input files is not the same, the -resolution flag.\n"
100 : " enable the user to control the way the output resolution is computed. average is the default.\n"
101 : " o Input files may be any valid GDAL dataset or a GDAL raster tile index.\n"
102 : " o For a GDAL raster tile index, all entries will be added to the VRT.\n"
103 : " o If one GDAL dataset is made of several subdatasets and has 0 raster bands, its\n"
104 : " datasets will be added to the VRT rather than the dataset itself.\n"
105 : " o By default, only datasets of same projection and band characteristics may be added to the VRT.\n"
106 0 : );
107 0 : exit( 1 );
108 : }
109 :
110 :
111 38 : int GetSrcDstWin(DatasetProperty* psDP,
112 : double we_res, double ns_res,
113 : double minX, double minY, double maxX, double maxY,
114 : int* pnSrcXOff, int* pnSrcYOff, int* pnSrcXSize, int* pnSrcYSize,
115 : int* pnDstXOff, int* pnDstYOff, int* pnDstXSize, int* pnDstYSize)
116 : {
117 : /* Check that the destination bounding box intersects the source bounding box */
118 76 : if ( psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_X] +
119 : psDP->nRasterXSize *
120 38 : psDP->adfGeoTransform[GEOTRSFRM_WE_RES] < minX )
121 0 : return FALSE;
122 38 : if ( psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_X] > maxX )
123 0 : return FALSE;
124 76 : if ( psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_Y] +
125 : psDP->nRasterYSize *
126 38 : psDP->adfGeoTransform[GEOTRSFRM_NS_RES] > maxY )
127 0 : return FALSE;
128 38 : if ( psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_Y] < minY )
129 0 : return FALSE;
130 :
131 38 : *pnSrcXSize = psDP->nRasterXSize;
132 38 : *pnSrcYSize = psDP->nRasterYSize;
133 38 : if ( psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_X] < minX )
134 : {
135 1 : *pnSrcXOff = (int)((minX - psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_X]) /
136 1 : psDP->adfGeoTransform[GEOTRSFRM_WE_RES] + 0.5);
137 1 : *pnDstXOff = 0;
138 : }
139 : else
140 : {
141 37 : *pnSrcXOff = 0;
142 : *pnDstXOff = (int)
143 37 : (0.5 + (psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_X] - minX) / we_res);
144 : }
145 38 : if ( maxY < psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_Y])
146 : {
147 1 : *pnSrcYOff = (int)((psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_Y] - maxY) /
148 1 : -psDP->adfGeoTransform[GEOTRSFRM_NS_RES] + 0.5);
149 1 : *pnDstYOff = 0;
150 : }
151 : else
152 : {
153 37 : *pnSrcYOff = 0;
154 : *pnDstYOff = (int)
155 37 : (0.5 + (maxY - psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_Y]) / -ns_res);
156 : }
157 : *pnDstXSize = (int)
158 : (0.5 + psDP->nRasterXSize *
159 38 : psDP->adfGeoTransform[GEOTRSFRM_WE_RES] / we_res);
160 : *pnDstYSize = (int)
161 : (0.5 + psDP->nRasterYSize *
162 38 : psDP->adfGeoTransform[GEOTRSFRM_NS_RES] / ns_res);
163 :
164 38 : return TRUE;
165 : }
166 :
167 : /************************************************************************/
168 : /* GDALBuildVRT() */
169 : /************************************************************************/
170 :
171 12 : CPLErr GDALBuildVRT( const char* pszOutputFilename,
172 : int* pnInputFiles, char*** pppszInputFilenames,
173 : ResolutionStrategy resolutionStrategy,
174 : double we_res, double ns_res,
175 : double minX, double minY, double maxX, double maxY,
176 : int bSeparate, int bAllowProjectionDifference,
177 : int bAddAlpha,
178 : const char* pszSrcNoData, const char* pszVRTNoData,
179 : GDALProgressFunc pfnProgress, void * pProgressData)
180 : {
181 12 : char* projectionRef = NULL;
182 12 : int nBands = 0;
183 12 : BandProperty* bandProperties = NULL;
184 : int i,j;
185 : int rasterXSize;
186 : int rasterYSize;
187 12 : int nCount = 0;
188 12 : int bFirst = TRUE;
189 12 : VRTDatasetH hVRTDS = NULL;
190 12 : CPLErr eErr = CE_None;
191 :
192 12 : if( pfnProgress == NULL )
193 0 : pfnProgress = GDALDummyProgress;
194 :
195 12 : int bUserExtent = (minX != 0 || minY != 0 || maxX != 0 || maxY != 0);
196 12 : if (bUserExtent)
197 : {
198 2 : if (minX >= maxX || minY >= maxY )
199 : {
200 0 : CPLError(CE_Failure, CPLE_IllegalArg, "Invalid user extent");
201 0 : return CE_Failure;
202 : }
203 : }
204 :
205 12 : if (resolutionStrategy == USER_RESOLUTION)
206 : {
207 2 : if (we_res <= 0 || ns_res <= 0)
208 : {
209 0 : CPLError(CE_Failure, CPLE_IllegalArg, "Invalid user resolution");
210 0 : return CE_Failure;
211 : }
212 :
213 : /* We work with negative north-south resolution in all the following code */
214 2 : ns_res = -ns_res;
215 : }
216 : else
217 : {
218 10 : we_res = ns_res = 0;
219 : }
220 :
221 12 : char** ppszInputFilenames = *pppszInputFilenames;
222 12 : int nInputFiles = *pnInputFiles;
223 :
224 : DatasetProperty* psDatasetProperties =
225 12 : (DatasetProperty*) CPLCalloc(nInputFiles, sizeof(DatasetProperty));
226 :
227 :
228 12 : int bAllowSrcNoData = TRUE;
229 12 : double* padfSrcNoData = NULL;
230 12 : int nSrcNoDataCount = 0;
231 12 : if (pszSrcNoData != NULL)
232 : {
233 1 : if (EQUAL(pszSrcNoData, "none"))
234 : {
235 0 : bAllowSrcNoData = FALSE;
236 : }
237 : else
238 : {
239 1 : char **papszTokens = CSLTokenizeString( pszSrcNoData );
240 1 : nSrcNoDataCount = CSLCount(papszTokens);
241 1 : padfSrcNoData = (double *) CPLMalloc(sizeof(double) * nSrcNoDataCount);
242 2 : for(i=0;i<nSrcNoDataCount;i++)
243 1 : padfSrcNoData[i] = atof(papszTokens[i]);
244 1 : CSLDestroy(papszTokens);
245 : }
246 : }
247 :
248 :
249 12 : int bAllowVRTNoData = TRUE;
250 12 : double* padfVRTNoData = NULL;
251 12 : int nVRTNoDataCount = 0;
252 12 : if (pszVRTNoData != NULL)
253 : {
254 1 : if (EQUAL(pszVRTNoData, "none"))
255 : {
256 0 : bAllowVRTNoData = FALSE;
257 : }
258 : else
259 : {
260 1 : char **papszTokens = CSLTokenizeString( pszVRTNoData );
261 1 : nVRTNoDataCount = CSLCount(papszTokens);
262 1 : padfVRTNoData = (double *) CPLMalloc(sizeof(double) * nVRTNoDataCount);
263 2 : for(i=0;i<nVRTNoDataCount;i++)
264 1 : padfVRTNoData[i] = atof(papszTokens[i]);
265 1 : CSLDestroy(papszTokens);
266 : }
267 : }
268 :
269 52 : for(i=0;i<nInputFiles;i++)
270 : {
271 40 : const char* dsFileName = ppszInputFilenames[i];
272 :
273 40 : if (!pfnProgress( 1.0 * (i+1) / nInputFiles, NULL, pProgressData))
274 : {
275 0 : eErr = CE_Failure;
276 0 : goto end;
277 : }
278 :
279 40 : GDALDatasetH hDS = GDALOpen(ppszInputFilenames[i], GA_ReadOnly );
280 40 : psDatasetProperties[i].isFileOK = FALSE;
281 :
282 40 : if (hDS)
283 : {
284 40 : char** papszMetadata = GDALGetMetadata( hDS, "SUBDATASETS" );
285 40 : if( CSLCount(papszMetadata) > 0 && GDALGetRasterCount(hDS) == 0 )
286 : {
287 : psDatasetProperties =
288 0 : (DatasetProperty*) CPLMalloc((nInputFiles+CSLCount(papszMetadata))*sizeof(DatasetProperty));
289 :
290 : ppszInputFilenames = (char**)CPLRealloc(ppszInputFilenames,
291 0 : sizeof(char*) * (nInputFiles+CSLCount(papszMetadata)));
292 0 : int count = 1;
293 : char subdatasetNameKey[256];
294 0 : sprintf(subdatasetNameKey, "SUBDATASET_%d_NAME", count);
295 0 : while(*papszMetadata != NULL)
296 : {
297 0 : if (EQUALN(*papszMetadata, subdatasetNameKey, strlen(subdatasetNameKey)))
298 : {
299 0 : ppszInputFilenames[nInputFiles++] =
300 0 : CPLStrdup(*papszMetadata+strlen(subdatasetNameKey)+1);
301 0 : count++;
302 0 : sprintf(subdatasetNameKey, "SUBDATASET_%d_NAME", count);
303 : }
304 0 : papszMetadata++;
305 : }
306 0 : GDALClose(hDS);
307 0 : continue;
308 : }
309 :
310 40 : const char* proj = GDALGetProjectionRef(hDS);
311 40 : GDALGetGeoTransform(hDS, psDatasetProperties[i].adfGeoTransform);
312 80 : if (psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_ROTATION_PARAM1] != 0 ||
313 40 : psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_ROTATION_PARAM2] != 0)
314 : {
315 : CPLError(CE_Warning, CPLE_NotSupported,
316 : "gdalbuildvrt does not support rotated geo transforms. Skipping %s",
317 0 : dsFileName);
318 0 : GDALClose(hDS);
319 0 : continue;
320 : }
321 40 : if (psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_NS_RES] >= 0)
322 : {
323 : CPLError(CE_Warning, CPLE_NotSupported,
324 : "gdalbuildvrt does not support positive NS resolution. Skipping %s",
325 0 : dsFileName);
326 0 : GDALClose(hDS);
327 0 : continue;
328 : }
329 40 : psDatasetProperties[i].nRasterXSize = GDALGetRasterXSize(hDS);
330 40 : psDatasetProperties[i].nRasterYSize = GDALGetRasterYSize(hDS);
331 40 : double product_minX = psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_TOPLEFT_X];
332 40 : double product_maxY = psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_TOPLEFT_Y];
333 : double product_maxX = product_minX +
334 : GDALGetRasterXSize(hDS) *
335 40 : psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_WE_RES];
336 : double product_minY = product_maxY +
337 : GDALGetRasterYSize(hDS) *
338 40 : psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_NS_RES];
339 :
340 : GDALGetBlockSize(GDALGetRasterBand( hDS, 1 ),
341 40 : &psDatasetProperties[i].nBlockXSize,
342 80 : &psDatasetProperties[i].nBlockYSize);
343 :
344 40 : int _nBands = GDALGetRasterCount(hDS);
345 40 : if (_nBands == 0)
346 : {
347 : CPLError(CE_Warning, CPLE_AppDefined,
348 0 : "Skipping %s as it has no bands", dsFileName);
349 0 : GDALClose(hDS);
350 0 : continue;
351 : }
352 40 : else if (_nBands > 1 && bSeparate)
353 : {
354 : CPLError(CE_Warning, CPLE_AppDefined, "%s has %d bands. Only the first one will "
355 : "be taken into account in the -separate case",
356 0 : dsFileName, _nBands);
357 0 : _nBands = 1;
358 : }
359 :
360 : /* For the -separate case */
361 40 : psDatasetProperties[i].firstBandType = GDALGetRasterDataType(GDALGetRasterBand(hDS, 1));
362 :
363 40 : psDatasetProperties[i].padfNoDataValues = (double*)CPLMalloc(sizeof(double) * _nBands);
364 40 : psDatasetProperties[i].panHasNoData = (int*)CPLMalloc(sizeof(int) * _nBands);
365 85 : for(j=0;j<_nBands;j++)
366 : {
367 45 : if (nSrcNoDataCount > 0)
368 : {
369 2 : psDatasetProperties[i].panHasNoData[j] = TRUE;
370 2 : if (j < nSrcNoDataCount)
371 2 : psDatasetProperties[i].padfNoDataValues[j] = padfSrcNoData[j];
372 : else
373 0 : psDatasetProperties[i].padfNoDataValues[j] = padfSrcNoData[nSrcNoDataCount - 1];
374 : }
375 : else
376 : {
377 43 : psDatasetProperties[i].padfNoDataValues[j] =
378 : GDALGetRasterNoDataValue(GDALGetRasterBand(hDS, j+1),
379 43 : &psDatasetProperties[i].panHasNoData[j]);
380 : }
381 : }
382 :
383 40 : if (bFirst)
384 : {
385 12 : if (proj)
386 12 : projectionRef = CPLStrdup(proj);
387 12 : if (!bUserExtent)
388 : {
389 10 : minX = product_minX;
390 10 : minY = product_minY;
391 10 : maxX = product_maxX;
392 10 : maxY = product_maxY;
393 : }
394 12 : nBands = _nBands;
395 :
396 12 : if (!bSeparate)
397 : {
398 11 : bandProperties = (BandProperty*)CPLMalloc(nBands*sizeof(BandProperty));
399 24 : for(j=0;j<nBands;j++)
400 : {
401 13 : GDALRasterBandH hRasterBand = GDALGetRasterBand( hDS, j+1 );
402 13 : bandProperties[j].colorInterpretation =
403 13 : GDALGetRasterColorInterpretation(hRasterBand);
404 13 : bandProperties[j].dataType = GDALGetRasterDataType(hRasterBand);
405 13 : if (bandProperties[j].colorInterpretation == GCI_PaletteIndex)
406 : {
407 1 : bandProperties[j].colorTable =
408 1 : GDALGetRasterColorTable( hRasterBand );
409 1 : if (bandProperties[j].colorTable)
410 : {
411 1 : bandProperties[j].colorTable =
412 1 : GDALCloneColorTable(bandProperties[j].colorTable);
413 : }
414 : }
415 : else
416 12 : bandProperties[j].colorTable = 0;
417 :
418 13 : if (nVRTNoDataCount > 0)
419 : {
420 1 : bandProperties[j].bHasNoData = TRUE;
421 1 : if (j < nVRTNoDataCount)
422 1 : bandProperties[j].noDataValue = padfVRTNoData[j];
423 : else
424 0 : bandProperties[j].noDataValue = padfVRTNoData[nVRTNoDataCount - 1];
425 : }
426 : else
427 : {
428 12 : bandProperties[j].noDataValue =
429 12 : GDALGetRasterNoDataValue(hRasterBand, &bandProperties[j].bHasNoData);
430 : }
431 : }
432 : }
433 : }
434 : else
435 : {
436 28 : if ((proj != NULL && projectionRef == NULL) ||
437 : (proj == NULL && projectionRef != NULL) ||
438 : (proj != NULL && projectionRef != NULL && EQUAL(proj, projectionRef) == FALSE))
439 : {
440 1 : if (!bAllowProjectionDifference)
441 : {
442 : CPLError(CE_Warning, CPLE_NotSupported,
443 : "gdalbuildvrt does not support heterogenous projection. Skipping %s",
444 1 : dsFileName);
445 1 : GDALClose(hDS);
446 1 : continue;
447 : }
448 : }
449 27 : if (!bSeparate)
450 : {
451 24 : if (nBands != _nBands)
452 : {
453 : CPLError(CE_Warning, CPLE_NotSupported,
454 : "gdalbuildvrt does not support heterogenous band numbers. Skipping %s",
455 1 : dsFileName);
456 1 : GDALClose(hDS);
457 1 : continue;
458 : }
459 48 : for(j=0;j<nBands;j++)
460 : {
461 25 : GDALRasterBandH hRasterBand = GDALGetRasterBand( hDS, j+1 );
462 50 : if (bandProperties[j].colorInterpretation != GDALGetRasterColorInterpretation(hRasterBand) ||
463 25 : bandProperties[j].dataType != GDALGetRasterDataType(hRasterBand))
464 : {
465 : CPLError(CE_Warning, CPLE_NotSupported,
466 : "gdalbuildvrt does not support heterogenous band characteristics. Skipping %s",
467 0 : dsFileName);
468 0 : GDALClose(hDS);
469 : }
470 25 : if (bandProperties[j].colorTable)
471 : {
472 1 : GDALColorTableH colorTable = GDALGetRasterColorTable( hRasterBand );
473 2 : if (colorTable == NULL ||
474 1 : GDALGetColorEntryCount(colorTable) != GDALGetColorEntryCount(bandProperties[j].colorTable))
475 : {
476 : CPLError(CE_Warning, CPLE_NotSupported,
477 : "gdalbuildvrt does not support heterogenous band characteristics. Skipping %s",
478 0 : dsFileName);
479 0 : GDALClose(hDS);
480 0 : break;
481 : }
482 : /* We should check that the palette are the same too ! */
483 : }
484 : }
485 23 : if (j != nBands)
486 0 : continue;
487 : }
488 26 : if (!bUserExtent)
489 : {
490 23 : if (product_minX < minX) minX = product_minX;
491 23 : if (product_minY < minY) minY = product_minY;
492 23 : if (product_maxX > maxX) maxX = product_maxX;
493 23 : if (product_maxY > maxY) maxY = product_maxY;
494 : }
495 : }
496 :
497 38 : if (resolutionStrategy == AVERAGE_RESOLUTION)
498 : {
499 33 : we_res += psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_WE_RES];
500 33 : ns_res += psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_NS_RES];
501 : }
502 5 : else if (resolutionStrategy != USER_RESOLUTION)
503 : {
504 0 : if (bFirst)
505 : {
506 0 : we_res = psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_WE_RES];
507 0 : ns_res = psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_NS_RES];
508 : }
509 0 : else if (resolutionStrategy == HIGHEST_RESOLUTION)
510 : {
511 0 : we_res = MIN(we_res, psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_WE_RES]);
512 : /* Yes : as ns_res is negative, the highest resolution is the max value */
513 0 : ns_res = MAX(ns_res, psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_NS_RES]);
514 : }
515 : else
516 : {
517 0 : we_res = MAX(we_res, psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_WE_RES]);
518 : /* Yes : as ns_res is negative, the lowest resolution is the min value */
519 0 : ns_res = MIN(ns_res, psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_NS_RES]);
520 : }
521 : }
522 :
523 38 : psDatasetProperties[i].isFileOK = 1;
524 38 : nCount ++;
525 38 : bFirst = FALSE;
526 38 : GDALClose(hDS);
527 : }
528 : else
529 : {
530 : CPLError(CE_Warning, CPLE_AppDefined,
531 0 : "Warning : can't open %s. Skipping it", dsFileName);
532 : }
533 : }
534 :
535 12 : *pppszInputFilenames = ppszInputFilenames;
536 12 : *pnInputFiles = nInputFiles;
537 :
538 12 : if (nCount == 0)
539 0 : goto end;
540 :
541 12 : if (resolutionStrategy == AVERAGE_RESOLUTION)
542 : {
543 10 : we_res /= nCount;
544 10 : ns_res /= nCount;
545 : }
546 :
547 12 : rasterXSize = (int)(0.5 + (maxX - minX) / we_res);
548 12 : rasterYSize = (int)(0.5 + (maxY - minY) / -ns_res);
549 :
550 12 : if (rasterXSize == 0 || rasterYSize == 0)
551 : {
552 : CPLError(CE_Failure, CPLE_AppDefined,
553 0 : "Computed VRT dimension is invalid. You've probably specified unappropriate resolution.");
554 0 : goto end;
555 : }
556 :
557 12 : hVRTDS = VRTCreate(rasterXSize, rasterYSize);
558 12 : GDALSetDescription(hVRTDS, pszOutputFilename);
559 :
560 12 : if (projectionRef)
561 : {
562 12 : GDALSetProjection(hVRTDS, projectionRef);
563 : }
564 :
565 : double adfGeoTransform[6];
566 12 : adfGeoTransform[GEOTRSFRM_TOPLEFT_X] = minX;
567 12 : adfGeoTransform[GEOTRSFRM_WE_RES] = we_res;
568 12 : adfGeoTransform[GEOTRSFRM_ROTATION_PARAM1] = 0;
569 12 : adfGeoTransform[GEOTRSFRM_TOPLEFT_Y] = maxY;
570 12 : adfGeoTransform[GEOTRSFRM_ROTATION_PARAM2] = 0;
571 12 : adfGeoTransform[GEOTRSFRM_NS_RES] = ns_res;
572 12 : GDALSetGeoTransform(hVRTDS, adfGeoTransform);
573 :
574 12 : if (bSeparate)
575 : {
576 1 : int iBand = 1;
577 5 : for(i=0;i<nInputFiles;i++)
578 : {
579 4 : if (psDatasetProperties[i].isFileOK == 0)
580 0 : continue;
581 :
582 : int nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
583 : nDstXOff, nDstYOff, nDstXSize, nDstYSize;
584 4 : if ( ! GetSrcDstWin(&psDatasetProperties[i],
585 : we_res, ns_res, minX, minY, maxX, maxY,
586 : &nSrcXOff, &nSrcYOff, &nSrcXSize, &nSrcYSize,
587 : &nDstXOff, &nDstYOff, &nDstXSize, &nDstYSize) )
588 0 : continue;
589 :
590 4 : const char* dsFileName = ppszInputFilenames[i];
591 :
592 4 : GDALAddBand(hVRTDS, psDatasetProperties[i].firstBandType, NULL);
593 :
594 : GDALProxyPoolDatasetH hProxyDS =
595 : GDALProxyPoolDatasetCreate(dsFileName,
596 4 : psDatasetProperties[i].nRasterXSize,
597 4 : psDatasetProperties[i].nRasterYSize,
598 : GA_ReadOnly, TRUE, projectionRef,
599 12 : psDatasetProperties[i].adfGeoTransform);
600 : GDALProxyPoolDatasetAddSrcBandDescription(hProxyDS,
601 4 : psDatasetProperties[i].firstBandType,
602 4 : psDatasetProperties[i].nBlockXSize,
603 12 : psDatasetProperties[i].nBlockYSize);
604 :
605 : VRTSourcedRasterBandH hVRTBand =
606 4 : (VRTSourcedRasterBandH)GDALGetRasterBand(hVRTDS, iBand);
607 4 : if (bAllowSrcNoData && psDatasetProperties[i].panHasNoData[0])
608 : {
609 0 : GDALSetRasterNoDataValue(hVRTBand, psDatasetProperties[i].padfNoDataValues[0]);
610 : VRTAddComplexSource(hVRTBand, GDALGetRasterBand((GDALDatasetH)hProxyDS, 1),
611 : nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
612 : nDstXOff, nDstYOff, nDstXSize, nDstYSize,
613 0 : 0, 1, psDatasetProperties[i].padfNoDataValues[0]);
614 : }
615 : else
616 : /* Place the raster band at the right position in the VRT */
617 : VRTAddSimpleSource(hVRTBand, GDALGetRasterBand((GDALDatasetH)hProxyDS, 1),
618 : nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
619 : nDstXOff, nDstYOff, nDstXSize, nDstYSize,
620 4 : "near", VRT_NODATA_UNSET);
621 :
622 4 : GDALDereferenceDataset(hProxyDS);
623 :
624 4 : iBand ++;
625 : }
626 : }
627 : else
628 : {
629 24 : for(j=0;j<nBands;j++)
630 : {
631 : GDALRasterBandH hBand;
632 13 : GDALAddBand(hVRTDS, bandProperties[j].dataType, NULL);
633 13 : hBand = GDALGetRasterBand(hVRTDS, j+1);
634 13 : GDALSetRasterColorInterpretation(hBand, bandProperties[j].colorInterpretation);
635 13 : if (bandProperties[j].colorInterpretation == GCI_PaletteIndex)
636 : {
637 1 : GDALSetRasterColorTable(hBand, bandProperties[j].colorTable);
638 : }
639 13 : if (bAllowVRTNoData && bandProperties[j].bHasNoData)
640 4 : GDALSetRasterNoDataValue(hBand, bandProperties[j].noDataValue);
641 : }
642 :
643 11 : if (bAddAlpha)
644 : {
645 : GDALRasterBandH hBand;
646 0 : GDALAddBand(hVRTDS, GDT_Byte, NULL);
647 0 : hBand = GDALGetRasterBand(hVRTDS, nBands + 1);
648 0 : GDALSetRasterColorInterpretation(hBand, GCI_AlphaBand);
649 : }
650 :
651 47 : for(i=0;i<nInputFiles;i++)
652 : {
653 36 : if (psDatasetProperties[i].isFileOK == 0)
654 2 : continue;
655 :
656 : int nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
657 : nDstXOff, nDstYOff, nDstXSize, nDstYSize;
658 34 : if ( ! GetSrcDstWin(&psDatasetProperties[i],
659 : we_res, ns_res, minX, minY, maxX, maxY,
660 : &nSrcXOff, &nSrcYOff, &nSrcXSize, &nSrcYSize,
661 : &nDstXOff, &nDstYOff, &nDstXSize, &nDstYSize) )
662 0 : continue;
663 :
664 34 : const char* dsFileName = ppszInputFilenames[i];
665 :
666 : GDALProxyPoolDatasetH hProxyDS =
667 : GDALProxyPoolDatasetCreate(dsFileName,
668 34 : psDatasetProperties[i].nRasterXSize,
669 34 : psDatasetProperties[i].nRasterYSize,
670 : GA_ReadOnly, TRUE, projectionRef,
671 102 : psDatasetProperties[i].adfGeoTransform);
672 :
673 72 : for(j=0;j<nBands;j++)
674 : {
675 : GDALProxyPoolDatasetAddSrcBandDescription(hProxyDS,
676 38 : bandProperties[j].dataType,
677 38 : psDatasetProperties[i].nBlockXSize,
678 114 : psDatasetProperties[i].nBlockYSize);
679 : }
680 :
681 72 : for(j=0;j<nBands;j++)
682 : {
683 : VRTSourcedRasterBandH hVRTBand =
684 38 : (VRTSourcedRasterBandH)GDALGetRasterBand(hVRTDS, j + 1);
685 :
686 : /* Place the raster band at the right position in the VRT */
687 46 : if (bAllowSrcNoData && psDatasetProperties[i].panHasNoData[j])
688 : VRTAddComplexSource(hVRTBand, GDALGetRasterBand((GDALDatasetH)hProxyDS, j + 1),
689 : nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
690 : nDstXOff, nDstYOff, nDstXSize, nDstYSize,
691 8 : 0, 1, psDatasetProperties[i].padfNoDataValues[j]);
692 : else
693 : VRTAddSimpleSource(hVRTBand, GDALGetRasterBand((GDALDatasetH)hProxyDS, j + 1),
694 : nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
695 : nDstXOff, nDstYOff, nDstXSize, nDstYSize,
696 30 : "near", VRT_NODATA_UNSET);
697 : }
698 :
699 34 : if (bAddAlpha)
700 : {
701 : VRTSourcedRasterBandH hVRTBand =
702 0 : (VRTSourcedRasterBandH)GDALGetRasterBand(hVRTDS, nBands + 1);
703 : /* Little trick : we use an offset of 255 and a scaling of 0, so that in areas covered */
704 : /* by the source, the value of the alpha band will be 255, otherwise it will be 0 */
705 : VRTAddComplexSource(hVRTBand, GDALGetRasterBand((GDALDatasetH)hProxyDS, 1),
706 : nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
707 : nDstXOff, nDstYOff, nDstXSize, nDstYSize,
708 0 : 255, 0, VRT_NODATA_UNSET);
709 : }
710 :
711 34 : GDALDereferenceDataset(hProxyDS);
712 : }
713 : }
714 12 : GDALClose(hVRTDS);
715 :
716 : end:
717 52 : for(i=0;i<nInputFiles;i++)
718 : {
719 40 : CPLFree(psDatasetProperties[i].padfNoDataValues);
720 40 : CPLFree(psDatasetProperties[i].panHasNoData);
721 : }
722 12 : CPLFree(psDatasetProperties);
723 12 : if (!bSeparate && bandProperties != NULL)
724 : {
725 24 : for(j=0;j<nBands;j++)
726 : {
727 13 : GDALDestroyColorTable(bandProperties[j].colorTable);
728 : }
729 : }
730 12 : CPLFree(bandProperties);
731 12 : CPLFree(projectionRef);
732 12 : CPLFree(padfSrcNoData);
733 12 : CPLFree(padfVRTNoData);
734 :
735 12 : return eErr;
736 : }
737 :
738 : /************************************************************************/
739 : /* add_file_to_list() */
740 : /************************************************************************/
741 :
742 37 : static void add_file_to_list(const char* filename, const char* tile_index,
743 : int* pnInputFiles, char*** pppszInputFilenames)
744 : {
745 : #ifndef OGR_ENABLED
746 : CPLError(CE_Failure, CPLE_AppDefined, "OGR support needed to read tileindex");
747 : *pnInputFiles = 0;
748 : *pppszInputFilenames = NULL;
749 : #else
750 : OGRDataSourceH hDS;
751 : OGRLayerH hLayer;
752 : OGRFeatureDefnH hFDefn;
753 : int j, ti_field;
754 :
755 37 : int nInputFiles = *pnInputFiles;
756 37 : char** ppszInputFilenames = *pppszInputFilenames;
757 :
758 37 : if (EQUAL(CPLGetExtension(filename), "SHP"))
759 : {
760 1 : OGRRegisterAll();
761 :
762 : /* Handle GDALTIndex Shapefile as a special case */
763 1 : hDS = OGROpen( filename, FALSE, NULL );
764 1 : if( hDS == NULL )
765 : {
766 : fprintf( stderr, "Unable to open shapefile `%s'.\n",
767 0 : filename );
768 0 : exit(2);
769 : }
770 :
771 1 : hLayer = OGR_DS_GetLayer(hDS, 0);
772 :
773 1 : hFDefn = OGR_L_GetLayerDefn(hLayer);
774 :
775 1 : for( ti_field = 0; ti_field < OGR_FD_GetFieldCount(hFDefn); ti_field++ )
776 : {
777 1 : OGRFieldDefnH hFieldDefn = OGR_FD_GetFieldDefn( hFDefn, ti_field );
778 1 : const char* pszName = OGR_Fld_GetNameRef(hFieldDefn);
779 :
780 1 : if (strcmp(pszName, "LOCATION") == 0 && strcmp("LOCATION", tile_index) != 0 )
781 : {
782 : fprintf( stderr, "This shapefile seems to be a tile index of "
783 0 : "OGR features and not GDAL products.\n");
784 : }
785 1 : if( strcmp(pszName, tile_index) == 0 )
786 1 : break;
787 : }
788 :
789 1 : if( ti_field == OGR_FD_GetFieldCount(hFDefn) )
790 : {
791 : fprintf( stderr, "Unable to find field `%s' in DBF file `%s'.\n",
792 0 : tile_index, filename );
793 0 : return;
794 : }
795 :
796 : /* Load in memory existing file names in SHP */
797 1 : int nTileIndexFiles = OGR_L_GetFeatureCount(hLayer, TRUE);
798 1 : if (nTileIndexFiles == 0)
799 : {
800 0 : fprintf( stderr, "Tile index %s is empty. Skipping it.\n", filename);
801 0 : return;
802 : }
803 :
804 : ppszInputFilenames = (char**)CPLRealloc(ppszInputFilenames,
805 1 : sizeof(char*) * (nInputFiles+nTileIndexFiles));
806 5 : for(j=0;j<nTileIndexFiles;j++)
807 : {
808 4 : OGRFeatureH hFeat = OGR_L_GetNextFeature(hLayer);
809 8 : ppszInputFilenames[nInputFiles++] =
810 4 : CPLStrdup(OGR_F_GetFieldAsString(hFeat, ti_field ));
811 4 : OGR_F_Destroy(hFeat);
812 : }
813 :
814 1 : OGR_DS_Destroy( hDS );
815 : }
816 : else
817 : {
818 : ppszInputFilenames = (char**)CPLRealloc(ppszInputFilenames,
819 36 : sizeof(char*) * (nInputFiles+1));
820 36 : ppszInputFilenames[nInputFiles++] = CPLStrdup(filename);
821 : }
822 :
823 37 : *pnInputFiles = nInputFiles;
824 37 : *pppszInputFilenames = ppszInputFilenames;
825 : #endif
826 : }
827 :
828 : /************************************************************************/
829 : /* main() */
830 : /************************************************************************/
831 :
832 13 : int main( int nArgc, char ** papszArgv )
833 :
834 : {
835 13 : const char *tile_index = "location";
836 13 : const char *resolution = NULL;
837 13 : int nInputFiles = 0;
838 13 : char ** ppszInputFilenames = NULL;
839 13 : const char * pszOutputFilename = NULL;
840 : int i, iArg;
841 13 : int bSeparate = FALSE;
842 13 : int bAllowProjectionDifference = FALSE;
843 13 : int bQuiet = FALSE;
844 13 : GDALProgressFunc pfnProgress = NULL;
845 13 : double we_res = 0, ns_res = 0;
846 13 : double xmin = 0, ymin = 0, xmax = 0, ymax = 0;
847 13 : int bAddAlpha = FALSE;
848 13 : int bForceOverwrite = FALSE;
849 13 : const char* pszSrcNoData = NULL;
850 13 : const char* pszVRTNoData = NULL;
851 :
852 13 : GDALAllRegister();
853 :
854 13 : nArgc = GDALGeneralCmdLineProcessor( nArgc, &papszArgv, 0 );
855 13 : if( nArgc < 1 )
856 0 : exit( -nArgc );
857 :
858 : /* -------------------------------------------------------------------- */
859 : /* Parse commandline. */
860 : /* -------------------------------------------------------------------- */
861 65 : for( iArg = 1; iArg < nArgc; iArg++ )
862 : {
863 53 : if( EQUAL(papszArgv[iArg], "--utility_version") )
864 : {
865 : printf("%s was compiled against GDAL %s and is running against GDAL %s\n",
866 1 : papszArgv[0], GDAL_RELEASE_NAME, GDALVersionInfo("RELEASE_NAME"));
867 1 : return 0;
868 : }
869 52 : else if( EQUAL(papszArgv[iArg],"-tileindex") &&
870 : iArg + 1 < nArgc)
871 : {
872 0 : tile_index = papszArgv[++iArg];
873 : }
874 52 : else if( EQUAL(papszArgv[iArg],"-resolution") &&
875 : iArg + 1 < nArgc)
876 : {
877 0 : resolution = papszArgv[++iArg];
878 : }
879 53 : else if( EQUAL(papszArgv[iArg],"-input_file_list") &&
880 : iArg + 1 < nArgc)
881 : {
882 1 : const char* input_file_list = papszArgv[++iArg];
883 1 : FILE* f = VSIFOpen(input_file_list, "r");
884 1 : if (f)
885 : {
886 4 : while(1)
887 : {
888 5 : const char* filename = CPLReadLine(f);
889 5 : if (filename == NULL)
890 : break;
891 : add_file_to_list(filename, tile_index,
892 4 : &nInputFiles, &ppszInputFilenames);
893 : }
894 1 : VSIFClose(f);
895 : }
896 : }
897 51 : else if ( EQUAL(papszArgv[iArg],"-separate") )
898 : {
899 1 : bSeparate = TRUE;
900 : }
901 50 : else if ( EQUAL(papszArgv[iArg],"-allow_projection_difference") )
902 : {
903 0 : bAllowProjectionDifference = TRUE;
904 : }
905 : /* Alternate syntax for output file */
906 50 : else if( EQUAL(papszArgv[iArg],"-o") &&
907 : iArg + 1 < nArgc)
908 : {
909 0 : pszOutputFilename = papszArgv[++iArg];
910 : }
911 50 : else if ( EQUAL(papszArgv[iArg],"-q") || EQUAL(papszArgv[iArg],"-quiet") )
912 : {
913 0 : bQuiet = TRUE;
914 : }
915 52 : else if ( EQUAL(papszArgv[iArg],"-tr") && iArg + 2 < nArgc)
916 : {
917 2 : we_res = atof(papszArgv[++iArg]);
918 2 : ns_res = atof(papszArgv[++iArg]);
919 : }
920 50 : else if ( EQUAL(papszArgv[iArg],"-te") && iArg + 4 < nArgc)
921 : {
922 2 : xmin = atof(papszArgv[++iArg]);
923 2 : ymin = atof(papszArgv[++iArg]);
924 2 : xmax = atof(papszArgv[++iArg]);
925 2 : ymax = atof(papszArgv[++iArg]);
926 : }
927 46 : else if ( EQUAL(papszArgv[iArg],"-addalpha") )
928 : {
929 0 : bAddAlpha = TRUE;
930 : }
931 46 : else if ( EQUAL(papszArgv[iArg],"-overwrite") )
932 : {
933 0 : bForceOverwrite = TRUE;
934 : }
935 47 : else if ( EQUAL(papszArgv[iArg],"-srcnodata") && iArg + 1 < nArgc)
936 : {
937 1 : pszSrcNoData = papszArgv[++iArg];
938 : }
939 45 : else if ( EQUAL(papszArgv[iArg],"-vrtnodata") && iArg + 1 < nArgc)
940 : {
941 0 : pszVRTNoData = papszArgv[++iArg];
942 : }
943 45 : else if ( papszArgv[iArg][0] == '-' )
944 : {
945 0 : printf("Unrecognized option : %s\n", papszArgv[iArg]);
946 0 : exit(-1);
947 : }
948 45 : else if( pszOutputFilename == NULL )
949 : {
950 12 : pszOutputFilename = papszArgv[iArg];
951 : }
952 : else
953 : {
954 33 : add_file_to_list(papszArgv[iArg], tile_index,
955 66 : &nInputFiles, &ppszInputFilenames);
956 : }
957 : }
958 :
959 12 : if( pszOutputFilename == NULL || nInputFiles == 0 )
960 0 : Usage();
961 :
962 12 : if (!bQuiet)
963 12 : pfnProgress = GDALTermProgress;
964 :
965 : /* Avoid overwriting a non VRT dataset if the user did not put the */
966 : /* filenames in the right order */
967 : VSIStatBuf sBuf;
968 12 : if (!bForceOverwrite)
969 : {
970 12 : int bExists = (VSIStat(pszOutputFilename, &sBuf) == 0);
971 12 : if (bExists)
972 : {
973 6 : GDALDriverH hDriver = GDALIdentifyDriver( pszOutputFilename, NULL );
974 6 : if (hDriver && !EQUAL(GDALGetDriverShortName(hDriver), "VRT"))
975 : {
976 : fprintf(stderr,
977 : "'%s' is an existing GDAL dataset managed by %s driver.\n"
978 : "There is an high chance you did not put filenames in the right order.\n"
979 : "If you want to overwrite %s, add -overwrite option to the command line.\n\n",
980 0 : pszOutputFilename, GDALGetDriverShortName(hDriver), pszOutputFilename);
981 0 : Usage();
982 : }
983 : }
984 : }
985 :
986 12 : if (we_res != 0 && ns_res != 0 &&
987 : resolution != NULL && !EQUAL(resolution, "user"))
988 : {
989 0 : fprintf(stderr, "-tr option is not compatible with -resolution %s\n", resolution);
990 0 : Usage();
991 : }
992 :
993 12 : if (bAddAlpha && bSeparate)
994 : {
995 0 : fprintf(stderr, "-addalpha option is not compatible with -separate\n");
996 0 : Usage();
997 : }
998 :
999 12 : ResolutionStrategy eStrategy = AVERAGE_RESOLUTION;
1000 24 : if ( resolution == NULL || EQUAL(resolution, "user") )
1001 : {
1002 14 : if ( we_res != 0 || ns_res != 0)
1003 2 : eStrategy = USER_RESOLUTION;
1004 10 : else if ( resolution != NULL && EQUAL(resolution, "user") )
1005 : {
1006 0 : fprintf(stderr, "-tr option must be used with -resolution user\n");
1007 0 : Usage();
1008 : }
1009 : }
1010 0 : else if ( EQUAL(resolution, "average") )
1011 0 : eStrategy = AVERAGE_RESOLUTION;
1012 0 : else if ( EQUAL(resolution, "highest") )
1013 0 : eStrategy = HIGHEST_RESOLUTION;
1014 0 : else if ( EQUAL(resolution, "lowest") )
1015 0 : eStrategy = LOWEST_RESOLUTION;
1016 : else
1017 : {
1018 0 : fprintf(stderr, "invalid value (%s) for -resolution\n", resolution);
1019 0 : Usage();
1020 : }
1021 :
1022 : /* If -srcnodata is specified, use it as the -vrtnodata if the latter is not */
1023 : /* specified */
1024 12 : if (pszSrcNoData != NULL && pszVRTNoData == NULL)
1025 1 : pszVRTNoData = pszSrcNoData;
1026 :
1027 : GDALBuildVRT(pszOutputFilename, &nInputFiles, &ppszInputFilenames,
1028 : eStrategy, we_res, ns_res, xmin, ymin, xmax, ymax,
1029 : bSeparate, bAllowProjectionDifference, bAddAlpha,
1030 12 : pszSrcNoData, pszVRTNoData, pfnProgress, NULL);
1031 :
1032 52 : for(i=0;i<nInputFiles;i++)
1033 : {
1034 40 : CPLFree(ppszInputFilenames[i]);
1035 : }
1036 12 : CPLFree(ppszInputFilenames);
1037 :
1038 :
1039 12 : CSLDestroy( papszArgv );
1040 12 : GDALDumpOpenDatasets( stderr );
1041 12 : GDALDestroyDriverManager();
1042 :
1043 12 : return 0;
1044 : }
|
