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