1 : /******************************************************************************
2 : * $Id: gdaltindex.c 19806 2010-06-06 11:56:19Z 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 19806 2010-06-06 11:56:19Z rouault $");
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] index_file [gdal_file]*\n"
50 : "\n"
51 : "eg.\n"
52 : " % gdaltindex doq_index.shp doq/*.tif\n"
53 : "\n"
54 : "NOTES:\n"
55 : " o The shapefile (index_file) will be created if it doesn't already exist.\n"
56 : " o The default tile index field is 'location'.\n"
57 : " o Raster filenames will be put in the file exactly as they are specified\n"
58 : " on the commandline unless the option -write_absolute_path is used.\n"
59 : " o If -skip_different_projection is specified, only files with same projection ref\n"
60 : " as files already inserted in the tileindex will be inserted.\n"
61 : " o Simple rectangular polygons are generated in the same\n"
62 : " coordinate system as the rasters.\n" );
63 0 : exit(1);
64 : }
65 :
66 : /************************************************************************/
67 : /* main() */
68 : /************************************************************************/
69 :
70 12 : int main(int argc, char *argv[])
71 : {
72 12 : const char *index_filename = NULL;
73 12 : const char *tile_index = "location";
74 : int i_arg, ti_field;
75 : OGRDataSourceH hTileIndexDS;
76 12 : OGRLayerH hLayer = NULL;
77 : OGRFeatureDefnH hFDefn;
78 12 : int write_absolute_path = FALSE;
79 12 : char* current_path = NULL;
80 : int i;
81 : int nExistingFiles;
82 12 : int skip_different_projection = FALSE;
83 12 : char** existingFilesTab = NULL;
84 12 : int alreadyExistingProjectionRefValid = FALSE;
85 12 : char* alreadyExistingProjectionRef = NULL;
86 : char* index_filename_mod;
87 : int bExists;
88 : VSIStatBuf sStatBuf;
89 :
90 : /* Check that we are running against at least GDAL 1.4 */
91 : /* Note to developers : if we use newer API, please change the requirement */
92 12 : if (atoi(GDALVersionInfo("VERSION_NUM")) < 1400)
93 : {
94 0 : fprintf(stderr, "At least, GDAL >= 1.4.0 is required for this version of %s, "
95 : "which was compiled against GDAL %s\n", argv[0], GDAL_RELEASE_NAME);
96 0 : exit(1);
97 : }
98 :
99 12 : GDALAllRegister();
100 12 : OGRRegisterAll();
101 :
102 12 : argc = GDALGeneralCmdLineProcessor( argc, &argv, 0 );
103 12 : if( argc < 1 )
104 0 : exit( -argc );
105 :
106 : /* -------------------------------------------------------------------- */
107 : /* Get commandline arguments other than the GDAL raster filenames. */
108 : /* -------------------------------------------------------------------- */
109 14 : for( i_arg = 1; i_arg < argc; i_arg++ )
110 : {
111 14 : if( EQUAL(argv[i_arg], "--utility_version") )
112 : {
113 2 : printf("%s was compiled against GDAL %s and is running against GDAL %s\n",
114 : argv[0], GDAL_RELEASE_NAME, GDALVersionInfo("RELEASE_NAME"));
115 2 : return 0;
116 : }
117 12 : else if( strcmp(argv[i_arg],"-tileindex") == 0 )
118 : {
119 0 : tile_index = argv[++i_arg];
120 : }
121 12 : else if ( strcmp(argv[i_arg],"-write_absolute_path") == 0 )
122 : {
123 0 : write_absolute_path = TRUE;
124 : }
125 12 : else if ( strcmp(argv[i_arg],"-skip_different_projection") == 0 )
126 : {
127 2 : skip_different_projection = TRUE;
128 : }
129 10 : else if( argv[i_arg][0] == '-' )
130 0 : Usage();
131 10 : else if( index_filename == NULL )
132 : {
133 10 : index_filename = argv[i_arg];
134 10 : i_arg++;
135 10 : break;
136 : }
137 : }
138 :
139 10 : if( index_filename == NULL || i_arg == argc )
140 0 : Usage();
141 :
142 : /* -------------------------------------------------------------------- */
143 : /* Open or create the target shapefile and DBF file. */
144 : /* -------------------------------------------------------------------- */
145 10 : index_filename_mod = CPLStrdup(CPLResetExtension(index_filename, "shp"));
146 :
147 10 : bExists = (VSIStat(index_filename_mod, &sStatBuf) == 0);
148 10 : if (!bExists)
149 : {
150 4 : CPLFree(index_filename_mod);
151 4 : index_filename_mod = CPLStrdup(CPLResetExtension(index_filename, "SHP"));
152 4 : bExists = (VSIStat(index_filename_mod, &sStatBuf) == 0);
153 : }
154 10 : CPLFree(index_filename_mod);
155 :
156 10 : if (bExists)
157 : {
158 6 : hTileIndexDS = OGROpen( index_filename, TRUE, NULL );
159 6 : if (hTileIndexDS != NULL)
160 : {
161 6 : hLayer = OGR_DS_GetLayer(hTileIndexDS, 0);
162 : }
163 : }
164 : else
165 : {
166 : OGRSFDriverH hDriver;
167 4 : const char* pszDriverName = "ESRI Shapefile";
168 :
169 4 : printf( "Creating new index file...\n" );
170 4 : hDriver = OGRGetDriverByName( pszDriverName );
171 4 : if( hDriver == NULL )
172 : {
173 0 : printf( "%s driver not available.\n", pszDriverName );
174 0 : exit( 1 );
175 : }
176 :
177 4 : hTileIndexDS = OGR_Dr_CreateDataSource( hDriver, index_filename, NULL );
178 4 : if (hTileIndexDS)
179 : {
180 4 : char* pszLayerName = CPLStrdup(CPLGetBasename(index_filename));
181 4 : OGRSpatialReferenceH hSpatialRef = NULL;
182 4 : GDALDatasetH hDS = GDALOpen( argv[i_arg], GA_ReadOnly );
183 4 : if (hDS)
184 : {
185 4 : const char* pszWKT = GDALGetProjectionRef(hDS);
186 4 : if (pszWKT != NULL && pszWKT[0] != '\0')
187 : {
188 4 : hSpatialRef = OSRNewSpatialReference(pszWKT);
189 : }
190 4 : GDALClose(hDS);
191 : }
192 :
193 4 : hLayer = OGR_DS_CreateLayer( hTileIndexDS, pszLayerName, hSpatialRef, wkbPolygon, NULL );
194 4 : CPLFree(pszLayerName);
195 4 : if (hSpatialRef)
196 4 : OSRRelease(hSpatialRef);
197 :
198 4 : if (hLayer)
199 : {
200 4 : OGRFieldDefnH hFieldDefn = OGR_Fld_Create( tile_index, OFTString );
201 4 : OGR_Fld_SetWidth( hFieldDefn, 255);
202 4 : OGR_L_CreateField( hLayer, hFieldDefn, TRUE );
203 4 : OGR_Fld_Destroy(hFieldDefn);
204 : }
205 : }
206 : }
207 :
208 10 : if( hTileIndexDS == NULL || hLayer == NULL )
209 : {
210 0 : fprintf( stderr, "Unable to open/create shapefile `%s'.\n",
211 : index_filename );
212 0 : exit(2);
213 : }
214 :
215 10 : hFDefn = OGR_L_GetLayerDefn(hLayer);
216 :
217 10 : for( ti_field = 0; ti_field < OGR_FD_GetFieldCount(hFDefn); ti_field++ )
218 : {
219 10 : OGRFieldDefnH hFieldDefn = OGR_FD_GetFieldDefn( hFDefn, ti_field );
220 10 : if( strcmp(OGR_Fld_GetNameRef(hFieldDefn), tile_index) == 0 )
221 10 : break;
222 : }
223 :
224 10 : if( ti_field == OGR_FD_GetFieldCount(hFDefn) )
225 : {
226 0 : fprintf( stderr, "Unable to find field `%s' in DBF file `%s'.\n",
227 : tile_index, index_filename );
228 0 : exit(2);
229 : }
230 :
231 : /* Load in memory existing file names in SHP */
232 10 : nExistingFiles = OGR_L_GetFeatureCount(hLayer, FALSE);
233 10 : if (nExistingFiles)
234 : {
235 : OGRFeatureH hFeature;
236 6 : existingFilesTab = (char**)CPLMalloc(nExistingFiles * sizeof(char*));
237 26 : for(i=0;i<nExistingFiles;i++)
238 : {
239 20 : hFeature = OGR_L_GetNextFeature(hLayer);
240 20 : existingFilesTab[i] = CPLStrdup(OGR_F_GetFieldAsString( hFeature, ti_field ));
241 20 : if (i == 0)
242 : {
243 6 : GDALDatasetH hDS = GDALOpen(existingFilesTab[i], GA_ReadOnly );
244 6 : if (hDS)
245 : {
246 6 : alreadyExistingProjectionRefValid = TRUE;
247 6 : alreadyExistingProjectionRef = CPLStrdup(GDALGetProjectionRef(hDS));
248 6 : GDALClose(hDS);
249 : }
250 : }
251 20 : OGR_F_Destroy( hFeature );
252 : }
253 : }
254 :
255 10 : if (write_absolute_path)
256 : {
257 0 : current_path = CPLGetCurrentDir();
258 0 : if (current_path == NULL)
259 : {
260 0 : fprintf( stderr, "This system does not support the CPLGetCurrentDir call. "
261 : "The option -write_absolute_path will have no effect\n");
262 0 : write_absolute_path = FALSE;
263 : }
264 : }
265 :
266 : /* -------------------------------------------------------------------- */
267 : /* loop over GDAL files, processing. */
268 : /* -------------------------------------------------------------------- */
269 36 : for( ; i_arg < argc; i_arg++ )
270 : {
271 : GDALDatasetH hDS;
272 : double adfGeoTransform[6];
273 : double adfX[5], adfY[5];
274 : int nXSize, nYSize;
275 : char* fileNameToWrite;
276 : const char* projectionRef;
277 : VSIStatBuf sStatBuf;
278 : int k;
279 : OGRFeatureH hFeature;
280 : OGRGeometryH hPoly, hRing;
281 :
282 : /* Make sure it is a file before building absolute path name */
283 26 : if (write_absolute_path && CPLIsFilenameRelative( argv[i_arg] ) &&
284 0 : VSIStat( argv[i_arg], &sStatBuf ) == 0)
285 : {
286 0 : fileNameToWrite = CPLStrdup(CPLProjectRelativeFilename(current_path, argv[i_arg]));
287 : }
288 : else
289 : {
290 26 : fileNameToWrite = CPLStrdup(argv[i_arg]);
291 : }
292 :
293 : /* Checks that file is not already in tileindex */
294 54 : for(i=0;i<nExistingFiles;i++)
295 : {
296 36 : if (EQUAL(fileNameToWrite, existingFilesTab[i]))
297 : {
298 8 : fprintf(stderr, "File %s is already in tileindex. Skipping it.\n",
299 : fileNameToWrite);
300 8 : break;
301 : }
302 : }
303 26 : if (i != nExistingFiles)
304 : {
305 8 : CPLFree(fileNameToWrite);
306 8 : continue;
307 : }
308 :
309 18 : hDS = GDALOpen( argv[i_arg], GA_ReadOnly );
310 18 : if( hDS == NULL )
311 : {
312 0 : fprintf( stderr, "Unable to open %s, skipping.\n",
313 0 : argv[i_arg] );
314 0 : CPLFree(fileNameToWrite);
315 0 : continue;
316 : }
317 :
318 18 : GDALGetGeoTransform( hDS, adfGeoTransform );
319 18 : if( adfGeoTransform[0] == 0.0
320 0 : && adfGeoTransform[1] == 1.0
321 0 : && adfGeoTransform[3] == 0.0
322 0 : && ABS(adfGeoTransform[5]) == 1.0 )
323 : {
324 0 : fprintf( stderr,
325 : "It appears no georeferencing is available for\n"
326 : "`%s', skipping.\n",
327 0 : argv[i_arg] );
328 0 : GDALClose( hDS );
329 0 : CPLFree(fileNameToWrite);
330 0 : continue;
331 : }
332 :
333 18 : projectionRef = GDALGetProjectionRef(hDS);
334 18 : if (alreadyExistingProjectionRefValid)
335 : {
336 : int projectionRefNotNull, alreadyExistingProjectionRefNotNull;
337 14 : projectionRefNotNull = projectionRef && projectionRef[0];
338 14 : alreadyExistingProjectionRefNotNull = alreadyExistingProjectionRef && alreadyExistingProjectionRef[0];
339 28 : if ((projectionRefNotNull &&
340 : alreadyExistingProjectionRefNotNull &&
341 14 : EQUAL(projectionRef, alreadyExistingProjectionRef) == 0) ||
342 : (projectionRefNotNull != alreadyExistingProjectionRefNotNull))
343 : {
344 4 : fprintf(stderr, "Warning : %s is not using the same projection system as "
345 : "other files in the tileindex. This may cause problems when "
346 2 : "using it in MapServer for example.%s\n", argv[i_arg],
347 : (skip_different_projection) ? " Skipping it" : "");
348 2 : if (skip_different_projection)
349 : {
350 2 : CPLFree(fileNameToWrite);
351 2 : GDALClose( hDS );
352 2 : continue;
353 : }
354 : }
355 : }
356 : else
357 : {
358 4 : alreadyExistingProjectionRefValid = TRUE;
359 4 : alreadyExistingProjectionRef = CPLStrdup(projectionRef);
360 : }
361 :
362 16 : nXSize = GDALGetRasterXSize( hDS );
363 16 : nYSize = GDALGetRasterYSize( hDS );
364 :
365 32 : adfX[0] = adfGeoTransform[0]
366 16 : + 0 * adfGeoTransform[1]
367 16 : + 0 * adfGeoTransform[2];
368 32 : adfY[0] = adfGeoTransform[3]
369 16 : + 0 * adfGeoTransform[4]
370 16 : + 0 * adfGeoTransform[5];
371 :
372 32 : adfX[1] = adfGeoTransform[0]
373 16 : + nXSize * adfGeoTransform[1]
374 32 : + 0 * adfGeoTransform[2];
375 32 : adfY[1] = adfGeoTransform[3]
376 16 : + nXSize * adfGeoTransform[4]
377 32 : + 0 * adfGeoTransform[5];
378 :
379 32 : adfX[2] = adfGeoTransform[0]
380 16 : + nXSize * adfGeoTransform[1]
381 32 : + nYSize * adfGeoTransform[2];
382 32 : adfY[2] = adfGeoTransform[3]
383 16 : + nXSize * adfGeoTransform[4]
384 32 : + nYSize * adfGeoTransform[5];
385 :
386 32 : adfX[3] = adfGeoTransform[0]
387 16 : + 0 * adfGeoTransform[1]
388 16 : + nYSize * adfGeoTransform[2];
389 32 : adfY[3] = adfGeoTransform[3]
390 16 : + 0 * adfGeoTransform[4]
391 16 : + nYSize * adfGeoTransform[5];
392 :
393 32 : adfX[4] = adfGeoTransform[0]
394 16 : + 0 * adfGeoTransform[1]
395 16 : + 0 * adfGeoTransform[2];
396 32 : adfY[4] = adfGeoTransform[3]
397 16 : + 0 * adfGeoTransform[4]
398 16 : + 0 * adfGeoTransform[5];
399 :
400 16 : hFeature = OGR_F_Create( OGR_L_GetLayerDefn( hLayer ) );
401 16 : OGR_F_SetFieldString( hFeature, ti_field, fileNameToWrite );
402 :
403 16 : hPoly = OGR_G_CreateGeometry(wkbPolygon);
404 16 : hRing = OGR_G_CreateGeometry(wkbLinearRing);
405 96 : for(k=0;k<5;k++)
406 80 : OGR_G_SetPoint_2D(hRing, k, adfX[k], adfY[k]);
407 16 : OGR_G_AddGeometryDirectly( hPoly, hRing );
408 16 : OGR_F_SetGeometryDirectly( hFeature, hPoly );
409 :
410 16 : if( OGR_L_CreateFeature( hLayer, hFeature ) != OGRERR_NONE )
411 : {
412 0 : printf( "Failed to create feature in shapefile.\n" );
413 0 : break;
414 : }
415 :
416 16 : OGR_F_Destroy( hFeature );
417 :
418 :
419 16 : CPLFree(fileNameToWrite);
420 :
421 16 : GDALClose( hDS );
422 : }
423 :
424 10 : CPLFree(current_path);
425 :
426 10 : if (nExistingFiles)
427 : {
428 26 : for(i=0;i<nExistingFiles;i++)
429 : {
430 20 : CPLFree(existingFilesTab[i]);
431 : }
432 6 : CPLFree(existingFilesTab);
433 : }
434 10 : CPLFree(alreadyExistingProjectionRef);
435 :
436 10 : OGR_DS_Destroy( hTileIndexDS );
437 :
438 10 : GDALDestroyDriverManager();
439 10 : OGRCleanupAll();
440 10 : CSLDestroy(argv);
441 :
442 10 : exit( 0 );
443 : }
|
