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