1 : /******************************************************************************
2 : * $Id: gdal_rasterize.cpp 19218 2010-03-27 21:38:28Z rouault $
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 19218 2010-03-27 21:38:28Z rouault $");
39 :
40 : /************************************************************************/
41 : /* ArgIsNumeric() */
42 : /************************************************************************/
43 :
44 24 : static int ArgIsNumeric( const char *pszArg )
45 :
46 : {
47 24 : if( pszArg[0] == '-' )
48 12 : pszArg++;
49 :
50 24 : if( *pszArg == '\0' )
51 0 : return FALSE;
52 :
53 72 : while( *pszArg != '\0' )
54 : {
55 36 : if( (*pszArg < '0' || *pszArg > '9') && *pszArg != '.' )
56 12 : return FALSE;
57 24 : pszArg++;
58 : }
59 :
60 12 : return TRUE;
61 : }
62 :
63 : /************************************************************************/
64 : /* Usage() */
65 : /************************************************************************/
66 :
67 0 : static void Usage()
68 :
69 : {
70 : printf(
71 : "Usage: gdal_rasterize [-b band]* [-i] [-at]\n"
72 : " [-burn value]* | [-a attribute_name] [-3d]\n"
73 : " [-l layername]* [-where expression] [-sql select_statement]\n"
74 : " [-of format] [-a_srs srs_def] [-co \"NAME=VALUE\"]*\n"
75 : " [-a_nodata value] [-init value]*\n"
76 : " [-te xmin ymin xmax ymax] [-tr xres yres] [-ts width height]\n"
77 : " [-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/\n"
78 : " CInt16/CInt32/CFloat32/CFloat64}] [-q]\n"
79 0 : " <src_datasource> <dst_filename>\n" );
80 0 : exit( 1 );
81 : }
82 :
83 : /************************************************************************/
84 : /* InvertGeometries() */
85 : /************************************************************************/
86 :
87 : static void InvertGeometries( GDALDatasetH hDstDS,
88 0 : std::vector<OGRGeometryH> &ahGeometries )
89 :
90 : {
91 : OGRGeometryH hCollection =
92 0 : OGR_G_CreateGeometry( wkbGeometryCollection );
93 :
94 : /* -------------------------------------------------------------------- */
95 : /* Create a ring that is a bit outside the raster dataset. */
96 : /* -------------------------------------------------------------------- */
97 : OGRGeometryH hUniversePoly, hUniverseRing;
98 : double adfGeoTransform[6];
99 0 : int brx = GDALGetRasterXSize( hDstDS ) + 2;
100 0 : int bry = GDALGetRasterYSize( hDstDS ) + 2;
101 :
102 0 : GDALGetGeoTransform( hDstDS, adfGeoTransform );
103 :
104 0 : hUniverseRing = OGR_G_CreateGeometry( wkbLinearRing );
105 :
106 : OGR_G_AddPoint_2D(
107 : hUniverseRing,
108 : adfGeoTransform[0] + -2*adfGeoTransform[1] + -2*adfGeoTransform[2],
109 0 : adfGeoTransform[3] + -2*adfGeoTransform[4] + -2*adfGeoTransform[5] );
110 :
111 : OGR_G_AddPoint_2D(
112 : hUniverseRing,
113 : adfGeoTransform[0] + brx*adfGeoTransform[1] + -2*adfGeoTransform[2],
114 0 : adfGeoTransform[3] + brx*adfGeoTransform[4] + -2*adfGeoTransform[5] );
115 :
116 : OGR_G_AddPoint_2D(
117 : hUniverseRing,
118 : adfGeoTransform[0] + brx*adfGeoTransform[1] + bry*adfGeoTransform[2],
119 0 : adfGeoTransform[3] + brx*adfGeoTransform[4] + bry*adfGeoTransform[5] );
120 :
121 : OGR_G_AddPoint_2D(
122 : hUniverseRing,
123 : adfGeoTransform[0] + -2*adfGeoTransform[1] + bry*adfGeoTransform[2],
124 0 : adfGeoTransform[3] + -2*adfGeoTransform[4] + bry*adfGeoTransform[5] );
125 :
126 : OGR_G_AddPoint_2D(
127 : hUniverseRing,
128 : adfGeoTransform[0] + -2*adfGeoTransform[1] + -2*adfGeoTransform[2],
129 0 : adfGeoTransform[3] + -2*adfGeoTransform[4] + -2*adfGeoTransform[5] );
130 :
131 0 : hUniversePoly = OGR_G_CreateGeometry( wkbPolygon );
132 0 : OGR_G_AddGeometryDirectly( hUniversePoly, hUniverseRing );
133 :
134 0 : OGR_G_AddGeometryDirectly( hCollection, hUniversePoly );
135 :
136 : /* -------------------------------------------------------------------- */
137 : /* Add the rest of the geometries into our collection. */
138 : /* -------------------------------------------------------------------- */
139 : unsigned int iGeom;
140 :
141 0 : for( iGeom = 0; iGeom < ahGeometries.size(); iGeom++ )
142 0 : OGR_G_AddGeometryDirectly( hCollection, ahGeometries[iGeom] );
143 :
144 0 : ahGeometries.resize(1);
145 0 : ahGeometries[0] = hCollection;
146 0 : }
147 :
148 : /************************************************************************/
149 : /* ProcessLayer() */
150 : /* */
151 : /* Process all the features in a layer selection, collecting */
152 : /* geometries and burn values. */
153 : /************************************************************************/
154 :
155 : static void ProcessLayer(
156 : OGRLayerH hSrcLayer, int bSRSIsSet,
157 : GDALDatasetH hDstDS, std::vector<int> anBandList,
158 : std::vector<double> &adfBurnValues, int b3D, int bInverse,
159 : const char *pszBurnAttribute, char **papszRasterizeOptions,
160 4 : GDALProgressFunc pfnProgress, void* pProgressData )
161 :
162 : {
163 : /* -------------------------------------------------------------------- */
164 : /* Checkout that SRS are the same. */
165 : /* If -a_srs is specified, skip the test */
166 : /* -------------------------------------------------------------------- */
167 4 : if (!bSRSIsSet)
168 : {
169 3 : OGRSpatialReferenceH hDstSRS = NULL;
170 3 : if( GDALGetProjectionRef( hDstDS ) != NULL )
171 : {
172 : char *pszProjection;
173 :
174 3 : pszProjection = (char *) GDALGetProjectionRef( hDstDS );
175 :
176 3 : hDstSRS = OSRNewSpatialReference(NULL);
177 3 : if( OSRImportFromWkt( hDstSRS, &pszProjection ) != CE_None )
178 : {
179 1 : OSRDestroySpatialReference(hDstSRS);
180 1 : hDstSRS = NULL;
181 : }
182 : }
183 :
184 3 : OGRSpatialReferenceH hSrcSRS = OGR_L_GetSpatialRef(hSrcLayer);
185 5 : if( hDstSRS != NULL && hSrcSRS != NULL )
186 : {
187 2 : if( OSRIsSame(hSrcSRS, hDstSRS) == FALSE )
188 : {
189 : fprintf(stderr,
190 : "Warning : the output raster dataset and the input vector layer do not have the same SRS.\n"
191 0 : "Results might be incorrect (no on-the-fly reprojection of input data).\n");
192 : }
193 : }
194 1 : else if( hDstSRS != NULL && hSrcSRS == NULL )
195 : {
196 : fprintf(stderr,
197 : "Warning : the output raster dataset has a SRS, but the input vector layer SRS is unknown.\n"
198 0 : "Ensure input vector has the same SRS, otherwise results might be incorrect.\n");
199 : }
200 1 : else if( hDstSRS == NULL && hSrcLayer != NULL )
201 : {
202 : fprintf(stderr,
203 : "Warning : the input vector layer has a SRS, but the output raster dataset SRS is unknown.\n"
204 1 : "Ensure output raster dataset has the same SRS, otherwise results might be incorrect.\n");
205 : }
206 :
207 3 : if( hDstSRS != NULL )
208 : {
209 2 : OSRDestroySpatialReference(hDstSRS);
210 : }
211 : }
212 :
213 : /* -------------------------------------------------------------------- */
214 : /* Get field index, and check. */
215 : /* -------------------------------------------------------------------- */
216 4 : int iBurnField = -1;
217 :
218 4 : if( pszBurnAttribute )
219 : {
220 : iBurnField = OGR_FD_GetFieldIndex( OGR_L_GetLayerDefn( hSrcLayer ),
221 0 : pszBurnAttribute );
222 0 : if( iBurnField == -1 )
223 : {
224 : printf( "Failed to find field %s on layer %s, skipping.\n",
225 : pszBurnAttribute,
226 0 : OGR_FD_GetName( OGR_L_GetLayerDefn( hSrcLayer ) ) );
227 0 : return;
228 : }
229 : }
230 :
231 : /* -------------------------------------------------------------------- */
232 : /* Collect the geometries from this layer, and build list of */
233 : /* burn values. */
234 : /* -------------------------------------------------------------------- */
235 : OGRFeatureH hFeat;
236 4 : std::vector<OGRGeometryH> ahGeometries;
237 8 : std::vector<double> adfFullBurnValues;
238 :
239 4 : OGR_L_ResetReading( hSrcLayer );
240 :
241 1250 : while( (hFeat = OGR_L_GetNextFeature( hSrcLayer )) != NULL )
242 : {
243 : OGRGeometryH hGeom;
244 :
245 1242 : if( OGR_F_GetGeometryRef( hFeat ) == NULL )
246 : {
247 1 : OGR_F_Destroy( hFeat );
248 1 : continue;
249 : }
250 :
251 1241 : hGeom = OGR_G_Clone( OGR_F_GetGeometryRef( hFeat ) );
252 1241 : ahGeometries.push_back( hGeom );
253 :
254 2492 : for( unsigned int iBand = 0; iBand < anBandList.size(); iBand++ )
255 : {
256 1251 : if( adfBurnValues.size() > 0 )
257 : adfFullBurnValues.push_back(
258 15 : adfBurnValues[MIN(iBand,adfBurnValues.size()-1)] );
259 1236 : else if( pszBurnAttribute )
260 : {
261 0 : adfFullBurnValues.push_back( OGR_F_GetFieldAsDouble( hFeat, iBurnField ) );
262 : }
263 : /* I have made the 3D option exclusive to other options since it
264 : can be used to modify the value from "-burn value" or
265 : "-a attribute_name" */
266 1251 : if( b3D )
267 : {
268 : // TODO: get geometry "z" value
269 : /* Points and Lines will have their "z" values collected at the
270 : point and line levels respectively. However filled polygons
271 : (GDALdllImageFilledPolygon) can use some help by getting
272 : their "z" values here. */
273 1236 : adfFullBurnValues.push_back( 0.0 );
274 : }
275 : }
276 :
277 1241 : OGR_F_Destroy( hFeat );
278 : }
279 :
280 : /* -------------------------------------------------------------------- */
281 : /* If we are in inverse mode, we add one extra ring around the */
282 : /* whole dataset to invert the concept of insideness and then */
283 : /* merge everything into one geomtry collection. */
284 : /* -------------------------------------------------------------------- */
285 4 : if( bInverse )
286 : {
287 0 : InvertGeometries( hDstDS, ahGeometries );
288 : }
289 :
290 : /* -------------------------------------------------------------------- */
291 : /* Perform the burn. */
292 : /* -------------------------------------------------------------------- */
293 : GDALRasterizeGeometries( hDstDS, anBandList.size(), &(anBandList[0]),
294 : ahGeometries.size(), &(ahGeometries[0]),
295 : NULL, NULL, &(adfFullBurnValues[0]),
296 : papszRasterizeOptions,
297 4 : pfnProgress, pProgressData );
298 :
299 : /* -------------------------------------------------------------------- */
300 : /* Cleanup geometries. */
301 : /* -------------------------------------------------------------------- */
302 : int iGeom;
303 :
304 1245 : for( iGeom = ahGeometries.size()-1; iGeom >= 0; iGeom-- )
305 1245 : OGR_G_DestroyGeometry( ahGeometries[iGeom] );
306 : }
307 :
308 : /************************************************************************/
309 : /* CreateOutputDataset() */
310 : /************************************************************************/
311 :
312 : static
313 : GDALDatasetH CreateOutputDataset(std::vector<OGRLayerH> ahLayers,
314 : OGRSpatialReferenceH hSRS,
315 : int bGotBounds, OGREnvelope sEnvelop,
316 : GDALDriverH hDriver, const char* pszDstFilename,
317 : int nXSize, int nYSize, double dfXRes, double dfYRes,
318 : int nBandCount, GDALDataType eOutputType,
319 : char** papszCreateOptions, std::vector<double> adfInitVals,
320 2 : int bNoDataSet, double dfNoData)
321 : {
322 2 : int bFirstLayer = TRUE;
323 2 : char* pszWKT = NULL;
324 2 : GDALDatasetH hDstDS = NULL;
325 : unsigned int i;
326 :
327 4 : for( i = 0; i < ahLayers.size(); i++ )
328 : {
329 2 : OGRLayerH hLayer = ahLayers[i];
330 :
331 2 : if (!bGotBounds)
332 : {
333 2 : OGREnvelope sLayerEnvelop;
334 :
335 2 : if (OGR_L_GetExtent(hLayer, &sLayerEnvelop, FALSE) != OGRERR_NONE)
336 : {
337 0 : fprintf(stderr, "Cannot get layer extent\n");
338 0 : exit(2);
339 : }
340 :
341 2 : if (bFirstLayer)
342 : {
343 2 : sEnvelop.MinX = sLayerEnvelop.MinX;
344 2 : sEnvelop.MinY = sLayerEnvelop.MinY;
345 2 : sEnvelop.MaxX = sLayerEnvelop.MaxX;
346 2 : sEnvelop.MaxY = sLayerEnvelop.MaxY;
347 :
348 2 : if (hSRS == NULL)
349 1 : hSRS = OGR_L_GetSpatialRef(hLayer);
350 :
351 2 : bFirstLayer = FALSE;
352 : }
353 : else
354 : {
355 0 : sEnvelop.MinX = MIN(sEnvelop.MinX, sLayerEnvelop.MinX);
356 0 : sEnvelop.MinY = MIN(sEnvelop.MinY, sLayerEnvelop.MinY);
357 0 : sEnvelop.MaxX = MAX(sEnvelop.MaxX, sLayerEnvelop.MaxX);
358 0 : sEnvelop.MaxY = MAX(sEnvelop.MaxY, sLayerEnvelop.MaxY);
359 : }
360 : }
361 : else
362 : {
363 0 : if (bFirstLayer)
364 : {
365 0 : if (hSRS == NULL)
366 0 : hSRS = OGR_L_GetSpatialRef(hLayer);
367 :
368 0 : bFirstLayer = FALSE;
369 : }
370 : }
371 : }
372 :
373 : double adfProjection[6];
374 2 : if (dfXRes == 0 && dfYRes == 0)
375 : {
376 1 : dfXRes = (sEnvelop.MaxX - sEnvelop.MinX) / nXSize;
377 1 : dfYRes = (sEnvelop.MaxY - sEnvelop.MinY) / nYSize;
378 : }
379 :
380 2 : adfProjection[0] = sEnvelop.MinX;
381 2 : adfProjection[1] = dfXRes;
382 2 : adfProjection[2] = 0;
383 2 : adfProjection[3] = sEnvelop.MaxY;
384 2 : adfProjection[4] = 0;
385 2 : adfProjection[5] = -dfYRes;
386 :
387 2 : if (nXSize == 0 && nYSize == 0)
388 : {
389 1 : nXSize = (int)(0.5 + (sEnvelop.MaxX - sEnvelop.MinX) / dfXRes);
390 1 : nYSize = (int)(0.5 + (sEnvelop.MaxY - sEnvelop.MinY) / dfYRes);
391 : }
392 :
393 : hDstDS = GDALCreate(hDriver, pszDstFilename, nXSize, nYSize,
394 2 : nBandCount, eOutputType, papszCreateOptions);
395 2 : if (hDstDS == NULL)
396 : {
397 0 : fprintf(stderr, "Cannot create %s\n", pszDstFilename);
398 0 : exit(2);
399 : }
400 :
401 2 : GDALSetGeoTransform(hDstDS, adfProjection);
402 :
403 2 : if (hSRS)
404 2 : OSRExportToWkt(hSRS, &pszWKT);
405 2 : if (pszWKT)
406 2 : GDALSetProjection(hDstDS, pszWKT);
407 2 : CPLFree(pszWKT);
408 :
409 : int iBand;
410 : /*if( nBandCount == 3 || nBandCount == 4 )
411 : {
412 : for(iBand = 0; iBand < nBandCount; iBand++)
413 : {
414 : GDALRasterBandH hBand = GDALGetRasterBand(hDstDS, iBand + 1);
415 : GDALSetRasterColorInterpretation(hBand, (GDALColorInterp)(GCI_RedBand + iBand));
416 : }
417 : }*/
418 :
419 2 : if (bNoDataSet)
420 : {
421 4 : for(iBand = 0; iBand < nBandCount; iBand++)
422 : {
423 2 : GDALRasterBandH hBand = GDALGetRasterBand(hDstDS, iBand + 1);
424 2 : GDALSetRasterNoDataValue(hBand, dfNoData);
425 : }
426 : }
427 :
428 2 : if (adfInitVals.size() != 0)
429 : {
430 0 : for(iBand = 0; iBand < MIN(nBandCount,(int)adfInitVals.size()); iBand++)
431 : {
432 0 : GDALRasterBandH hBand = GDALGetRasterBand(hDstDS, iBand + 1);
433 0 : GDALFillRaster(hBand, adfInitVals[iBand], 0);
434 : }
435 : }
436 :
437 2 : return hDstDS;
438 : }
439 :
440 : /************************************************************************/
441 : /* main() */
442 : /************************************************************************/
443 :
444 5 : int main( int argc, char ** argv )
445 :
446 : {
447 5 : int i, b3D = FALSE;
448 5 : int bInverse = FALSE;
449 5 : const char *pszSrcFilename = NULL;
450 5 : const char *pszDstFilename = NULL;
451 5 : char **papszLayers = NULL;
452 5 : const char *pszSQL = NULL;
453 5 : const char *pszBurnAttribute = NULL;
454 5 : const char *pszWHERE = NULL;
455 5 : std::vector<int> anBandList;
456 10 : std::vector<double> adfBurnValues;
457 5 : char **papszRasterizeOptions = NULL;
458 5 : double dfXRes = 0, dfYRes = 0;
459 5 : int bCreateOutput = FALSE;
460 5 : const char* pszFormat = "GTiff";
461 5 : char **papszCreateOptions = NULL;
462 5 : GDALDriverH hDriver = NULL;
463 5 : GDALDataType eOutputType = GDT_Float64;
464 10 : std::vector<double> adfInitVals;
465 5 : int bNoDataSet = FALSE;
466 5 : double dfNoData = 0;
467 5 : OGREnvelope sEnvelop;
468 5 : int bGotBounds = FALSE;
469 5 : int nXSize = 0, nYSize = 0;
470 5 : int bQuiet = FALSE;
471 5 : GDALProgressFunc pfnProgress = GDALTermProgress;
472 5 : OGRSpatialReferenceH hSRS = NULL;
473 :
474 :
475 : /* Check that we are running against at least GDAL 1.4 */
476 : /* Note to developers : if we use newer API, please change the requirement */
477 5 : if (atoi(GDALVersionInfo("VERSION_NUM")) < 1400)
478 : {
479 : fprintf(stderr, "At least, GDAL >= 1.4.0 is required for this version of %s, "
480 0 : "which was compiled against GDAL %s\n", argv[0], GDAL_RELEASE_NAME);
481 0 : exit(1);
482 : }
483 :
484 5 : GDALAllRegister();
485 5 : OGRRegisterAll();
486 :
487 5 : argc = GDALGeneralCmdLineProcessor( argc, &argv, 0 );
488 5 : if( argc < 1 )
489 0 : exit( -argc );
490 :
491 : /* -------------------------------------------------------------------- */
492 : /* Parse arguments. */
493 : /* -------------------------------------------------------------------- */
494 38 : for( i = 1; i < argc; i++ )
495 : {
496 34 : if( EQUAL(argv[i], "--utility_version") )
497 : {
498 : printf("%s was compiled against GDAL %s and is running against GDAL %s\n",
499 1 : argv[0], GDAL_RELEASE_NAME, GDALVersionInfo("RELEASE_NAME"));
500 1 : return 0;
501 : }
502 34 : else if( EQUAL(argv[i],"-q") || EQUAL(argv[i],"-quiet") )
503 : {
504 1 : bQuiet = TRUE;
505 1 : pfnProgress = GDALDummyProgress;
506 : }
507 32 : else if( EQUAL(argv[i],"-a") && i < argc-1 )
508 : {
509 0 : pszBurnAttribute = argv[++i];
510 : }
511 38 : else if( EQUAL(argv[i],"-b") && i < argc-1 )
512 : {
513 6 : if (strchr(argv[i+1], ' '))
514 : {
515 0 : char** papszTokens = CSLTokenizeString( argv[i+1] );
516 0 : char** papszIter = papszTokens;
517 0 : while(papszIter && *papszIter)
518 : {
519 0 : anBandList.push_back(atoi(*papszIter));
520 0 : papszIter ++;
521 : }
522 0 : CSLDestroy(papszTokens);
523 0 : i += 1;
524 : }
525 : else
526 : {
527 18 : while(i < argc-1 && ArgIsNumeric(argv[i+1]))
528 : {
529 6 : anBandList.push_back(atoi(argv[i+1]));
530 6 : i += 1;
531 : }
532 : }
533 : }
534 26 : else if( EQUAL(argv[i],"-3d") )
535 : {
536 2 : b3D = TRUE;
537 : papszRasterizeOptions =
538 2 : CSLSetNameValue( papszRasterizeOptions, "BURN_VALUE_FROM", "Z");
539 : }
540 24 : else if( EQUAL(argv[i],"-i") )
541 : {
542 0 : bInverse = TRUE;
543 : }
544 24 : else if( EQUAL(argv[i],"-at") )
545 : {
546 : papszRasterizeOptions =
547 1 : CSLSetNameValue( papszRasterizeOptions, "ALL_TOUCHED", "TRUE" );
548 : }
549 29 : else if( EQUAL(argv[i],"-burn") && i < argc-1 )
550 : {
551 6 : if (strchr(argv[i+1], ' '))
552 : {
553 0 : char** papszTokens = CSLTokenizeString( argv[i+1] );
554 0 : char** papszIter = papszTokens;
555 0 : while(papszIter && *papszIter)
556 : {
557 0 : adfBurnValues.push_back(atof(*papszIter));
558 0 : papszIter ++;
559 : }
560 0 : CSLDestroy(papszTokens);
561 0 : i += 1;
562 : }
563 : else
564 : {
565 18 : while(i < argc-1 && ArgIsNumeric(argv[i+1]))
566 : {
567 6 : adfBurnValues.push_back(atof(argv[i+1]));
568 6 : i += 1;
569 : }
570 : }
571 : }
572 17 : else if( EQUAL(argv[i],"-where") && i < argc-1 )
573 : {
574 0 : pszWHERE = argv[++i];
575 : }
576 21 : else if( EQUAL(argv[i],"-l") && i < argc-1 )
577 : {
578 4 : papszLayers = CSLAddString( papszLayers, argv[++i] );
579 : }
580 13 : else if( EQUAL(argv[i],"-sql") && i < argc-1 )
581 : {
582 0 : pszSQL = argv[++i];
583 : }
584 13 : else if( EQUAL(argv[i],"-of") && i < argc-1 )
585 : {
586 0 : pszFormat = argv[++i];
587 0 : bCreateOutput = TRUE;
588 : }
589 13 : else if( EQUAL(argv[i],"-init") && i < argc - 1 )
590 : {
591 0 : if (strchr(argv[i+1], ' '))
592 : {
593 0 : char** papszTokens = CSLTokenizeString( argv[i+1] );
594 0 : char** papszIter = papszTokens;
595 0 : while(papszIter && *papszIter)
596 : {
597 0 : adfInitVals.push_back(atof(*papszIter));
598 0 : papszIter ++;
599 : }
600 0 : CSLDestroy(papszTokens);
601 0 : i += 1;
602 : }
603 : else
604 : {
605 0 : while(i < argc-1 && ArgIsNumeric(argv[i+1]))
606 : {
607 0 : adfInitVals.push_back(atof(argv[i+1]));
608 0 : i += 1;
609 : }
610 : }
611 0 : bCreateOutput = TRUE;
612 : }
613 15 : else if( EQUAL(argv[i],"-a_nodata") && i < argc - 1 )
614 : {
615 2 : dfNoData = atof(argv[i+1]);
616 2 : bNoDataSet = TRUE;
617 2 : i += 1;
618 2 : bCreateOutput = TRUE;
619 : }
620 12 : else if( EQUAL(argv[i],"-a_srs") && i < argc-1 )
621 : {
622 1 : hSRS = OSRNewSpatialReference( NULL );
623 :
624 1 : if( OSRSetFromUserInput(hSRS, argv[i+1]) != OGRERR_NONE )
625 : {
626 : fprintf( stderr, "Failed to process SRS definition: %s\n",
627 0 : argv[i+1] );
628 0 : exit( 1 );
629 : }
630 :
631 1 : i++;
632 1 : bCreateOutput = TRUE;
633 : }
634 :
635 10 : else if( EQUAL(argv[i],"-te") && i < argc - 4 )
636 : {
637 0 : sEnvelop.MinX = atof(argv[++i]);
638 0 : sEnvelop.MinY = atof(argv[++i]);
639 0 : sEnvelop.MaxX = atof(argv[++i]);
640 0 : sEnvelop.MaxY = atof(argv[++i]);
641 0 : bGotBounds = TRUE;
642 0 : bCreateOutput = TRUE;
643 : }
644 10 : else if( EQUAL(argv[i],"-a_ullr") && i < argc - 4 )
645 : {
646 0 : sEnvelop.MinX = atof(argv[++i]);
647 0 : sEnvelop.MaxY = atof(argv[++i]);
648 0 : sEnvelop.MaxX = atof(argv[++i]);
649 0 : sEnvelop.MinY = atof(argv[++i]);
650 0 : bGotBounds = TRUE;
651 0 : bCreateOutput = TRUE;
652 : }
653 10 : else if( EQUAL(argv[i],"-co") && i < argc-1 )
654 : {
655 0 : papszCreateOptions = CSLAddString( papszCreateOptions, argv[++i] );
656 0 : bCreateOutput = TRUE;
657 : }
658 10 : else if( EQUAL(argv[i],"-ot") && i < argc-1 )
659 : {
660 : int iType;
661 :
662 0 : for( iType = 1; iType < GDT_TypeCount; iType++ )
663 : {
664 0 : if( GDALGetDataTypeName((GDALDataType)iType) != NULL
665 : && EQUAL(GDALGetDataTypeName((GDALDataType)iType),
666 : argv[i+1]) )
667 : {
668 0 : eOutputType = (GDALDataType) iType;
669 : }
670 : }
671 :
672 0 : if( eOutputType == GDT_Unknown )
673 : {
674 0 : printf( "Unknown output pixel type: %s\n", argv[i+1] );
675 0 : Usage();
676 0 : GDALDestroyDriverManager();
677 0 : exit( 2 );
678 : }
679 0 : i++;
680 0 : bCreateOutput = TRUE;
681 : }
682 11 : else if( (EQUAL(argv[i],"-ts") || EQUAL(argv[i],"-outsize")) && i < argc-2 )
683 : {
684 1 : nXSize = atoi(argv[++i]);
685 1 : nYSize = atoi(argv[++i]);
686 1 : if (nXSize <= 0 || nYSize <= 0)
687 : {
688 0 : printf( "Wrong value for -outsize parameters\n");
689 0 : Usage();
690 0 : exit( 2 );
691 : }
692 1 : bCreateOutput = TRUE;
693 : }
694 10 : else if( EQUAL(argv[i],"-tr") && i < argc-2 )
695 : {
696 1 : dfXRes = atof(argv[++i]);
697 1 : dfYRes = fabs(atof(argv[++i]));
698 1 : if( dfXRes == 0 || dfYRes == 0 )
699 : {
700 0 : printf( "Wrong value for -tr parameters\n");
701 0 : Usage();
702 0 : exit( 2 );
703 : }
704 1 : bCreateOutput = TRUE;
705 : }
706 :
707 8 : else if( pszSrcFilename == NULL )
708 : {
709 4 : pszSrcFilename = argv[i];
710 : }
711 4 : else if( pszDstFilename == NULL )
712 : {
713 4 : pszDstFilename = argv[i];
714 : }
715 : else
716 0 : Usage();
717 : }
718 :
719 4 : if( pszSrcFilename == NULL || pszDstFilename == NULL )
720 : {
721 0 : fprintf( stderr, "Missing source or destination.\n\n" );
722 0 : Usage();
723 : }
724 :
725 4 : if( pszSQL == NULL && papszLayers == NULL )
726 : {
727 0 : fprintf( stderr, "At least one of -l or -sql required.\n\n" );
728 0 : Usage();
729 : }
730 :
731 4 : if( adfBurnValues.size() == 0 && pszBurnAttribute == NULL && !b3D )
732 : {
733 0 : fprintf( stderr, "At least one of -3d, -burn or -a required.\n\n" );
734 0 : Usage();
735 : }
736 :
737 4 : if( bCreateOutput )
738 : {
739 2 : if( dfXRes == 0 && dfYRes == 0 && nXSize == 0 && nYSize == 0 )
740 : {
741 0 : fprintf( stderr, "'-tr xres yes' or '-ts xsize ysize' is required.\n\n" );
742 0 : Usage();
743 : }
744 :
745 2 : if( anBandList.size() != 0 )
746 : {
747 0 : fprintf( stderr, "-b option cannot be used when creating a GDAL dataset.\n\n" );
748 0 : Usage();
749 : }
750 :
751 2 : int nBandCount = 1;
752 :
753 2 : if (adfBurnValues.size() != 0)
754 0 : nBandCount = adfBurnValues.size();
755 :
756 2 : if ((int)adfInitVals.size() > nBandCount)
757 0 : nBandCount = adfInitVals.size();
758 :
759 2 : if (adfInitVals.size() == 1)
760 : {
761 0 : for(i=1;i<=nBandCount - 1;i++)
762 0 : adfInitVals.push_back( adfInitVals[0] );
763 : }
764 :
765 : int i;
766 4 : for(i=1;i<=nBandCount;i++)
767 2 : anBandList.push_back( i );
768 : }
769 : else
770 : {
771 2 : if( anBandList.size() == 0 )
772 0 : anBandList.push_back( 1 );
773 : }
774 :
775 : /* -------------------------------------------------------------------- */
776 : /* Open source vector dataset. */
777 : /* -------------------------------------------------------------------- */
778 : OGRDataSourceH hSrcDS;
779 :
780 4 : hSrcDS = OGROpen( pszSrcFilename, FALSE, NULL );
781 4 : if( hSrcDS == NULL )
782 : {
783 : fprintf( stderr, "Failed to open feature source: %s\n",
784 0 : pszSrcFilename);
785 0 : exit( 1 );
786 : }
787 :
788 : /* -------------------------------------------------------------------- */
789 : /* Open target raster file. Eventually we will add optional */
790 : /* creation. */
791 : /* -------------------------------------------------------------------- */
792 4 : GDALDatasetH hDstDS = NULL;
793 :
794 4 : if (bCreateOutput)
795 : {
796 : /* -------------------------------------------------------------------- */
797 : /* Find the output driver. */
798 : /* -------------------------------------------------------------------- */
799 2 : hDriver = GDALGetDriverByName( pszFormat );
800 2 : if( hDriver == NULL
801 : || GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) == NULL )
802 : {
803 : int iDr;
804 :
805 : printf( "Output driver `%s' not recognised or does not support\n",
806 0 : pszFormat );
807 : printf( "direct output file creation. The following format drivers are configured\n"
808 0 : "and support direct output:\n" );
809 :
810 0 : for( iDr = 0; iDr < GDALGetDriverCount(); iDr++ )
811 : {
812 0 : GDALDriverH hDriver = GDALGetDriver(iDr);
813 :
814 0 : if( GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL) != NULL )
815 : {
816 : printf( " %s: %s\n",
817 : GDALGetDriverShortName( hDriver ),
818 0 : GDALGetDriverLongName( hDriver ) );
819 : }
820 : }
821 0 : printf( "\n" );
822 0 : exit( 1 );
823 : }
824 : }
825 : else
826 : {
827 2 : hDstDS = GDALOpen( pszDstFilename, GA_Update );
828 2 : if( hDstDS == NULL )
829 0 : exit( 2 );
830 : }
831 :
832 : /* -------------------------------------------------------------------- */
833 : /* Process SQL request. */
834 : /* -------------------------------------------------------------------- */
835 4 : if( pszSQL != NULL )
836 : {
837 : OGRLayerH hLayer;
838 :
839 0 : hLayer = OGR_DS_ExecuteSQL( hSrcDS, pszSQL, NULL, NULL );
840 0 : if( hLayer != NULL )
841 : {
842 0 : if (bCreateOutput)
843 : {
844 0 : std::vector<OGRLayerH> ahLayers;
845 0 : ahLayers.push_back(hLayer);
846 :
847 : hDstDS = CreateOutputDataset(ahLayers, hSRS,
848 : bGotBounds, sEnvelop,
849 : hDriver, pszDstFilename,
850 : nXSize, nYSize, dfXRes, dfYRes,
851 : anBandList.size(), eOutputType,
852 : papszCreateOptions, adfInitVals,
853 0 : bNoDataSet, dfNoData);
854 : }
855 :
856 : ProcessLayer( hLayer, hSRS != NULL, hDstDS, anBandList,
857 : adfBurnValues, b3D, bInverse, pszBurnAttribute,
858 0 : papszRasterizeOptions, pfnProgress, NULL );
859 :
860 0 : OGR_DS_ReleaseResultSet( hSrcDS, hLayer );
861 : }
862 : }
863 :
864 : /* -------------------------------------------------------------------- */
865 : /* Create output file if necessary. */
866 : /* -------------------------------------------------------------------- */
867 4 : int nLayerCount = CSLCount(papszLayers);
868 :
869 4 : if (bCreateOutput && hDstDS == NULL)
870 : {
871 2 : std::vector<OGRLayerH> ahLayers;
872 :
873 4 : for( i = 0; i < nLayerCount; i++ )
874 : {
875 2 : OGRLayerH hLayer = OGR_DS_GetLayerByName( hSrcDS, papszLayers[i] );
876 2 : if( hLayer == NULL )
877 : {
878 0 : continue;
879 : }
880 2 : ahLayers.push_back(hLayer);
881 : }
882 :
883 : hDstDS = CreateOutputDataset(ahLayers, hSRS,
884 : bGotBounds, sEnvelop,
885 : hDriver, pszDstFilename,
886 : nXSize, nYSize, dfXRes, dfYRes,
887 : anBandList.size(), eOutputType,
888 : papszCreateOptions, adfInitVals,
889 2 : bNoDataSet, dfNoData);
890 : }
891 :
892 : /* -------------------------------------------------------------------- */
893 : /* Process each layer. */
894 : /* -------------------------------------------------------------------- */
895 :
896 8 : for( i = 0; i < nLayerCount; i++ )
897 : {
898 4 : OGRLayerH hLayer = OGR_DS_GetLayerByName( hSrcDS, papszLayers[i] );
899 4 : if( hLayer == NULL )
900 : {
901 : fprintf( stderr, "Unable to find layer %s, skipping.\n",
902 0 : papszLayers[i] );
903 0 : continue;
904 : }
905 :
906 4 : if( pszWHERE )
907 : {
908 0 : if( OGR_L_SetAttributeFilter( hLayer, pszWHERE ) != OGRERR_NONE )
909 0 : break;
910 : }
911 :
912 : void *pScaledProgress;
913 : pScaledProgress =
914 : GDALCreateScaledProgress( 0.0, 1.0 * (i + 1) / nLayerCount,
915 4 : pfnProgress, NULL );
916 :
917 : ProcessLayer( hLayer, hSRS != NULL, hDstDS, anBandList,
918 : adfBurnValues, b3D, bInverse, pszBurnAttribute,
919 4 : papszRasterizeOptions, GDALScaledProgress, pScaledProgress );
920 :
921 4 : GDALDestroyScaledProgress( pScaledProgress );
922 : }
923 :
924 : /* -------------------------------------------------------------------- */
925 : /* Cleanup */
926 : /* -------------------------------------------------------------------- */
927 :
928 4 : OGR_DS_Destroy( hSrcDS );
929 4 : GDALClose( hDstDS );
930 :
931 4 : OSRDestroySpatialReference(hSRS);
932 :
933 4 : CSLDestroy( argv );
934 4 : CSLDestroy( papszRasterizeOptions );
935 4 : CSLDestroy( papszLayers );
936 4 : CSLDestroy( papszCreateOptions );
937 :
938 4 : GDALDestroyDriverManager();
939 4 : OGRCleanupAll();
940 :
941 4 : return 0;
942 : }
|
