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