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