1 : /******************************************************************************
2 : * $Id: gdaltindex.c 25582 2013-01-29 21:13:43Z rouault $
3 : *
4 : * Project: MapServer
5 : * Purpose: Commandline App to build tile index for raster files.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2001, Frank Warmerdam, DM Solutions Group Inc
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 "ogr_api.h"
31 : #include "ogr_srs_api.h"
32 : #include "gdal.h"
33 : #include "cpl_port.h"
34 : #include "cpl_conv.h"
35 : #include "cpl_string.h"
36 :
37 : CPL_CVSID("$Id: gdaltindex.c 25582 2013-01-29 21:13:43Z rouault $");
38 :
39 : /************************************************************************/
40 : /* Usage() */
41 : /************************************************************************/
42 :
43 0 : static void Usage(const char* pszErrorMsg)
44 :
45 : {
46 0 : fprintf(stdout, "%s",
47 : "\n"
48 : "Usage: gdaltindex [-tileindex field_name] [-write_absolute_path] \n"
49 : " [-skip_different_projection] [-t_srs target_srs]\n"
50 : " index_file [gdal_file]*\n"
51 : "\n"
52 : "eg.\n"
53 : " % gdaltindex doq_index.shp doq/*.tif\n"
54 : "\n"
55 : "NOTES:\n"
56 : " o The shapefile (index_file) will be created if it doesn't already exist.\n"
57 : " o The default tile index field is 'location'.\n"
58 : " o Raster filenames will be put in the file exactly as they are specified\n"
59 : " on the commandline unless the option -write_absolute_path is used.\n"
60 : " o If -skip_different_projection is specified, only files with same projection ref\n"
61 : " as files already inserted in the tileindex will be inserted (unless t_srs is specified).\n"
62 : " o If -t_srs is specified, geometries of input files will be transformed to the desired\n"
63 : " target coordinate reference system.\n"
64 : " Note that using this option generates files that are NOT compatible with MapServer.\n"
65 : " o Simple rectangular polygons are generated in the same coordinate reference system\n"
66 : " as the rasters, or in target reference system if the -t_srs option is used.\n");
67 :
68 0 : if( pszErrorMsg != NULL )
69 0 : fprintf(stderr, "\nFAILURE: %s\n", pszErrorMsg);
70 :
71 0 : exit(1);
72 : }
73 :
74 : /************************************************************************/
75 : /* main() */
76 : /************************************************************************/
77 :
78 : #define CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(nExtraArg) \
79 : do { if (i_arg + nExtraArg >= argc) \
80 : Usage(CPLSPrintf("%s option requires %d argument(s)", argv[i_arg], nExtraArg)); } while(0)
81 :
82 7 : int main(int argc, char *argv[])
83 : {
84 7 : const char *index_filename = NULL;
85 7 : const char *tile_index = "location";
86 : int i_arg, ti_field;
87 : OGRDataSourceH hTileIndexDS;
88 7 : OGRLayerH hLayer = NULL;
89 : OGRFeatureDefnH hFDefn;
90 7 : int write_absolute_path = FALSE;
91 7 : char* current_path = NULL;
92 : int i;
93 : int nExistingFiles;
94 7 : int skip_different_projection = FALSE;
95 7 : char** existingFilesTab = NULL;
96 7 : int alreadyExistingProjectionRefValid = FALSE;
97 7 : char* alreadyExistingProjectionRef = NULL;
98 : char* index_filename_mod;
99 : int bExists;
100 : VSIStatBuf sStatBuf;
101 7 : const char *pszTargetSRS = "";
102 7 : int bSetTargetSRS = FALSE;
103 7 : OGRSpatialReferenceH hTargetSRS = NULL;
104 :
105 : /* Check that we are running against at least GDAL 1.4 */
106 : /* Note to developers : if we use newer API, please change the requirement */
107 7 : if (atoi(GDALVersionInfo("VERSION_NUM")) < 1400)
108 : {
109 0 : fprintf(stderr, "At least, GDAL >= 1.4.0 is required for this version of %s, "
110 : "which was compiled against GDAL %s\n", argv[0], GDAL_RELEASE_NAME);
111 0 : exit(1);
112 : }
113 :
114 7 : GDALAllRegister();
115 7 : OGRRegisterAll();
116 :
117 7 : argc = GDALGeneralCmdLineProcessor( argc, &argv, 0 );
118 7 : if( argc < 1 )
119 0 : exit( -argc );
120 :
121 : /* -------------------------------------------------------------------- */
122 : /* Get commandline arguments other than the GDAL raster filenames. */
123 : /* -------------------------------------------------------------------- */
124 9 : for( i_arg = 1; i_arg < argc; i_arg++ )
125 : {
126 9 : if( EQUAL(argv[i_arg], "--utility_version") )
127 : {
128 1 : printf("%s was compiled against GDAL %s and is running against GDAL %s\n",
129 : argv[0], GDAL_RELEASE_NAME, GDALVersionInfo("RELEASE_NAME"));
130 1 : return 0;
131 : }
132 8 : else if( EQUAL(argv[i_arg],"--help") )
133 0 : Usage(NULL);
134 8 : else if( strcmp(argv[i_arg],"-tileindex") == 0 )
135 : {
136 0 : CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
137 0 : tile_index = argv[++i_arg];
138 : }
139 8 : else if( strcmp(argv[i_arg],"-t_srs") == 0 )
140 : {
141 1 : CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
142 1 : pszTargetSRS = argv[++i_arg];
143 1 : bSetTargetSRS = TRUE;
144 : }
145 7 : else if ( strcmp(argv[i_arg],"-write_absolute_path") == 0 )
146 : {
147 0 : write_absolute_path = TRUE;
148 : }
149 7 : else if ( strcmp(argv[i_arg],"-skip_different_projection") == 0 )
150 : {
151 1 : skip_different_projection = TRUE;
152 : }
153 6 : else if( argv[i_arg][0] == '-' )
154 0 : Usage(CPLSPrintf("Unkown option name '%s'", argv[i_arg]));
155 6 : else if( index_filename == NULL )
156 : {
157 6 : index_filename = argv[i_arg];
158 6 : i_arg++;
159 6 : break;
160 : }
161 : }
162 :
163 6 : if( index_filename == NULL )
164 0 : Usage("No index filename specified.");
165 6 : if( i_arg == argc )
166 0 : Usage("No file to index specified.");
167 :
168 : /* -------------------------------------------------------------------- */
169 : /* Create and validate target SRS if given. */
170 : /* -------------------------------------------------------------------- */
171 6 : if( bSetTargetSRS )
172 : {
173 1 : if ( skip_different_projection )
174 : {
175 0 : fprintf( stderr,
176 : "Warning : -skip_different_projection does not apply "
177 : "when -t_srs is requested.\n" );
178 : }
179 1 : hTargetSRS = OSRNewSpatialReference("");
180 1 : if( OSRSetFromUserInput( hTargetSRS, pszTargetSRS ) != CE_None )
181 : {
182 0 : OSRDestroySpatialReference( hTargetSRS );
183 0 : fprintf( stderr, "Invalid target SRS `%s'.\n",
184 : pszTargetSRS );
185 0 : exit(1);
186 : }
187 : }
188 :
189 : /* -------------------------------------------------------------------- */
190 : /* Open or create the target shapefile and DBF file. */
191 : /* -------------------------------------------------------------------- */
192 6 : index_filename_mod = CPLStrdup(CPLResetExtension(index_filename, "shp"));
193 :
194 6 : bExists = (VSIStat(index_filename_mod, &sStatBuf) == 0);
195 6 : if (!bExists)
196 : {
197 2 : CPLFree(index_filename_mod);
198 2 : index_filename_mod = CPLStrdup(CPLResetExtension(index_filename, "SHP"));
199 2 : bExists = (VSIStat(index_filename_mod, &sStatBuf) == 0);
200 : }
201 6 : CPLFree(index_filename_mod);
202 :
203 6 : if (bExists)
204 : {
205 4 : hTileIndexDS = OGROpen( index_filename, TRUE, NULL );
206 4 : if (hTileIndexDS != NULL)
207 : {
208 4 : hLayer = OGR_DS_GetLayer(hTileIndexDS, 0);
209 : }
210 : }
211 : else
212 : {
213 : OGRSFDriverH hDriver;
214 2 : const char* pszDriverName = "ESRI Shapefile";
215 :
216 2 : printf( "Creating new index file...\n" );
217 2 : hDriver = OGRGetDriverByName( pszDriverName );
218 2 : if( hDriver == NULL )
219 : {
220 0 : printf( "%s driver not available.\n", pszDriverName );
221 0 : exit( 1 );
222 : }
223 :
224 2 : hTileIndexDS = OGR_Dr_CreateDataSource( hDriver, index_filename, NULL );
225 2 : if (hTileIndexDS)
226 : {
227 2 : char* pszLayerName = CPLStrdup(CPLGetBasename(index_filename));
228 :
229 : /* get spatial reference for output file from target SRS (if set) */
230 : /* or from first input file */
231 2 : OGRSpatialReferenceH hSpatialRef = NULL;
232 2 : if( bSetTargetSRS )
233 : {
234 0 : hSpatialRef = OSRClone( hTargetSRS );
235 : }
236 : else
237 : {
238 2 : GDALDatasetH hDS = GDALOpen( argv[i_arg], GA_ReadOnly );
239 2 : if (hDS)
240 : {
241 2 : const char* pszWKT = GDALGetProjectionRef(hDS);
242 2 : if (pszWKT != NULL && pszWKT[0] != '\0')
243 : {
244 2 : hSpatialRef = OSRNewSpatialReference(pszWKT);
245 : }
246 2 : GDALClose(hDS);
247 : }
248 : }
249 :
250 2 : hLayer = OGR_DS_CreateLayer( hTileIndexDS, pszLayerName, hSpatialRef, wkbPolygon, NULL );
251 2 : CPLFree(pszLayerName);
252 2 : if (hSpatialRef)
253 2 : OSRRelease(hSpatialRef);
254 :
255 2 : if (hLayer)
256 : {
257 2 : OGRFieldDefnH hFieldDefn = OGR_Fld_Create( tile_index, OFTString );
258 2 : OGR_Fld_SetWidth( hFieldDefn, 255);
259 2 : OGR_L_CreateField( hLayer, hFieldDefn, TRUE );
260 2 : OGR_Fld_Destroy(hFieldDefn);
261 : }
262 : }
263 : }
264 :
265 6 : if( hTileIndexDS == NULL || hLayer == NULL )
266 : {
267 0 : fprintf( stderr, "Unable to open/create shapefile `%s'.\n",
268 : index_filename );
269 0 : exit(2);
270 : }
271 :
272 6 : hFDefn = OGR_L_GetLayerDefn(hLayer);
273 :
274 6 : for( ti_field = 0; ti_field < OGR_FD_GetFieldCount(hFDefn); ti_field++ )
275 : {
276 6 : OGRFieldDefnH hFieldDefn = OGR_FD_GetFieldDefn( hFDefn, ti_field );
277 6 : if( strcmp(OGR_Fld_GetNameRef(hFieldDefn), tile_index) == 0 )
278 6 : break;
279 : }
280 :
281 6 : if( ti_field == OGR_FD_GetFieldCount(hFDefn) )
282 : {
283 0 : fprintf( stderr, "Unable to find field `%s' in DBF file `%s'.\n",
284 : tile_index, index_filename );
285 0 : exit(2);
286 : }
287 :
288 : /* Load in memory existing file names in SHP */
289 6 : nExistingFiles = OGR_L_GetFeatureCount(hLayer, FALSE);
290 6 : if (nExistingFiles)
291 : {
292 : OGRFeatureH hFeature;
293 4 : existingFilesTab = (char**)CPLMalloc(nExistingFiles * sizeof(char*));
294 18 : for(i=0;i<nExistingFiles;i++)
295 : {
296 14 : hFeature = OGR_L_GetNextFeature(hLayer);
297 14 : existingFilesTab[i] = CPLStrdup(OGR_F_GetFieldAsString( hFeature, ti_field ));
298 14 : if (i == 0)
299 : {
300 4 : GDALDatasetH hDS = GDALOpen(existingFilesTab[i], GA_ReadOnly );
301 4 : if (hDS)
302 : {
303 4 : alreadyExistingProjectionRefValid = TRUE;
304 4 : alreadyExistingProjectionRef = CPLStrdup(GDALGetProjectionRef(hDS));
305 4 : GDALClose(hDS);
306 : }
307 : }
308 14 : OGR_F_Destroy( hFeature );
309 : }
310 : }
311 :
312 6 : if (write_absolute_path)
313 : {
314 0 : current_path = CPLGetCurrentDir();
315 0 : if (current_path == NULL)
316 : {
317 0 : fprintf( stderr, "This system does not support the CPLGetCurrentDir call. "
318 : "The option -write_absolute_path will have no effect\n");
319 0 : write_absolute_path = FALSE;
320 : }
321 : }
322 :
323 : /* -------------------------------------------------------------------- */
324 : /* loop over GDAL files, processing. */
325 : /* -------------------------------------------------------------------- */
326 20 : for( ; i_arg < argc; i_arg++ )
327 : {
328 : GDALDatasetH hDS;
329 : double adfGeoTransform[6];
330 : double adfX[5], adfY[5];
331 : int nXSize, nYSize;
332 : char* fileNameToWrite;
333 : const char* projectionRef;
334 : VSIStatBuf sStatBuf;
335 : int k;
336 : OGRFeatureH hFeature;
337 : OGRGeometryH hPoly, hRing;
338 :
339 : /* Make sure it is a file before building absolute path name */
340 14 : if (write_absolute_path && CPLIsFilenameRelative( argv[i_arg] ) &&
341 0 : VSIStat( argv[i_arg], &sStatBuf ) == 0)
342 : {
343 0 : fileNameToWrite = CPLStrdup(CPLProjectRelativeFilename(current_path, argv[i_arg]));
344 : }
345 : else
346 : {
347 14 : fileNameToWrite = CPLStrdup(argv[i_arg]);
348 : }
349 :
350 : /* Checks that file is not already in tileindex */
351 32 : for(i=0;i<nExistingFiles;i++)
352 : {
353 22 : if (EQUAL(fileNameToWrite, existingFilesTab[i]))
354 : {
355 4 : fprintf(stderr, "File %s is already in tileindex. Skipping it.\n",
356 : fileNameToWrite);
357 4 : break;
358 : }
359 : }
360 14 : if (i != nExistingFiles)
361 : {
362 4 : CPLFree(fileNameToWrite);
363 4 : continue;
364 : }
365 :
366 10 : hDS = GDALOpen( argv[i_arg], GA_ReadOnly );
367 10 : if( hDS == NULL )
368 : {
369 0 : fprintf( stderr, "Unable to open %s, skipping.\n",
370 0 : argv[i_arg] );
371 0 : CPLFree(fileNameToWrite);
372 0 : continue;
373 : }
374 :
375 10 : GDALGetGeoTransform( hDS, adfGeoTransform );
376 10 : if( adfGeoTransform[0] == 0.0
377 0 : && adfGeoTransform[1] == 1.0
378 0 : && adfGeoTransform[3] == 0.0
379 0 : && ABS(adfGeoTransform[5]) == 1.0 )
380 : {
381 0 : fprintf( stderr,
382 : "It appears no georeferencing is available for\n"
383 : "`%s', skipping.\n",
384 0 : argv[i_arg] );
385 0 : GDALClose( hDS );
386 0 : CPLFree(fileNameToWrite);
387 0 : continue;
388 : }
389 :
390 10 : projectionRef = GDALGetProjectionRef(hDS);
391 :
392 : /* if not set target srs, test that the current file uses same projection as others */
393 10 : if( !bSetTargetSRS )
394 : {
395 9 : if (alreadyExistingProjectionRefValid)
396 : {
397 : int projectionRefNotNull, alreadyExistingProjectionRefNotNull;
398 7 : projectionRefNotNull = projectionRef && projectionRef[0];
399 7 : alreadyExistingProjectionRefNotNull = alreadyExistingProjectionRef && alreadyExistingProjectionRef[0];
400 14 : if ((projectionRefNotNull &&
401 : alreadyExistingProjectionRefNotNull &&
402 7 : EQUAL(projectionRef, alreadyExistingProjectionRef) == 0) ||
403 : (projectionRefNotNull != alreadyExistingProjectionRefNotNull))
404 : {
405 2 : fprintf(stderr, "Warning : %s is not using the same projection system as "
406 : "other files in the tileindex.\n"
407 : "This may cause problems when using it in MapServer for example.\n"
408 : "Use -t_srs option to set target projection system (not supported by MapServer).\n"
409 1 : "%s\n", argv[i_arg],
410 : (skip_different_projection) ? "Skipping this file." : "");
411 1 : if (skip_different_projection)
412 : {
413 1 : CPLFree(fileNameToWrite);
414 1 : GDALClose( hDS );
415 1 : continue;
416 : }
417 : }
418 : }
419 : else
420 : {
421 2 : alreadyExistingProjectionRefValid = TRUE;
422 2 : alreadyExistingProjectionRef = CPLStrdup(projectionRef);
423 : }
424 : }
425 :
426 9 : nXSize = GDALGetRasterXSize( hDS );
427 9 : nYSize = GDALGetRasterYSize( hDS );
428 :
429 18 : adfX[0] = adfGeoTransform[0]
430 9 : + 0 * adfGeoTransform[1]
431 9 : + 0 * adfGeoTransform[2];
432 18 : adfY[0] = adfGeoTransform[3]
433 9 : + 0 * adfGeoTransform[4]
434 9 : + 0 * adfGeoTransform[5];
435 :
436 18 : adfX[1] = adfGeoTransform[0]
437 9 : + nXSize * adfGeoTransform[1]
438 18 : + 0 * adfGeoTransform[2];
439 18 : adfY[1] = adfGeoTransform[3]
440 9 : + nXSize * adfGeoTransform[4]
441 18 : + 0 * adfGeoTransform[5];
442 :
443 18 : adfX[2] = adfGeoTransform[0]
444 9 : + nXSize * adfGeoTransform[1]
445 18 : + nYSize * adfGeoTransform[2];
446 18 : adfY[2] = adfGeoTransform[3]
447 9 : + nXSize * adfGeoTransform[4]
448 18 : + nYSize * adfGeoTransform[5];
449 :
450 18 : adfX[3] = adfGeoTransform[0]
451 9 : + 0 * adfGeoTransform[1]
452 9 : + nYSize * adfGeoTransform[2];
453 18 : adfY[3] = adfGeoTransform[3]
454 9 : + 0 * adfGeoTransform[4]
455 9 : + nYSize * adfGeoTransform[5];
456 :
457 18 : adfX[4] = adfGeoTransform[0]
458 9 : + 0 * adfGeoTransform[1]
459 9 : + 0 * adfGeoTransform[2];
460 18 : adfY[4] = adfGeoTransform[3]
461 9 : + 0 * adfGeoTransform[4]
462 9 : + 0 * adfGeoTransform[5];
463 :
464 : /* if set target srs, do the forward transformation of all points */
465 9 : if( bSetTargetSRS )
466 : {
467 1 : OGRSpatialReferenceH hSourceSRS = NULL;
468 1 : OGRCoordinateTransformationH hCT = NULL;
469 1 : hSourceSRS = OSRNewSpatialReference( projectionRef );
470 1 : if( hSourceSRS && !OSRIsSame( hSourceSRS, hTargetSRS ) )
471 : {
472 1 : hCT = OCTNewCoordinateTransformation( hSourceSRS, hTargetSRS );
473 1 : if( hCT == NULL || !OCTTransform( hCT, 5, adfX, adfY, NULL ) )
474 : {
475 0 : fprintf( stderr,
476 : "Warning : unable to transform points from source SRS `%s' to target SRS `%s'\n"
477 : "for file `%s' - file skipped\n",
478 : projectionRef, pszTargetSRS, fileNameToWrite );
479 0 : if ( hCT )
480 0 : OCTDestroyCoordinateTransformation( hCT );
481 0 : if ( hSourceSRS )
482 0 : OSRDestroySpatialReference( hSourceSRS );
483 0 : continue;
484 : }
485 1 : if ( hCT )
486 1 : OCTDestroyCoordinateTransformation( hCT );
487 : }
488 1 : if ( hSourceSRS )
489 1 : OSRDestroySpatialReference( hSourceSRS );
490 : }
491 :
492 9 : hFeature = OGR_F_Create( OGR_L_GetLayerDefn( hLayer ) );
493 9 : OGR_F_SetFieldString( hFeature, ti_field, fileNameToWrite );
494 :
495 9 : hPoly = OGR_G_CreateGeometry(wkbPolygon);
496 9 : hRing = OGR_G_CreateGeometry(wkbLinearRing);
497 54 : for(k=0;k<5;k++)
498 45 : OGR_G_SetPoint_2D(hRing, k, adfX[k], adfY[k]);
499 9 : OGR_G_AddGeometryDirectly( hPoly, hRing );
500 9 : OGR_F_SetGeometryDirectly( hFeature, hPoly );
501 :
502 9 : if( OGR_L_CreateFeature( hLayer, hFeature ) != OGRERR_NONE )
503 : {
504 0 : printf( "Failed to create feature in shapefile.\n" );
505 0 : break;
506 : }
507 :
508 9 : OGR_F_Destroy( hFeature );
509 :
510 :
511 9 : CPLFree(fileNameToWrite);
512 :
513 9 : GDALClose( hDS );
514 : }
515 :
516 6 : CPLFree(current_path);
517 :
518 6 : if (nExistingFiles)
519 : {
520 18 : for(i=0;i<nExistingFiles;i++)
521 : {
522 14 : CPLFree(existingFilesTab[i]);
523 : }
524 4 : CPLFree(existingFilesTab);
525 : }
526 6 : CPLFree(alreadyExistingProjectionRef);
527 :
528 6 : if ( hTargetSRS )
529 1 : OSRDestroySpatialReference( hTargetSRS );
530 :
531 6 : OGR_DS_Destroy( hTileIndexDS );
532 :
533 6 : GDALDestroyDriverManager();
534 6 : OGRCleanupAll();
535 6 : CSLDestroy(argv);
536 :
537 6 : exit( 0 );
538 : }
|
