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