1 : /******************************************************************************
2 : * $Id: gdal_rasterize.cpp 18164 2009-12-03 19:06:03Z chaitanya $
3 : *
4 : * Project: GDAL Utilities
5 : * Purpose: Rasterize OGR shapes into a GDAL raster.
6 : * Author: Frank Warmerdam <warmerdam@pobox.com>
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2005, Frank Warmerdam <warmerdam@pobox.com>
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.h"
31 : #include "gdal_alg.h"
32 : #include "cpl_conv.h"
33 : #include "ogr_api.h"
34 : #include "ogr_srs_api.h"
35 : #include "cpl_string.h"
36 : #include <vector>
37 :
38 : CPL_CVSID("$Id: gdal_rasterize.cpp 18164 2009-12-03 19:06:03Z chaitanya $");
39 :
40 : /************************************************************************/
41 : /* Usage() */
42 : /************************************************************************/
43 :
44 0 : static void Usage()
45 :
46 : {
47 : printf(
48 : "Usage: gdal_rasterize [-b band] [-i] [-at]\n"
49 : " [-burn value] | [-a attribute_name] [-3d]\n"
50 : // " [-of format_driver] [-co key=value]\n"
51 : // " [-te xmin ymin xmax ymax] [-tr xres yres] [-ts width height]\n"
52 : " [-l layername]* [-where expression] [-sql select_statement]\n"
53 0 : " <src_datasource> <dst_filename>\n" );
54 0 : exit( 1 );
55 : }
56 :
57 : /************************************************************************/
58 : /* InvertGeometries() */
59 : /************************************************************************/
60 :
61 0 : static void InvertGeometries( GDALDatasetH hDstDS,
62 : std::vector<OGRGeometryH> &ahGeometries )
63 :
64 : {
65 : OGRGeometryH hCollection =
66 0 : OGR_G_CreateGeometry( wkbGeometryCollection );
67 :
68 : /* -------------------------------------------------------------------- */
69 : /* Create a ring that is a bit outside the raster dataset. */
70 : /* -------------------------------------------------------------------- */
71 : OGRGeometryH hUniversePoly, hUniverseRing;
72 : double adfGeoTransform[6];
73 0 : int brx = GDALGetRasterXSize( hDstDS ) + 2;
74 0 : int bry = GDALGetRasterYSize( hDstDS ) + 2;
75 :
76 0 : GDALGetGeoTransform( hDstDS, adfGeoTransform );
77 :
78 0 : hUniverseRing = OGR_G_CreateGeometry( wkbLinearRing );
79 :
80 : OGR_G_AddPoint_2D(
81 : hUniverseRing,
82 0 : adfGeoTransform[0] + -2*adfGeoTransform[1] + -2*adfGeoTransform[2],
83 0 : adfGeoTransform[3] + -2*adfGeoTransform[4] + -2*adfGeoTransform[5] );
84 :
85 : OGR_G_AddPoint_2D(
86 : hUniverseRing,
87 0 : adfGeoTransform[0] + brx*adfGeoTransform[1] + -2*adfGeoTransform[2],
88 0 : adfGeoTransform[3] + brx*adfGeoTransform[4] + -2*adfGeoTransform[5] );
89 :
90 : OGR_G_AddPoint_2D(
91 : hUniverseRing,
92 0 : adfGeoTransform[0] + brx*adfGeoTransform[1] + bry*adfGeoTransform[2],
93 0 : adfGeoTransform[3] + brx*adfGeoTransform[4] + bry*adfGeoTransform[5] );
94 :
95 : OGR_G_AddPoint_2D(
96 : hUniverseRing,
97 0 : adfGeoTransform[0] + -2*adfGeoTransform[1] + bry*adfGeoTransform[2],
98 0 : adfGeoTransform[3] + -2*adfGeoTransform[4] + bry*adfGeoTransform[5] );
99 :
100 : OGR_G_AddPoint_2D(
101 : hUniverseRing,
102 0 : adfGeoTransform[0] + -2*adfGeoTransform[1] + -2*adfGeoTransform[2],
103 0 : adfGeoTransform[3] + -2*adfGeoTransform[4] + -2*adfGeoTransform[5] );
104 :
105 0 : hUniversePoly = OGR_G_CreateGeometry( wkbPolygon );
106 0 : OGR_G_AddGeometryDirectly( hUniversePoly, hUniverseRing );
107 :
108 0 : OGR_G_AddGeometryDirectly( hCollection, hUniversePoly );
109 :
110 : /* -------------------------------------------------------------------- */
111 : /* Add the rest of the geometries into our collection. */
112 : /* -------------------------------------------------------------------- */
113 : unsigned int iGeom;
114 :
115 0 : for( iGeom = 0; iGeom < ahGeometries.size(); iGeom++ )
116 0 : OGR_G_AddGeometryDirectly( hCollection, ahGeometries[iGeom] );
117 :
118 0 : ahGeometries.resize(1);
119 0 : ahGeometries[0] = hCollection;
120 0 : }
121 :
122 : /************************************************************************/
123 : /* ProcessLayer() */
124 : /* */
125 : /* Process all the features in a layer selection, collecting */
126 : /* geometries and burn values. */
127 : /************************************************************************/
128 :
129 2 : static void ProcessLayer(
130 : OGRLayerH hSrcLayer,
131 : GDALDatasetH hDstDS, std::vector<int> anBandList,
132 : std::vector<double> &adfBurnValues, int b3D, int bInverse,
133 : const char *pszBurnAttribute, char **papszRasterizeOptions )
134 :
135 : {
136 : /* -------------------------------------------------------------------- */
137 : /* Checkout that SRS are the same. */
138 : /* -------------------------------------------------------------------- */
139 2 : OGRSpatialReferenceH hDstSRS = NULL;
140 2 : if( GDALGetProjectionRef( hDstDS ) != NULL )
141 : {
142 : char *pszProjection;
143 :
144 2 : pszProjection = (char *) GDALGetProjectionRef( hDstDS );
145 :
146 2 : hDstSRS = OSRNewSpatialReference(NULL);
147 2 : if( OSRImportFromWkt( hDstSRS, &pszProjection ) != CE_None )
148 : {
149 1 : OSRDestroySpatialReference(hDstSRS);
150 1 : hDstSRS = NULL;
151 : }
152 : }
153 :
154 2 : OGRSpatialReferenceH hSrcSRS = OGR_L_GetSpatialRef(hSrcLayer);
155 3 : if( hDstSRS != NULL && hSrcSRS != NULL )
156 : {
157 1 : if( OSRIsSame(hSrcSRS, hDstSRS) == FALSE )
158 : {
159 : fprintf(stderr,
160 : "Warning : the output raster dataset and the input vector layer do not have the same SRS.\n"
161 0 : "Results might be incorrect (no on-the-fly reprojection of input data).\n");
162 : }
163 : }
164 1 : else if( hDstSRS != NULL && hSrcSRS == NULL )
165 : {
166 : fprintf(stderr,
167 : "Warning : the output raster dataset has a SRS, but the input vector layer SRS is unknown.\n"
168 0 : "Ensure input vector has the same SRS, otherwise results might be incorrect.\n");
169 : }
170 1 : else if( hDstSRS == NULL && hSrcLayer != NULL )
171 : {
172 : fprintf(stderr,
173 : "Warning : the input vector layer has a SRS, but the output raster dataset SRS is unknown.\n"
174 1 : "Ensure output raster dataset has the same SRS, otherwise results might be incorrect.\n");
175 : }
176 :
177 2 : if( hDstSRS != NULL )
178 : {
179 1 : OSRDestroySpatialReference(hDstSRS);
180 : }
181 :
182 : /* -------------------------------------------------------------------- */
183 : /* Get field index, and check. */
184 : /* -------------------------------------------------------------------- */
185 2 : int iBurnField = -1;
186 :
187 2 : if( pszBurnAttribute )
188 : {
189 : iBurnField = OGR_FD_GetFieldIndex( OGR_L_GetLayerDefn( hSrcLayer ),
190 0 : pszBurnAttribute );
191 0 : if( iBurnField == -1 )
192 : {
193 : printf( "Failed to find field %s on layer %s, skipping.\n",
194 : pszBurnAttribute,
195 0 : OGR_FD_GetName( OGR_L_GetLayerDefn( hSrcLayer ) ) );
196 0 : return;
197 : }
198 : }
199 :
200 : /* -------------------------------------------------------------------- */
201 : /* Collect the geometries from this layer, and build list of */
202 : /* burn values. */
203 : /* -------------------------------------------------------------------- */
204 : OGRFeatureH hFeat;
205 2 : std::vector<OGRGeometryH> ahGeometries;
206 2 : std::vector<double> adfFullBurnValues;
207 :
208 2 : OGR_L_ResetReading( hSrcLayer );
209 :
210 10 : while( (hFeat = OGR_L_GetNextFeature( hSrcLayer )) != NULL )
211 : {
212 : OGRGeometryH hGeom;
213 :
214 6 : if( OGR_F_GetGeometryRef( hFeat ) == NULL )
215 : {
216 1 : OGR_F_Destroy( hFeat );
217 1 : continue;
218 : }
219 :
220 5 : hGeom = OGR_G_Clone( OGR_F_GetGeometryRef( hFeat ) );
221 5 : ahGeometries.push_back( hGeom );
222 :
223 20 : for( unsigned int iBand = 0; iBand < anBandList.size(); iBand++ )
224 : {
225 15 : if( adfBurnValues.size() > 0 )
226 : adfFullBurnValues.push_back(
227 15 : adfBurnValues[MIN(iBand,adfBurnValues.size()-1)] );
228 0 : else if( pszBurnAttribute )
229 : {
230 0 : adfFullBurnValues.push_back( OGR_F_GetFieldAsDouble( hFeat, iBurnField ) );
231 : }
232 : /* I have made the 3D option exclusive to other options since it
233 : can be used to modify the value from "-burn value" or
234 : "-a attribute_name" */
235 15 : if( b3D )
236 : {
237 : // TODO: get geometry "z" value
238 : /* Points and Lines will have their "z" values collected at the
239 : point and line levels respectively. However filled polygons
240 : (GDALdllImageFilledPolygon) can use some help by getting
241 : their "z" values here. */
242 0 : adfFullBurnValues.push_back( 0.0 );
243 : }
244 : }
245 :
246 5 : OGR_F_Destroy( hFeat );
247 : }
248 :
249 : /* -------------------------------------------------------------------- */
250 : /* If we are in inverse mode, we add one extra ring around the */
251 : /* whole dataset to invert the concept of insideness and then */
252 : /* merge everything into one geomtry collection. */
253 : /* -------------------------------------------------------------------- */
254 2 : if( bInverse )
255 : {
256 0 : InvertGeometries( hDstDS, ahGeometries );
257 : }
258 :
259 : /* -------------------------------------------------------------------- */
260 : /* Perform the burn. */
261 : /* -------------------------------------------------------------------- */
262 : GDALRasterizeGeometries( hDstDS, anBandList.size(), &(anBandList[0]),
263 : ahGeometries.size(), &(ahGeometries[0]),
264 : NULL, NULL, &(adfFullBurnValues[0]),
265 : papszRasterizeOptions,
266 2 : GDALTermProgress, NULL );
267 :
268 : /* -------------------------------------------------------------------- */
269 : /* Cleanup geometries. */
270 : /* -------------------------------------------------------------------- */
271 : int iGeom;
272 :
273 7 : for( iGeom = ahGeometries.size()-1; iGeom >= 0; iGeom-- )
274 7 : OGR_G_DestroyGeometry( ahGeometries[iGeom] );
275 : }
276 :
277 : /************************************************************************/
278 : /* main() */
279 : /************************************************************************/
280 :
281 3 : int main( int argc, char ** argv )
282 :
283 : {
284 3 : int i, b3D = FALSE;
285 3 : int bInverse = FALSE;
286 3 : const char *pszSrcFilename = NULL;
287 3 : const char *pszDstFilename = NULL;
288 3 : char **papszLayers = NULL;
289 3 : const char *pszSQL = NULL;
290 3 : const char *pszBurnAttribute = NULL;
291 3 : const char *pszWHERE = NULL;
292 3 : std::vector<int> anBandList;
293 3 : std::vector<double> adfBurnValues;
294 3 : char **papszRasterizeOptions = NULL;
295 :
296 : /* Check that we are running against at least GDAL 1.4 */
297 : /* Note to developers : if we use newer API, please change the requirement */
298 3 : if (atoi(GDALVersionInfo("VERSION_NUM")) < 1400)
299 : {
300 : fprintf(stderr, "At least, GDAL >= 1.4.0 is required for this version of %s, "
301 0 : "which was compiled against GDAL %s\n", argv[0], GDAL_RELEASE_NAME);
302 0 : exit(1);
303 : }
304 :
305 3 : GDALAllRegister();
306 3 : OGRRegisterAll();
307 :
308 3 : argc = GDALGeneralCmdLineProcessor( argc, &argv, 0 );
309 3 : if( argc < 1 )
310 0 : exit( -argc );
311 :
312 : /* -------------------------------------------------------------------- */
313 : /* Parse arguments. */
314 : /* -------------------------------------------------------------------- */
315 22 : for( i = 1; i < argc; i++ )
316 : {
317 20 : if( EQUAL(argv[i], "--utility_version") )
318 : {
319 : printf("%s was compiled against GDAL %s and is running against GDAL %s\n",
320 1 : argv[0], GDAL_RELEASE_NAME, GDALVersionInfo("RELEASE_NAME"));
321 1 : return 0;
322 : }
323 19 : else if( EQUAL(argv[i],"-a") && i < argc-1 )
324 : {
325 0 : pszBurnAttribute = argv[++i];
326 : }
327 25 : else if( EQUAL(argv[i],"-b") && i < argc-1 )
328 : {
329 6 : anBandList.push_back( atoi(argv[++i]) );
330 : }
331 13 : else if( EQUAL(argv[i],"-3d") )
332 : {
333 0 : b3D = TRUE;
334 : papszRasterizeOptions =
335 0 : CSLSetNameValue( papszRasterizeOptions, "BURN_VALUE_FROM", "Z");
336 : }
337 13 : else if( EQUAL(argv[i],"-i") )
338 : {
339 0 : bInverse = TRUE;
340 : }
341 13 : else if( EQUAL(argv[i],"-at") )
342 : {
343 : papszRasterizeOptions =
344 1 : CSLSetNameValue( papszRasterizeOptions, "ALL_TOUCHED", "TRUE" );
345 : }
346 18 : else if( EQUAL(argv[i],"-burn") && i < argc-1 )
347 : {
348 6 : adfBurnValues.push_back( atof(argv[++i]) );
349 : }
350 6 : else if( EQUAL(argv[i],"-where") && i < argc-1 )
351 : {
352 0 : pszWHERE = argv[++i];
353 : }
354 8 : else if( EQUAL(argv[i],"-l") && i < argc-1 )
355 : {
356 2 : papszLayers = CSLAddString( papszLayers, argv[++i] );
357 : }
358 4 : else if( EQUAL(argv[i],"-sql") && i < argc-1 )
359 : {
360 0 : pszSQL = argv[++i];
361 : }
362 4 : else if( pszSrcFilename == NULL )
363 : {
364 2 : pszSrcFilename = argv[i];
365 : }
366 2 : else if( pszDstFilename == NULL )
367 : {
368 2 : pszDstFilename = argv[i];
369 : }
370 : else
371 0 : Usage();
372 : }
373 :
374 2 : if( pszSrcFilename == NULL || pszDstFilename == NULL )
375 : {
376 0 : fprintf( stderr, "Missing source or destination.\n\n" );
377 0 : Usage();
378 : }
379 :
380 2 : if( pszSQL == NULL && papszLayers == NULL )
381 : {
382 0 : fprintf( stderr, "At least one of -l or -sql required.\n\n" );
383 0 : Usage();
384 : }
385 :
386 2 : if( adfBurnValues.size() == 0 && pszBurnAttribute == NULL && !b3D )
387 : {
388 0 : fprintf( stderr, "At least one of -3d, -burn or -a required.\n\n" );
389 0 : Usage();
390 : }
391 :
392 2 : if( anBandList.size() == 0 )
393 0 : anBandList.push_back( 1 );
394 :
395 : /* -------------------------------------------------------------------- */
396 : /* Open source vector dataset. */
397 : /* -------------------------------------------------------------------- */
398 : OGRDataSourceH hSrcDS;
399 :
400 2 : hSrcDS = OGROpen( pszSrcFilename, FALSE, NULL );
401 2 : if( hSrcDS == NULL )
402 : {
403 : fprintf( stderr, "Failed to open feature source: %s\n",
404 0 : pszSrcFilename);
405 0 : exit( 1 );
406 : }
407 :
408 : /* -------------------------------------------------------------------- */
409 : /* Open target raster file. Eventually we will add optional */
410 : /* creation. */
411 : /* -------------------------------------------------------------------- */
412 : GDALDatasetH hDstDS;
413 :
414 2 : hDstDS = GDALOpen( pszDstFilename, GA_Update );
415 2 : if( hDstDS == NULL )
416 0 : exit( 2 );
417 :
418 : /* -------------------------------------------------------------------- */
419 : /* Process SQL request. */
420 : /* -------------------------------------------------------------------- */
421 2 : if( pszSQL != NULL )
422 : {
423 : OGRLayerH hLayer;
424 :
425 0 : hLayer = OGR_DS_ExecuteSQL( hSrcDS, pszSQL, NULL, NULL );
426 0 : if( hLayer != NULL )
427 : {
428 : ProcessLayer( hLayer, hDstDS, anBandList,
429 : adfBurnValues, b3D, bInverse, pszBurnAttribute,
430 0 : papszRasterizeOptions );
431 : }
432 : }
433 :
434 : /* -------------------------------------------------------------------- */
435 : /* Process each layer. */
436 : /* -------------------------------------------------------------------- */
437 2 : int nLayerCount = CSLCount(papszLayers);
438 4 : for( i = 0; i < nLayerCount; i++ )
439 : {
440 2 : OGRLayerH hLayer = OGR_DS_GetLayerByName( hSrcDS, papszLayers[i] );
441 2 : if( hLayer == NULL )
442 : {
443 : fprintf( stderr, "Unable to find layer %s, skipping.\n",
444 0 : papszLayers[i] );
445 0 : continue;
446 : }
447 :
448 2 : if( pszWHERE )
449 : {
450 0 : if( OGR_L_SetAttributeFilter( hLayer, pszWHERE ) != OGRERR_NONE )
451 0 : break;
452 : }
453 :
454 : ProcessLayer( hLayer, hDstDS, anBandList,
455 : adfBurnValues, b3D, bInverse, pszBurnAttribute,
456 2 : papszRasterizeOptions );
457 : }
458 :
459 : /* -------------------------------------------------------------------- */
460 : /* Cleanup */
461 : /* -------------------------------------------------------------------- */
462 :
463 2 : OGR_DS_Destroy( hSrcDS );
464 2 : GDALClose( hDstDS );
465 :
466 2 : CSLDestroy( argv );
467 2 : CSLDestroy( papszRasterizeOptions );
468 2 : CSLDestroy( papszLayers );
469 :
470 2 : GDALDestroyDriverManager();
471 2 : OGRCleanupAll();
472 :
473 2 : return 0;
474 : }
|
