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