1 /******************************************************************************
2 * $Id: gdalrasterize.cpp 28053 2014-12-04 09:31:07Z rouault $
3 *
4 * Project: GDAL
5 * Purpose: Vector rasterization.
6 * Author: Frank Warmerdam, warmerdam@pobox.com
7 *
8 ******************************************************************************
9 * Copyright (c) 2005, Frank Warmerdam <warmerdam@pobox.com>
10 * Copyright (c) 2008-2013, Even Rouault <even dot rouault at mines-paris dot org>
11 *
12 * Permission is hereby granted, free of charge, to any person obtaining a
13 * copy of this software and associated documentation files (the "Software"),
14 * to deal in the Software without restriction, including without limitation
15 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 * and/or sell copies of the Software, and to permit persons to whom the
17 * Software is furnished to do so, subject to the following conditions:
18 *
19 * The above copyright notice and this permission notice shall be included
20 * in all copies or substantial portions of the Software.
21 *
22 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 * DEALINGS IN THE SOFTWARE.
29 ****************************************************************************/
30
31 #include <vector>
32
33 #include "gdal_alg.h"
34 #include "gdal_alg_priv.h"
35 #include "gdal_priv.h"
36 #include "ogr_api.h"
37 #include "ogr_geometry.h"
38 #include "ogr_spatialref.h"
39
40 #ifdef OGR_ENABLED
41 #include "ogrsf_frmts.h"
42 #endif
43
44 /************************************************************************/
45 /* gvBurnScanline() */
46 /************************************************************************/
47
gvBurnScanline(void * pCBData,int nY,int nXStart,int nXEnd,double dfVariant)48 void gvBurnScanline( void *pCBData, int nY, int nXStart, int nXEnd,
49 double dfVariant )
50
51 {
52 GDALRasterizeInfo *psInfo = (GDALRasterizeInfo *) pCBData;
53 int iBand;
54
55 if( nXStart > nXEnd )
56 return;
57
58 CPLAssert( nY >= 0 && nY < psInfo->nYSize );
59 CPLAssert( nXStart <= nXEnd );
60 CPLAssert( nXStart < psInfo->nXSize );
61 CPLAssert( nXEnd >= 0 );
62
63 if( nXStart < 0 )
64 nXStart = 0;
65 if( nXEnd >= psInfo->nXSize )
66 nXEnd = psInfo->nXSize - 1;
67
68 if( psInfo->eType == GDT_Byte )
69 {
70 for( iBand = 0; iBand < psInfo->nBands; iBand++ )
71 {
72 unsigned char *pabyInsert;
73 unsigned char nBurnValue = (unsigned char)
74 ( psInfo->padfBurnValue[iBand] +
75 ( (psInfo->eBurnValueSource == GBV_UserBurnValue)?
76 0 : dfVariant ) );
77
78 pabyInsert = psInfo->pabyChunkBuf
79 + iBand * psInfo->nXSize * psInfo->nYSize
80 + nY * psInfo->nXSize + nXStart;
81
82 if( psInfo->eMergeAlg == GRMA_Add ) {
83 int nPixels = nXEnd - nXStart + 1;
84 while( nPixels-- > 0 )
85 *(pabyInsert++) += nBurnValue;
86 } else {
87 memset( pabyInsert, nBurnValue, nXEnd - nXStart + 1 );
88 }
89 }
90 }
91 else if( psInfo->eType == GDT_Float64 )
92 {
93 for( iBand = 0; iBand < psInfo->nBands; iBand++ )
94 {
95 int nPixels = nXEnd - nXStart + 1;
96 double *padfInsert;
97 double dfBurnValue =
98 ( psInfo->padfBurnValue[iBand] +
99 ( (psInfo->eBurnValueSource == GBV_UserBurnValue)?
100 0 : dfVariant ) );
101
102 padfInsert = ((double *) psInfo->pabyChunkBuf)
103 + iBand * psInfo->nXSize * psInfo->nYSize
104 + nY * psInfo->nXSize + nXStart;
105
106 if( psInfo->eMergeAlg == GRMA_Add ) {
107 while( nPixels-- > 0 )
108 *(padfInsert++) += dfBurnValue;
109 } else {
110 while( nPixels-- > 0 )
111 *(padfInsert++) = dfBurnValue;
112 }
113 }
114 }
115 else {
116 CPLAssert(0);
117 }
118 }
119
120 /************************************************************************/
121 /* gvBurnPoint() */
122 /************************************************************************/
123
gvBurnPoint(void * pCBData,int nY,int nX,double dfVariant)124 void gvBurnPoint( void *pCBData, int nY, int nX, double dfVariant )
125
126 {
127 GDALRasterizeInfo *psInfo = (GDALRasterizeInfo *) pCBData;
128 int iBand;
129
130 CPLAssert( nY >= 0 && nY < psInfo->nYSize );
131 CPLAssert( nX >= 0 && nX < psInfo->nXSize );
132
133 if( psInfo->eType == GDT_Byte )
134 {
135 for( iBand = 0; iBand < psInfo->nBands; iBand++ )
136 {
137 unsigned char *pbyInsert = psInfo->pabyChunkBuf
138 + iBand * psInfo->nXSize * psInfo->nYSize
139 + nY * psInfo->nXSize + nX;
140
141 if( psInfo->eMergeAlg == GRMA_Add ) {
142 *pbyInsert += (unsigned char)( psInfo->padfBurnValue[iBand] +
143 ( (psInfo->eBurnValueSource == GBV_UserBurnValue)?
144 0 : dfVariant ) );
145 } else {
146 *pbyInsert = (unsigned char)( psInfo->padfBurnValue[iBand] +
147 ( (psInfo->eBurnValueSource == GBV_UserBurnValue)?
148 0 : dfVariant ) );
149 }
150 }
151 }
152 else if( psInfo->eType == GDT_Float64 )
153 {
154 for( iBand = 0; iBand < psInfo->nBands; iBand++ )
155 {
156 double *pdfInsert = ((double *) psInfo->pabyChunkBuf)
157 + iBand * psInfo->nXSize * psInfo->nYSize
158 + nY * psInfo->nXSize + nX;
159
160 if( psInfo->eMergeAlg == GRMA_Add ) {
161 *pdfInsert += ( psInfo->padfBurnValue[iBand] +
162 ( (psInfo->eBurnValueSource == GBV_UserBurnValue)?
163 0 : dfVariant ) );
164 } else {
165 *pdfInsert = ( psInfo->padfBurnValue[iBand] +
166 ( (psInfo->eBurnValueSource == GBV_UserBurnValue)?
167 0 : dfVariant ) );
168 }
169 }
170 }
171 else {
172 CPLAssert(0);
173 }
174 }
175
176 /************************************************************************/
177 /* GDALCollectRingsFromGeometry() */
178 /************************************************************************/
179
GDALCollectRingsFromGeometry(OGRGeometry * poShape,std::vector<double> & aPointX,std::vector<double> & aPointY,std::vector<double> & aPointVariant,std::vector<int> & aPartSize,GDALBurnValueSrc eBurnValueSrc)180 static void GDALCollectRingsFromGeometry(
181 OGRGeometry *poShape,
182 std::vector<double> &aPointX, std::vector<double> &aPointY,
183 std::vector<double> &aPointVariant,
184 std::vector<int> &aPartSize, GDALBurnValueSrc eBurnValueSrc)
185
186 {
187 if( poShape == NULL )
188 return;
189
190 OGRwkbGeometryType eFlatType = wkbFlatten(poShape->getGeometryType());
191 int i;
192
193 if ( eFlatType == wkbPoint )
194 {
195 OGRPoint *poPoint = (OGRPoint *) poShape;
196 int nNewCount = aPointX.size() + 1;
197
198 aPointX.reserve( nNewCount );
199 aPointY.reserve( nNewCount );
200 aPointX.push_back( poPoint->getX() );
201 aPointY.push_back( poPoint->getY() );
202 aPartSize.push_back( 1 );
203 if( eBurnValueSrc != GBV_UserBurnValue )
204 {
205 /*switch( eBurnValueSrc )
206 {
207 case GBV_Z:*/
208 aPointVariant.reserve( nNewCount );
209 aPointVariant.push_back( poPoint->getZ() );
210 /*break;
211 case GBV_M:
212 aPointVariant.reserve( nNewCount );
213 aPointVariant.push_back( poPoint->getM() );
214 }*/
215 }
216 }
217 else if ( eFlatType == wkbLineString )
218 {
219 OGRLineString *poLine = (OGRLineString *) poShape;
220 int nCount = poLine->getNumPoints();
221 int nNewCount = aPointX.size() + nCount;
222
223 aPointX.reserve( nNewCount );
224 aPointY.reserve( nNewCount );
225 if( eBurnValueSrc != GBV_UserBurnValue )
226 aPointVariant.reserve( nNewCount );
227 for ( i = nCount - 1; i >= 0; i-- )
228 {
229 aPointX.push_back( poLine->getX(i) );
230 aPointY.push_back( poLine->getY(i) );
231 if( eBurnValueSrc != GBV_UserBurnValue )
232 {
233 /*switch( eBurnValueSrc )
234 {
235 case GBV_Z:*/
236 aPointVariant.push_back( poLine->getZ(i) );
237 /*break;
238 case GBV_M:
239 aPointVariant.push_back( poLine->getM(i) );
240 }*/
241 }
242 }
243 aPartSize.push_back( nCount );
244 }
245 else if ( EQUAL(poShape->getGeometryName(),"LINEARRING") )
246 {
247 OGRLinearRing *poRing = (OGRLinearRing *) poShape;
248 int nCount = poRing->getNumPoints();
249 int nNewCount = aPointX.size() + nCount;
250
251 aPointX.reserve( nNewCount );
252 aPointY.reserve( nNewCount );
253 if( eBurnValueSrc != GBV_UserBurnValue )
254 aPointVariant.reserve( nNewCount );
255 for ( i = nCount - 1; i >= 0; i-- )
256 {
257 aPointX.push_back( poRing->getX(i) );
258 aPointY.push_back( poRing->getY(i) );
259 }
260 if( eBurnValueSrc != GBV_UserBurnValue )
261 {
262 /*switch( eBurnValueSrc )
263 {
264 case GBV_Z:*/
265 aPointVariant.push_back( poRing->getZ(i) );
266 /*break;
267 case GBV_M:
268 aPointVariant.push_back( poRing->getM(i) );
269 }*/
270 }
271 aPartSize.push_back( nCount );
272 }
273 else if( eFlatType == wkbPolygon )
274 {
275 OGRPolygon *poPolygon = (OGRPolygon *) poShape;
276
277 GDALCollectRingsFromGeometry( poPolygon->getExteriorRing(),
278 aPointX, aPointY, aPointVariant,
279 aPartSize, eBurnValueSrc );
280
281 for( i = 0; i < poPolygon->getNumInteriorRings(); i++ )
282 GDALCollectRingsFromGeometry( poPolygon->getInteriorRing(i),
283 aPointX, aPointY, aPointVariant,
284 aPartSize, eBurnValueSrc );
285 }
286
287 else if( eFlatType == wkbMultiPoint
288 || eFlatType == wkbMultiLineString
289 || eFlatType == wkbMultiPolygon
290 || eFlatType == wkbGeometryCollection )
291 {
292 OGRGeometryCollection *poGC = (OGRGeometryCollection *) poShape;
293
294 for( i = 0; i < poGC->getNumGeometries(); i++ )
295 GDALCollectRingsFromGeometry( poGC->getGeometryRef(i),
296 aPointX, aPointY, aPointVariant,
297 aPartSize, eBurnValueSrc );
298 }
299 else
300 {
301 CPLDebug( "GDAL", "Rasterizer ignoring non-polygonal geometry." );
302 }
303 }
304
305 /************************************************************************/
306 /* gv_rasterize_one_shape() */
307 /************************************************************************/
308 static void
gv_rasterize_one_shape(unsigned char * pabyChunkBuf,int nYOff,int nXSize,int nYSize,int nBands,GDALDataType eType,int bAllTouched,OGRGeometry * poShape,double * padfBurnValue,GDALBurnValueSrc eBurnValueSrc,GDALRasterMergeAlg eMergeAlg,GDALTransformerFunc pfnTransformer,void * pTransformArg)309 gv_rasterize_one_shape( unsigned char *pabyChunkBuf, int nYOff,
310 int nXSize, int nYSize,
311 int nBands, GDALDataType eType, int bAllTouched,
312 OGRGeometry *poShape, double *padfBurnValue,
313 GDALBurnValueSrc eBurnValueSrc,
314 GDALRasterMergeAlg eMergeAlg,
315 GDALTransformerFunc pfnTransformer,
316 void *pTransformArg )
317
318 {
319 GDALRasterizeInfo sInfo;
320
321 if (poShape == NULL)
322 return;
323
324 sInfo.nXSize = nXSize;
325 sInfo.nYSize = nYSize;
326 sInfo.nBands = nBands;
327 sInfo.pabyChunkBuf = pabyChunkBuf;
328 sInfo.eType = eType;
329 sInfo.padfBurnValue = padfBurnValue;
330 sInfo.eBurnValueSource = eBurnValueSrc;
331 sInfo.eMergeAlg = eMergeAlg;
332
333 /* -------------------------------------------------------------------- */
334 /* Transform polygon geometries into a set of rings and a part */
335 /* size list. */
336 /* -------------------------------------------------------------------- */
337 std::vector<double> aPointX;
338 std::vector<double> aPointY;
339 std::vector<double> aPointVariant;
340 std::vector<int> aPartSize;
341
342 GDALCollectRingsFromGeometry( poShape, aPointX, aPointY, aPointVariant,
343 aPartSize, eBurnValueSrc );
344
345 /* -------------------------------------------------------------------- */
346 /* Transform points if needed. */
347 /* -------------------------------------------------------------------- */
348 if( pfnTransformer != NULL )
349 {
350 int *panSuccess = (int *) CPLCalloc(sizeof(int),aPointX.size());
351
352 // TODO: we need to add all appropriate error checking at some point.
353 pfnTransformer( pTransformArg, FALSE, aPointX.size(),
354 &(aPointX[0]), &(aPointY[0]), NULL, panSuccess );
355 CPLFree( panSuccess );
356 }
357
358 /* -------------------------------------------------------------------- */
359 /* Shift to account for the buffer offset of this buffer. */
360 /* -------------------------------------------------------------------- */
361 unsigned int i;
362
363 for( i = 0; i < aPointY.size(); i++ )
364 aPointY[i] -= nYOff;
365
366 /* -------------------------------------------------------------------- */
367 /* Perform the rasterization. */
368 /* According to the C++ Standard/23.2.4, elements of a vector are */
369 /* stored in continuous memory block. */
370 /* -------------------------------------------------------------------- */
371
372 // TODO - mloskot: Check if vectors are empty, otherwise it may
373 // lead to undefined behavior by returning non-referencable pointer.
374 // if (!aPointX.empty())
375 // /* fill polygon */
376 // else
377 // /* How to report this problem? */
378 switch ( wkbFlatten(poShape->getGeometryType()) )
379 {
380 case wkbPoint:
381 case wkbMultiPoint:
382 GDALdllImagePoint( sInfo.nXSize, nYSize,
383 aPartSize.size(), &(aPartSize[0]),
384 &(aPointX[0]), &(aPointY[0]),
385 (eBurnValueSrc == GBV_UserBurnValue)?
386 NULL : &(aPointVariant[0]),
387 gvBurnPoint, &sInfo );
388 break;
389 case wkbLineString:
390 case wkbMultiLineString:
391 {
392 if( bAllTouched )
393 GDALdllImageLineAllTouched( sInfo.nXSize, nYSize,
394 aPartSize.size(), &(aPartSize[0]),
395 &(aPointX[0]), &(aPointY[0]),
396 (eBurnValueSrc == GBV_UserBurnValue)?
397 NULL : &(aPointVariant[0]),
398 gvBurnPoint, &sInfo );
399 else
400 GDALdllImageLine( sInfo.nXSize, nYSize,
401 aPartSize.size(), &(aPartSize[0]),
402 &(aPointX[0]), &(aPointY[0]),
403 (eBurnValueSrc == GBV_UserBurnValue)?
404 NULL : &(aPointVariant[0]),
405 gvBurnPoint, &sInfo );
406 }
407 break;
408
409 default:
410 {
411 GDALdllImageFilledPolygon( sInfo.nXSize, nYSize,
412 aPartSize.size(), &(aPartSize[0]),
413 &(aPointX[0]), &(aPointY[0]),
414 (eBurnValueSrc == GBV_UserBurnValue)?
415 NULL : &(aPointVariant[0]),
416 gvBurnScanline, &sInfo );
417 if( bAllTouched )
418 {
419 /* Reverting the variants to the first value because the
420 polygon is filled using the variant from the first point of
421 the first segment. Should be removed when the code to full
422 polygons more appropriately is added. */
423 if(eBurnValueSrc == GBV_UserBurnValue)
424 {
425 GDALdllImageLineAllTouched( sInfo.nXSize, nYSize,
426 aPartSize.size(), &(aPartSize[0]),
427 &(aPointX[0]), &(aPointY[0]),
428 NULL,
429 gvBurnPoint, &sInfo );
430 }
431 else
432 {
433 unsigned int n;
434 for ( i = 0, n = 0; i < aPartSize.size(); i++ )
435 {
436 int j;
437 for ( j = 0; j < aPartSize[i]; j++ )
438 aPointVariant[n++] = aPointVariant[0];
439 }
440
441 GDALdllImageLineAllTouched( sInfo.nXSize, nYSize,
442 aPartSize.size(), &(aPartSize[0]),
443 &(aPointX[0]), &(aPointY[0]),
444 &(aPointVariant[0]),
445 gvBurnPoint, &sInfo );
446 }
447 }
448 }
449 break;
450 }
451 }
452
453 /************************************************************************/
454 /* GDALRasterizeOptions() */
455 /* */
456 /* Recognise a few rasterize options used by all three entry */
457 /* points. */
458 /************************************************************************/
459
GDALRasterizeOptions(char ** papszOptions,int * pbAllTouched,GDALBurnValueSrc * peBurnValueSource,GDALRasterMergeAlg * peMergeAlg)460 static CPLErr GDALRasterizeOptions(char **papszOptions,
461 int *pbAllTouched,
462 GDALBurnValueSrc *peBurnValueSource,
463 GDALRasterMergeAlg *peMergeAlg)
464 {
465 *pbAllTouched = CSLFetchBoolean( papszOptions, "ALL_TOUCHED", FALSE );
466
467 const char *pszOpt = CSLFetchNameValue( papszOptions, "BURN_VALUE_FROM" );
468 *peBurnValueSource = GBV_UserBurnValue;
469 if( pszOpt )
470 {
471 if( EQUAL(pszOpt,"Z"))
472 *peBurnValueSource = GBV_Z;
473 /*else if( EQUAL(pszOpt,"M"))
474 eBurnValueSource = GBV_M;*/
475 else
476 {
477 CPLError( CE_Failure, CPLE_AppDefined,
478 "Unrecognised value '%s' for BURN_VALUE_FROM.",
479 pszOpt );
480 return CE_Failure;
481 }
482 }
483
484 /* -------------------------------------------------------------------- */
485 /* MERGE_ALG=[REPLACE]/ADD */
486 /* -------------------------------------------------------------------- */
487 *peMergeAlg = GRMA_Replace;
488 pszOpt = CSLFetchNameValue( papszOptions, "MERGE_ALG" );
489 if( pszOpt )
490 {
491 if( EQUAL(pszOpt,"ADD"))
492 *peMergeAlg = GRMA_Add;
493 else if( EQUAL(pszOpt,"REPLACE"))
494 *peMergeAlg = GRMA_Replace;
495 else
496 {
497 CPLError( CE_Failure, CPLE_AppDefined,
498 "Unrecognised value '%s' for MERGE_ALG.",
499 pszOpt );
500 return CE_Failure;
501 }
502 }
503
504 return CE_None;
505 }
506
507 /************************************************************************/
508 /* GDALRasterizeGeometries() */
509 /************************************************************************/
510
511 /**
512 * Burn geometries into raster.
513 *
514 * Rasterize a list of geometric objects into a raster dataset. The
515 * geometries are passed as an array of OGRGeometryH handlers.
516 *
517 * If the geometries are in the georferenced coordinates of the raster
518 * dataset, then the pfnTransform may be passed in NULL and one will be
519 * derived internally from the geotransform of the dataset. The transform
520 * needs to transform the geometry locations into pixel/line coordinates
521 * on the raster dataset.
522 *
523 * The output raster may be of any GDAL supported datatype, though currently
524 * internally the burning is done either as GDT_Byte or GDT_Float32. This
525 * may be improved in the future. An explicit list of burn values for
526 * each geometry for each band must be passed in.
527 *
528 * The papszOption list of options currently only supports one option. The
529 * "ALL_TOUCHED" option may be enabled by setting it to "TRUE".
530 *
531 * @param hDS output data, must be opened in update mode.
532 * @param nBandCount the number of bands to be updated.
533 * @param panBandList the list of bands to be updated.
534 * @param nGeomCount the number of geometries being passed in pahGeometries.
535 * @param pahGeometries the array of geometries to burn in.
536 * @param pfnTransformer transformation to apply to geometries to put into
537 * pixel/line coordinates on raster. If NULL a geotransform based one will
538 * be created internally.
539 * @param pTransformArg callback data for transformer.
540 * @param padfGeomBurnValue the array of values to burn into the raster.
541 * There should be nBandCount values for each geometry.
542 * @param papszOptions special options controlling rasterization
543 * <dl>
544 * <dt>"ALL_TOUCHED":</dt> <dd>May be set to TRUE to set all pixels touched
545 * by the line or polygons, not just those whose center is within the polygon
546 * or that are selected by brezenhams line algorithm. Defaults to FALSE.</dd>
547 * <dt>"BURN_VALUE_FROM":</dt> <dd>May be set to "Z" to use the Z values of the
548 * geometries. dfBurnValue is added to this before burning.
549 * Defaults to GDALBurnValueSrc.GBV_UserBurnValue in which case just the
550 * dfBurnValue is burned. This is implemented only for points and lines for
551 * now. The M value may be supported in the future.</dd>
552 * <dt>"MERGE_ALG":</dt> <dd>May be REPLACE (the default) or ADD. REPLACE results in overwriting of value, while ADD adds the new value to the existing raster, suitable for heatmaps for instance.</dd>
553 * </dl>
554 * @param pfnProgress the progress function to report completion.
555 * @param pProgressArg callback data for progress function.
556 *
557 * @return CE_None on success or CE_Failure on error.
558 */
559
GDALRasterizeGeometries(GDALDatasetH hDS,int nBandCount,int * panBandList,int nGeomCount,OGRGeometryH * pahGeometries,GDALTransformerFunc pfnTransformer,void * pTransformArg,double * padfGeomBurnValue,char ** papszOptions,GDALProgressFunc pfnProgress,void * pProgressArg)560 CPLErr GDALRasterizeGeometries( GDALDatasetH hDS,
561 int nBandCount, int *panBandList,
562 int nGeomCount, OGRGeometryH *pahGeometries,
563 GDALTransformerFunc pfnTransformer,
564 void *pTransformArg,
565 double *padfGeomBurnValue,
566 char **papszOptions,
567 GDALProgressFunc pfnProgress,
568 void *pProgressArg )
569
570 {
571 GDALDataType eType;
572 int nYChunkSize, nScanlineBytes;
573 unsigned char *pabyChunkBuf;
574 int iY;
575 GDALDataset *poDS = (GDALDataset *) hDS;
576
577 if( pfnProgress == NULL )
578 pfnProgress = GDALDummyProgress;
579
580 /* -------------------------------------------------------------------- */
581 /* Do some rudimentary arg checking. */
582 /* -------------------------------------------------------------------- */
583 if( nBandCount == 0 || nGeomCount == 0 )
584 {
585 pfnProgress(1.0, "", pProgressArg );
586 return CE_None;
587 }
588
589 // prototype band.
590 GDALRasterBand *poBand = poDS->GetRasterBand( panBandList[0] );
591 if (poBand == NULL)
592 return CE_Failure;
593
594 /* -------------------------------------------------------------------- */
595 /* Options */
596 /* -------------------------------------------------------------------- */
597 int bAllTouched;
598 GDALBurnValueSrc eBurnValueSource;
599 GDALRasterMergeAlg eMergeAlg;
600 if( GDALRasterizeOptions(papszOptions, &bAllTouched,
601 &eBurnValueSource, &eMergeAlg) == CE_Failure) {
602 return CE_Failure;
603 }
604
605 /* -------------------------------------------------------------------- */
606 /* If we have no transformer, assume the geometries are in file */
607 /* georeferenced coordinates, and create a transformer to */
608 /* convert that to pixel/line coordinates. */
609 /* */
610 /* We really just need to apply an affine transform, but for */
611 /* simplicity we use the more general GenImgProjTransformer. */
612 /* -------------------------------------------------------------------- */
613 int bNeedToFreeTransformer = FALSE;
614
615 if( pfnTransformer == NULL )
616 {
617 bNeedToFreeTransformer = TRUE;
618
619 pTransformArg =
620 GDALCreateGenImgProjTransformer( NULL, NULL, hDS, NULL,
621 FALSE, 0.0, 0);
622 pfnTransformer = GDALGenImgProjTransform;
623 }
624
625 /* -------------------------------------------------------------------- */
626 /* Establish a chunksize to operate on. The larger the chunk */
627 /* size the less times we need to make a pass through all the */
628 /* shapes. */
629 /* -------------------------------------------------------------------- */
630 if( poBand->GetRasterDataType() == GDT_Byte )
631 eType = GDT_Byte;
632 else
633 eType = GDT_Float64;
634
635 nScanlineBytes = nBandCount * poDS->GetRasterXSize()
636 * (GDALGetDataTypeSize(eType)/8);
637
638 const char *pszYChunkSize = CSLFetchNameValue(papszOptions, "CHUNKYSIZE");
639 if( pszYChunkSize == NULL || ((nYChunkSize = atoi(pszYChunkSize))) == 0)
640 {
641 nYChunkSize = 10000000 / nScanlineBytes;
642 }
643
644 if( nYChunkSize > poDS->GetRasterYSize() )
645 nYChunkSize = poDS->GetRasterYSize();
646
647 CPLDebug( "GDAL", "Rasterizer operating on %d swaths of %d scanlines.",
648 (poDS->GetRasterYSize()+nYChunkSize-1) / nYChunkSize,
649 nYChunkSize );
650
651 pabyChunkBuf = (unsigned char *) VSIMalloc(nYChunkSize * nScanlineBytes);
652 if( pabyChunkBuf == NULL )
653 {
654 CPLError( CE_Failure, CPLE_OutOfMemory,
655 "Unable to allocate rasterization buffer." );
656 return CE_Failure;
657 }
658
659 /* ==================================================================== */
660 /* Loop over image in designated chunks. */
661 /* ==================================================================== */
662 CPLErr eErr = CE_None;
663
664 pfnProgress( 0.0, NULL, pProgressArg );
665
666 for( iY = 0;
667 iY < poDS->GetRasterYSize() && eErr == CE_None;
668 iY += nYChunkSize )
669 {
670 int nThisYChunkSize;
671 int iShape;
672
673 nThisYChunkSize = nYChunkSize;
674 if( nThisYChunkSize + iY > poDS->GetRasterYSize() )
675 nThisYChunkSize = poDS->GetRasterYSize() - iY;
676
677 eErr =
678 poDS->RasterIO(GF_Read,
679 0, iY, poDS->GetRasterXSize(), nThisYChunkSize,
680 pabyChunkBuf,poDS->GetRasterXSize(),nThisYChunkSize,
681 eType, nBandCount, panBandList,
682 0, 0, 0, NULL );
683 if( eErr != CE_None )
684 break;
685
686 for( iShape = 0; iShape < nGeomCount; iShape++ )
687 {
688 gv_rasterize_one_shape( pabyChunkBuf, iY,
689 poDS->GetRasterXSize(), nThisYChunkSize,
690 nBandCount, eType, bAllTouched,
691 (OGRGeometry *) pahGeometries[iShape],
692 padfGeomBurnValue + iShape*nBandCount,
693 eBurnValueSource, eMergeAlg,
694 pfnTransformer, pTransformArg );
695 }
696
697 eErr =
698 poDS->RasterIO( GF_Write, 0, iY,
699 poDS->GetRasterXSize(), nThisYChunkSize,
700 pabyChunkBuf,
701 poDS->GetRasterXSize(), nThisYChunkSize,
702 eType, nBandCount, panBandList, 0, 0, 0, NULL);
703
704 if( !pfnProgress((iY+nThisYChunkSize)/((double)poDS->GetRasterYSize()),
705 "", pProgressArg ) )
706 {
707 CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
708 eErr = CE_Failure;
709 }
710 }
711
712 /* -------------------------------------------------------------------- */
713 /* cleanup */
714 /* -------------------------------------------------------------------- */
715 VSIFree( pabyChunkBuf );
716
717 if( bNeedToFreeTransformer )
718 GDALDestroyTransformer( pTransformArg );
719
720 return eErr;
721 }
722
723 /************************************************************************/
724 /* GDALRasterizeLayers() */
725 /************************************************************************/
726
727 /**
728 * Burn geometries from the specified list of layers into raster.
729 *
730 * Rasterize all the geometric objects from a list of layers into a raster
731 * dataset. The layers are passed as an array of OGRLayerH handlers.
732 *
733 * If the geometries are in the georferenced coordinates of the raster
734 * dataset, then the pfnTransform may be passed in NULL and one will be
735 * derived internally from the geotransform of the dataset. The transform
736 * needs to transform the geometry locations into pixel/line coordinates
737 * on the raster dataset.
738 *
739 * The output raster may be of any GDAL supported datatype, though currently
740 * internally the burning is done either as GDT_Byte or GDT_Float32. This
741 * may be improved in the future. An explicit list of burn values for
742 * each layer for each band must be passed in.
743 *
744 * @param hDS output data, must be opened in update mode.
745 * @param nBandCount the number of bands to be updated.
746 * @param panBandList the list of bands to be updated.
747 * @param nLayerCount the number of layers being passed in pahLayers array.
748 * @param pahLayers the array of layers to burn in.
749 * @param pfnTransformer transformation to apply to geometries to put into
750 * pixel/line coordinates on raster. If NULL a geotransform based one will
751 * be created internally.
752 * @param pTransformArg callback data for transformer.
753 * @param padfLayerBurnValues the array of values to burn into the raster.
754 * There should be nBandCount values for each layer.
755 * @param papszOptions special options controlling rasterization:
756 * <dl>
757 * <dt>"ATTRIBUTE":</dt> <dd>Identifies an attribute field on the features to be
758 * used for a burn in value. The value will be burned into all output
759 * bands. If specified, padfLayerBurnValues will not be used and can be a NULL
760 * pointer.</dd>
761 * <dt>"CHUNKYSIZE":</dt> <dd>The height in lines of the chunk to operate on.
762 * The larger the chunk size the less times we need to make a pass through all
763 * the shapes. If it is not set or set to zero the default chunk size will be
764 * used. Default size will be estimated based on the GDAL cache buffer size
765 * using formula: cache_size_bytes/scanline_size_bytes, so the chunk will
766 * not exceed the cache.</dd>
767 * <dt>"ALL_TOUCHED":</dt> <dd>May be set to TRUE to set all pixels touched
768 * by the line or polygons, not just those whose center is within the polygon
769 * or that are selected by brezenhams line algorithm. Defaults to FALSE.</dd>
770 * <dt>"BURN_VALUE_FROM":</dt> <dd>May be set to "Z" to use the Z values of the
771 * geometries. The value from padfLayerBurnValues or the attribute field value
772 * is added to this before burning. In default case dfBurnValue is burned as it
773 * is. This is implemented properly only for points and lines for now. Polygons
774 * will be burned using the Z value from the first point. The M value may be
775 * supported in the future.</dd>
776 * <dt>"MERGE_ALG":</dt> <dd>May be REPLACE (the default) or ADD. REPLACE results in overwriting of value, while ADD adds the new value to the existing raster, suitable for heatmaps for instance.</dd>
777 * </dl>
778 * @param pfnProgress the progress function to report completion.
779 * @param pProgressArg callback data for progress function.
780 *
781 * @return CE_None on success or CE_Failure on error.
782 */
783
GDALRasterizeLayers(GDALDatasetH hDS,int nBandCount,int * panBandList,int nLayerCount,OGRLayerH * pahLayers,GDALTransformerFunc pfnTransformer,void * pTransformArg,double * padfLayerBurnValues,char ** papszOptions,GDALProgressFunc pfnProgress,void * pProgressArg)784 CPLErr GDALRasterizeLayers( GDALDatasetH hDS,
785 int nBandCount, int *panBandList,
786 int nLayerCount, OGRLayerH *pahLayers,
787 GDALTransformerFunc pfnTransformer,
788 void *pTransformArg,
789 double *padfLayerBurnValues,
790 char **papszOptions,
791 GDALProgressFunc pfnProgress,
792 void *pProgressArg )
793
794 {
795 #ifndef OGR_ENABLED
796 CPLError(CE_Failure, CPLE_NotSupported, "GDALRasterizeLayers() unimplemented in a non OGR build");
797 return CE_Failure;
798 #else
799 GDALDataType eType;
800 unsigned char *pabyChunkBuf;
801 GDALDataset *poDS = (GDALDataset *) hDS;
802
803 if( pfnProgress == NULL )
804 pfnProgress = GDALDummyProgress;
805
806 /* -------------------------------------------------------------------- */
807 /* Do some rudimentary arg checking. */
808 /* -------------------------------------------------------------------- */
809 if( nBandCount == 0 || nLayerCount == 0 )
810 return CE_None;
811
812 // prototype band.
813 GDALRasterBand *poBand = poDS->GetRasterBand( panBandList[0] );
814 if (poBand == NULL)
815 return CE_Failure;
816
817 /* -------------------------------------------------------------------- */
818 /* Options */
819 /* -------------------------------------------------------------------- */
820 int bAllTouched;
821 GDALBurnValueSrc eBurnValueSource;
822 GDALRasterMergeAlg eMergeAlg;
823 if( GDALRasterizeOptions(papszOptions, &bAllTouched,
824 &eBurnValueSource, &eMergeAlg) == CE_Failure) {
825 return CE_Failure;
826 }
827
828 /* -------------------------------------------------------------------- */
829 /* Establish a chunksize to operate on. The larger the chunk */
830 /* size the less times we need to make a pass through all the */
831 /* shapes. */
832 /* -------------------------------------------------------------------- */
833 int nYChunkSize, nScanlineBytes;
834 const char *pszYChunkSize =
835 CSLFetchNameValue( papszOptions, "CHUNKYSIZE" );
836
837 if( poBand->GetRasterDataType() == GDT_Byte )
838 eType = GDT_Byte;
839 else
840 eType = GDT_Float64;
841
842 nScanlineBytes = nBandCount * poDS->GetRasterXSize()
843 * (GDALGetDataTypeSize(eType)/8);
844
845 if ( pszYChunkSize && ((nYChunkSize = atoi(pszYChunkSize))) != 0 )
846 ;
847 else
848 {
849 GIntBig nYChunkSize64 = GDALGetCacheMax64() / nScanlineBytes;
850 if (nYChunkSize64 > INT_MAX)
851 nYChunkSize = INT_MAX;
852 else
853 nYChunkSize = (int)nYChunkSize64;
854 }
855
856 if( nYChunkSize < 1 )
857 nYChunkSize = 1;
858 if( nYChunkSize > poDS->GetRasterYSize() )
859 nYChunkSize = poDS->GetRasterYSize();
860
861 CPLDebug( "GDAL", "Rasterizer operating on %d swaths of %d scanlines.",
862 (poDS->GetRasterYSize()+nYChunkSize-1) / nYChunkSize,
863 nYChunkSize );
864 pabyChunkBuf = (unsigned char *) VSIMalloc(nYChunkSize * nScanlineBytes);
865 if( pabyChunkBuf == NULL )
866 {
867 CPLError( CE_Failure, CPLE_OutOfMemory,
868 "Unable to allocate rasterization buffer." );
869 return CE_Failure;
870 }
871
872 /* -------------------------------------------------------------------- */
873 /* Read the image once for all layers if user requested to render */
874 /* the whole raster in single chunk. */
875 /* -------------------------------------------------------------------- */
876 if ( nYChunkSize == poDS->GetRasterYSize() )
877 {
878 if ( poDS->RasterIO( GF_Read, 0, 0, poDS->GetRasterXSize(),
879 nYChunkSize, pabyChunkBuf,
880 poDS->GetRasterXSize(), nYChunkSize,
881 eType, nBandCount, panBandList, 0, 0, 0, NULL )
882 != CE_None )
883 {
884 CPLError( CE_Failure, CPLE_OutOfMemory,
885 "Unable to read buffer." );
886 CPLFree( pabyChunkBuf );
887 return CE_Failure;
888 }
889 }
890
891 /* ==================================================================== */
892 /* Read the specified layers transfoming and rasterizing */
893 /* geometries. */
894 /* ==================================================================== */
895 CPLErr eErr = CE_None;
896 int iLayer;
897 const char *pszBurnAttribute =
898 CSLFetchNameValue( papszOptions, "ATTRIBUTE" );
899
900 pfnProgress( 0.0, NULL, pProgressArg );
901
902 for( iLayer = 0; iLayer < nLayerCount; iLayer++ )
903 {
904 int iBurnField = -1;
905 double *padfBurnValues = NULL;
906 OGRLayer *poLayer = (OGRLayer *) pahLayers[iLayer];
907
908 if ( !poLayer )
909 {
910 CPLError( CE_Warning, CPLE_AppDefined,
911 "Layer element number %d is NULL, skipping.\n", iLayer );
912 continue;
913 }
914
915 /* -------------------------------------------------------------------- */
916 /* If the layer does not contain any features just skip it. */
917 /* Do not force the feature count, so if driver doesn't know */
918 /* exact number of features, go down the normal way. */
919 /* -------------------------------------------------------------------- */
920 if ( poLayer->GetFeatureCount(FALSE) == 0 )
921 continue;
922
923 if ( pszBurnAttribute )
924 {
925 iBurnField =
926 poLayer->GetLayerDefn()->GetFieldIndex( pszBurnAttribute );
927 if ( iBurnField == -1 )
928 {
929 CPLError( CE_Warning, CPLE_AppDefined,
930 "Failed to find field %s on layer %s, skipping.\n",
931 pszBurnAttribute,
932 poLayer->GetLayerDefn()->GetName() );
933 continue;
934 }
935 }
936 else
937 padfBurnValues = padfLayerBurnValues + iLayer * nBandCount;
938
939 /* -------------------------------------------------------------------- */
940 /* If we have no transformer, create the one from input file */
941 /* projection. Note that each layer can be georefernced */
942 /* separately. */
943 /* -------------------------------------------------------------------- */
944 int bNeedToFreeTransformer = FALSE;
945
946 if( pfnTransformer == NULL )
947 {
948 char *pszProjection = NULL;
949 bNeedToFreeTransformer = TRUE;
950
951 OGRSpatialReference *poSRS = poLayer->GetSpatialRef();
952 if ( !poSRS )
953 {
954 CPLError( CE_Warning, CPLE_AppDefined,
955 "Failed to fetch spatial reference on layer %s "
956 "to build transformer, assuming matching coordinate systems.\n",
957 poLayer->GetLayerDefn()->GetName() );
958 }
959 else
960 poSRS->exportToWkt( &pszProjection );
961
962 pTransformArg =
963 GDALCreateGenImgProjTransformer( NULL, pszProjection,
964 hDS, NULL, FALSE, 0.0, 0 );
965 pfnTransformer = GDALGenImgProjTransform;
966
967 CPLFree( pszProjection );
968 }
969
970 OGRFeature *poFeat;
971
972 poLayer->ResetReading();
973
974 /* -------------------------------------------------------------------- */
975 /* Loop over image in designated chunks. */
976 /* -------------------------------------------------------------------- */
977 int iY;
978 for( iY = 0;
979 iY < poDS->GetRasterYSize() && eErr == CE_None;
980 iY += nYChunkSize )
981 {
982 int nThisYChunkSize;
983
984 nThisYChunkSize = nYChunkSize;
985 if( nThisYChunkSize + iY > poDS->GetRasterYSize() )
986 nThisYChunkSize = poDS->GetRasterYSize() - iY;
987
988 // Only re-read image if not a single chunk is being rendered
989 if ( nYChunkSize < poDS->GetRasterYSize() )
990 {
991 eErr =
992 poDS->RasterIO( GF_Read, 0, iY,
993 poDS->GetRasterXSize(), nThisYChunkSize,
994 pabyChunkBuf,
995 poDS->GetRasterXSize(), nThisYChunkSize,
996 eType, nBandCount, panBandList, 0, 0, 0, NULL );
997 if( eErr != CE_None )
998 break;
999 }
1000
1001 double *padfAttrValues = (double *) VSIMalloc(sizeof(double) * nBandCount);
1002 while( (poFeat = poLayer->GetNextFeature()) != NULL )
1003 {
1004 OGRGeometry *poGeom = poFeat->GetGeometryRef();
1005
1006 if ( pszBurnAttribute )
1007 {
1008 int iBand;
1009 double dfAttrValue;
1010
1011 dfAttrValue = poFeat->GetFieldAsDouble( iBurnField );
1012 for (iBand = 0 ; iBand < nBandCount ; iBand++)
1013 padfAttrValues[iBand] = dfAttrValue;
1014
1015 padfBurnValues = padfAttrValues;
1016 }
1017
1018 gv_rasterize_one_shape( pabyChunkBuf, iY,
1019 poDS->GetRasterXSize(),
1020 nThisYChunkSize,
1021 nBandCount, eType, bAllTouched, poGeom,
1022 padfBurnValues, eBurnValueSource,
1023 eMergeAlg,
1024 pfnTransformer, pTransformArg );
1025
1026 delete poFeat;
1027 }
1028 VSIFree( padfAttrValues );
1029
1030 // Only write image if not a single chunk is being rendered
1031 if ( nYChunkSize < poDS->GetRasterYSize() )
1032 {
1033 eErr =
1034 poDS->RasterIO( GF_Write, 0, iY,
1035 poDS->GetRasterXSize(), nThisYChunkSize,
1036 pabyChunkBuf,
1037 poDS->GetRasterXSize(), nThisYChunkSize,
1038 eType, nBandCount, panBandList, 0, 0, 0, NULL );
1039 }
1040
1041 poLayer->ResetReading();
1042
1043 if( !pfnProgress((iY+nThisYChunkSize)/((double)poDS->GetRasterYSize()),
1044 "", pProgressArg) )
1045 {
1046 CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
1047 eErr = CE_Failure;
1048 }
1049 }
1050
1051 if ( bNeedToFreeTransformer )
1052 {
1053 GDALDestroyTransformer( pTransformArg );
1054 pTransformArg = NULL;
1055 pfnTransformer = NULL;
1056 }
1057 }
1058
1059 /* -------------------------------------------------------------------- */
1060 /* Write out the image once for all layers if user requested */
1061 /* to render the whole raster in single chunk. */
1062 /* -------------------------------------------------------------------- */
1063 if ( nYChunkSize == poDS->GetRasterYSize() )
1064 {
1065 poDS->RasterIO( GF_Write, 0, 0,
1066 poDS->GetRasterXSize(), nYChunkSize,
1067 pabyChunkBuf,
1068 poDS->GetRasterXSize(), nYChunkSize,
1069 eType, nBandCount, panBandList, 0, 0, 0, NULL );
1070 }
1071
1072 /* -------------------------------------------------------------------- */
1073 /* cleanup */
1074 /* -------------------------------------------------------------------- */
1075 VSIFree( pabyChunkBuf );
1076
1077 return eErr;
1078 #endif /* def OGR_ENABLED */
1079 }
1080
1081 /************************************************************************/
1082 /* GDALRasterizeLayersBuf() */
1083 /************************************************************************/
1084
1085 /**
1086 * Burn geometries from the specified list of layer into raster.
1087 *
1088 * Rasterize all the geometric objects from a list of layers into supplied
1089 * raster buffer. The layers are passed as an array of OGRLayerH handlers.
1090 *
1091 * If the geometries are in the georferenced coordinates of the raster
1092 * dataset, then the pfnTransform may be passed in NULL and one will be
1093 * derived internally from the geotransform of the dataset. The transform
1094 * needs to transform the geometry locations into pixel/line coordinates
1095 * of the target raster.
1096 *
1097 * The output raster may be of any GDAL supported datatype, though currently
1098 * internally the burning is done either as GDT_Byte or GDT_Float32. This
1099 * may be improved in the future.
1100 *
1101 * @param pData pointer to the output data array.
1102 *
1103 * @param nBufXSize width of the output data array in pixels.
1104 *
1105 * @param nBufYSize height of the output data array in pixels.
1106 *
1107 * @param eBufType data type of the output data array.
1108 *
1109 * @param nPixelSpace The byte offset from the start of one pixel value in
1110 * pData to the start of the next pixel value within a scanline. If defaulted
1111 * (0) the size of the datatype eBufType is used.
1112 *
1113 * @param nLineSpace The byte offset from the start of one scanline in
1114 * pData to the start of the next. If defaulted the size of the datatype
1115 * eBufType * nBufXSize is used.
1116 *
1117 * @param nLayerCount the number of layers being passed in pahLayers array.
1118 *
1119 * @param pahLayers the array of layers to burn in.
1120 *
1121 * @param pszDstProjection WKT defining the coordinate system of the target
1122 * raster.
1123 *
1124 * @param padfDstGeoTransform geotransformation matrix of the target raster.
1125 *
1126 * @param pfnTransformer transformation to apply to geometries to put into
1127 * pixel/line coordinates on raster. If NULL a geotransform based one will
1128 * be created internally.
1129 *
1130 * @param pTransformArg callback data for transformer.
1131 *
1132 * @param dfBurnValue the value to burn into the raster.
1133 *
1134 * @param papszOptions special options controlling rasterization:
1135 * <dl>
1136 * <dt>"ATTRIBUTE":</dt> <dd>Identifies an attribute field on the features to be
1137 * used for a burn in value. The value will be burned into all output
1138 * bands. If specified, padfLayerBurnValues will not be used and can be a NULL
1139 * pointer.</dd>
1140 * <dt>"ALL_TOUCHED":</dt> <dd>May be set to TRUE to set all pixels touched
1141 * by the line or polygons, not just those whose center is within the polygon
1142 * or that are selected by brezenhams line algorithm. Defaults to FALSE.</dd>
1143 * </dl>
1144 * <dt>"BURN_VALUE_FROM":</dt> <dd>May be set to "Z" to use
1145 * the Z values of the geometries. dfBurnValue or the attribute field value is
1146 * added to this before burning. In default case dfBurnValue is burned as it
1147 * is. This is implemented properly only for points and lines for now. Polygons
1148 * will be burned using the Z value from the first point. The M value may
1149 * be supported in the future.</dd>
1150 * <dt>"MERGE_ALG":</dt> <dd>May be REPLACE (the default) or ADD. REPLACE results in overwriting of value, while ADD adds the new value to the existing raster, suitable for heatmaps for instance.</dd>
1151 * </dl>
1152 *
1153 * @param pfnProgress the progress function to report completion.
1154 *
1155 * @param pProgressArg callback data for progress function.
1156 *
1157 *
1158 * @return CE_None on success or CE_Failure on error.
1159 */
1160
GDALRasterizeLayersBuf(void * pData,int nBufXSize,int nBufYSize,GDALDataType eBufType,int nPixelSpace,int nLineSpace,int nLayerCount,OGRLayerH * pahLayers,const char * pszDstProjection,double * padfDstGeoTransform,GDALTransformerFunc pfnTransformer,void * pTransformArg,double dfBurnValue,char ** papszOptions,GDALProgressFunc pfnProgress,void * pProgressArg)1161 CPLErr GDALRasterizeLayersBuf( void *pData, int nBufXSize, int nBufYSize,
1162 GDALDataType eBufType,
1163 int nPixelSpace, int nLineSpace,
1164 int nLayerCount, OGRLayerH *pahLayers,
1165 const char *pszDstProjection,
1166 double *padfDstGeoTransform,
1167 GDALTransformerFunc pfnTransformer,
1168 void *pTransformArg, double dfBurnValue,
1169 char **papszOptions,
1170 GDALProgressFunc pfnProgress,
1171 void *pProgressArg )
1172
1173 {
1174 #ifndef OGR_ENABLED
1175 CPLError(CE_Failure, CPLE_NotSupported, "GDALRasterizeLayersBuf() unimplemented in a non OGR build");
1176 return CE_Failure;
1177 #else
1178 /* -------------------------------------------------------------------- */
1179 /* If pixel and line spaceing are defaulted assign reasonable */
1180 /* value assuming a packed buffer. */
1181 /* -------------------------------------------------------------------- */
1182 if( nPixelSpace == 0 )
1183 nPixelSpace = GDALGetDataTypeSize( eBufType ) / 8;
1184
1185 if( nLineSpace == 0 )
1186 nLineSpace = nPixelSpace * nBufXSize;
1187
1188 if( pfnProgress == NULL )
1189 pfnProgress = GDALDummyProgress;
1190
1191 /* -------------------------------------------------------------------- */
1192 /* Do some rudimentary arg checking. */
1193 /* -------------------------------------------------------------------- */
1194 if( nLayerCount == 0 )
1195 return CE_None;
1196
1197 /* -------------------------------------------------------------------- */
1198 /* Options */
1199 /* -------------------------------------------------------------------- */
1200 int bAllTouched;
1201 GDALBurnValueSrc eBurnValueSource;
1202 GDALRasterMergeAlg eMergeAlg;
1203 if( GDALRasterizeOptions(papszOptions, &bAllTouched,
1204 &eBurnValueSource, &eMergeAlg) == CE_Failure) {
1205 return CE_Failure;
1206 }
1207
1208 /* ==================================================================== */
1209 /* Read thes pecified layers transfoming and rasterizing */
1210 /* geometries. */
1211 /* ==================================================================== */
1212 CPLErr eErr = CE_None;
1213 int iLayer;
1214 const char *pszBurnAttribute =
1215 CSLFetchNameValue( papszOptions, "ATTRIBUTE" );
1216
1217 pfnProgress( 0.0, NULL, pProgressArg );
1218
1219 for( iLayer = 0; iLayer < nLayerCount; iLayer++ )
1220 {
1221 int iBurnField = -1;
1222 OGRLayer *poLayer = (OGRLayer *) pahLayers[iLayer];
1223
1224 if ( !poLayer )
1225 {
1226 CPLError( CE_Warning, CPLE_AppDefined,
1227 "Layer element number %d is NULL, skipping.\n", iLayer );
1228 continue;
1229 }
1230
1231 /* -------------------------------------------------------------------- */
1232 /* If the layer does not contain any features just skip it. */
1233 /* Do not force the feature count, so if driver doesn't know */
1234 /* exact number of features, go down the normal way. */
1235 /* -------------------------------------------------------------------- */
1236 if ( poLayer->GetFeatureCount(FALSE) == 0 )
1237 continue;
1238
1239 if ( pszBurnAttribute )
1240 {
1241 iBurnField =
1242 poLayer->GetLayerDefn()->GetFieldIndex( pszBurnAttribute );
1243 if ( iBurnField == -1 )
1244 {
1245 CPLError( CE_Warning, CPLE_AppDefined,
1246 "Failed to find field %s on layer %s, skipping.\n",
1247 pszBurnAttribute,
1248 poLayer->GetLayerDefn()->GetName() );
1249 continue;
1250 }
1251 }
1252
1253 /* -------------------------------------------------------------------- */
1254 /* If we have no transformer, create the one from input file */
1255 /* projection. Note that each layer can be georefernced */
1256 /* separately. */
1257 /* -------------------------------------------------------------------- */
1258 int bNeedToFreeTransformer = FALSE;
1259
1260 if( pfnTransformer == NULL )
1261 {
1262 char *pszProjection = NULL;
1263 bNeedToFreeTransformer = TRUE;
1264
1265 OGRSpatialReference *poSRS = poLayer->GetSpatialRef();
1266 if ( !poSRS )
1267 {
1268 CPLError( CE_Warning, CPLE_AppDefined,
1269 "Failed to fetch spatial reference on layer %s "
1270 "to build transformer, assuming matching coordinate systems.\n",
1271 poLayer->GetLayerDefn()->GetName() );
1272 }
1273 else
1274 poSRS->exportToWkt( &pszProjection );
1275
1276 pTransformArg =
1277 GDALCreateGenImgProjTransformer3( pszProjection, NULL,
1278 pszDstProjection,
1279 padfDstGeoTransform );
1280 pfnTransformer = GDALGenImgProjTransform;
1281
1282 CPLFree( pszProjection );
1283 }
1284
1285 OGRFeature *poFeat;
1286
1287 poLayer->ResetReading();
1288
1289 while( (poFeat = poLayer->GetNextFeature()) != NULL )
1290 {
1291 OGRGeometry *poGeom = poFeat->GetGeometryRef();
1292
1293 if ( pszBurnAttribute )
1294 dfBurnValue = poFeat->GetFieldAsDouble( iBurnField );
1295
1296 gv_rasterize_one_shape( (unsigned char *) pData, 0,
1297 nBufXSize, nBufYSize,
1298 1, eBufType, bAllTouched, poGeom,
1299 &dfBurnValue, eBurnValueSource, eMergeAlg,
1300 pfnTransformer, pTransformArg );
1301
1302 delete poFeat;
1303 }
1304
1305 poLayer->ResetReading();
1306
1307 if( !pfnProgress(1, "", pProgressArg) )
1308 {
1309 CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
1310 eErr = CE_Failure;
1311 }
1312
1313 if ( bNeedToFreeTransformer )
1314 {
1315 GDALDestroyTransformer( pTransformArg );
1316 pTransformArg = NULL;
1317 pfnTransformer = NULL;
1318 }
1319 }
1320
1321 return eErr;
1322 #endif /* def OGR_ENABLED */
1323 }
1324