Ticket #3202: gdalrasterize.cpp

File gdalrasterize.cpp, 32.2 KB (added by xenononarcticus, 7 years ago)
Line 
1/******************************************************************************
2 * $Id: gdalrasterize.cpp 15665 2008-10-31 18:18:35Z warmerdam $
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
43typedef struct {
44    unsigned char * pabyChunkBuf;
45    int nXSize;
46    int nYSize;
47    int nBands;
48    GDALDataType eType;
49    double *padfBurnValue;
50        int bBurnValueSource; // 0 indicates burn value comes from padfBurnValue, 1 indicates value from Variant (Z, or possibly M in the future)
51} GDALRasterizeInfo;
52
53/************************************************************************/
54/*                           gvBurnScanline()                           */
55/************************************************************************/
56
57void gvBurnScanline( void *pCBData, int nY, int nXStart, int nXEnd )
58
59{
60    GDALRasterizeInfo *psInfo = (GDALRasterizeInfo *) pCBData;
61    int iBand;
62
63    if( nXStart > nXEnd )
64        return;
65
66    CPLAssert( nY >= 0 && nY < psInfo->nYSize );
67    CPLAssert( nXStart <= nXEnd );
68    CPLAssert( nXStart < psInfo->nXSize );
69    CPLAssert( nXEnd >= 0 );
70
71    if( nXStart < 0 )
72        nXStart = 0;
73    if( nXEnd >= psInfo->nXSize )
74        nXEnd = psInfo->nXSize - 1;
75
76    if( psInfo->eType == GDT_Byte )
77    {
78        for( iBand = 0; iBand < psInfo->nBands; iBand++ )
79        {
80            unsigned char *pabyInsert;
81            unsigned char nBurnValue = (unsigned char) 
82                psInfo->padfBurnValue[iBand];
83           
84            pabyInsert = psInfo->pabyChunkBuf
85                + iBand * psInfo->nXSize * psInfo->nYSize
86                + nY * psInfo->nXSize + nXStart;
87               
88            memset( pabyInsert, nBurnValue, nXEnd - nXStart + 1 );
89        }
90    }
91    else
92    {
93        for( iBand = 0; iBand < psInfo->nBands; iBand++ )
94        {
95            int nPixels = nXEnd - nXStart + 1;
96            float   *pafInsert;
97            float   fBurnValue = (float) psInfo->padfBurnValue[iBand];
98           
99            pafInsert = ((float *) psInfo->pabyChunkBuf) 
100                + iBand * psInfo->nXSize * psInfo->nYSize
101                + nY * psInfo->nXSize + nXStart;
102
103            while( nPixels-- > 0 )
104                *(pafInsert++) = fBurnValue;
105        }
106    }
107}
108
109/************************************************************************/
110/*                            gvBurnPoint()                             */
111/************************************************************************/
112
113void gvBurnPoint( void *pCBData, int nY, int nX, double dfVariant )
114
115{
116    GDALRasterizeInfo *psInfo = (GDALRasterizeInfo *) pCBData;
117    int iBand;
118
119    CPLAssert( nY >= 0 && nY < psInfo->nYSize );
120    CPLAssert( nX >= 0 && nX < psInfo->nXSize );
121
122    if( psInfo->eType == GDT_Byte )
123    {
124        for( iBand = 0; iBand < psInfo->nBands; iBand++ )
125        {
126            unsigned char *pbyInsert = psInfo->pabyChunkBuf
127                + iBand * psInfo->nXSize * psInfo->nYSize
128                + nY * psInfo->nXSize + nX;
129
130                        if(psInfo->bBurnValueSource) // should we include the Variant as our Burn Value source?
131                        {
132                                *pbyInsert = (unsigned char) dfVariant + (unsigned char)psInfo->padfBurnValue[iBand];
133                        } // if
134                        else
135                        {
136                                *pbyInsert = (unsigned char) psInfo->padfBurnValue[iBand];
137                        } // else
138        }
139    }
140    else
141    {
142        for( iBand = 0; iBand < psInfo->nBands; iBand++ )
143        {
144            float   *pfInsert = ((float *) psInfo->pabyChunkBuf) 
145                + iBand * psInfo->nXSize * psInfo->nYSize
146                + nY * psInfo->nXSize + nX;
147
148                        if(psInfo->bBurnValueSource) // should we use the Variant as our Burn Value source?
149                        {
150                                *pfInsert = (float) dfVariant + (float) psInfo->padfBurnValue[iBand];
151                        } // if
152                        else
153                        {
154                                *pfInsert = (float) psInfo->padfBurnValue[iBand];
155                        } // else
156        }
157    }
158}
159
160/************************************************************************/
161/*                    GDALCollectRingsFromGeometry()                    */
162/************************************************************************/
163
164static void GDALCollectRingsFromGeometry(
165    OGRGeometry *poShape, 
166    std::vector<double> &aPointX, std::vector<double> &aPointY, std::vector<double> &aPointVariant, 
167    std::vector<int> &aPartSize )
168
169{
170    if( poShape == NULL )
171        return;
172
173    OGRwkbGeometryType eFlatType = wkbFlatten(poShape->getGeometryType());
174    int i;
175
176    if ( eFlatType == wkbPoint )
177    {
178        OGRPoint    *poPoint = (OGRPoint *) poShape;
179        int nNewCount = aPointX.size() + 1;
180
181        aPointX.reserve( nNewCount );
182        aPointY.reserve( nNewCount );
183                aPointVariant.reserve( nNewCount );
184        aPointX.push_back( poPoint->getX() );
185        aPointY.push_back( poPoint->getY() );
186                aPointVariant.push_back( poPoint->getZ() ); // this could be getM() if available
187        aPartSize.push_back( 1 );
188    }
189    else if ( eFlatType == wkbLineString )
190    {
191        OGRLineString   *poLine = (OGRLineString *) poShape;
192        int nCount = poLine->getNumPoints();
193        int nNewCount = aPointX.size() + nCount;
194
195        aPointX.reserve( nNewCount );
196        aPointY.reserve( nNewCount );
197                aPointVariant.reserve( nNewCount );
198        for ( i = nCount - 1; i >= 0; i-- )
199        {
200            aPointX.push_back( poLine->getX(i) );
201            aPointY.push_back( poLine->getY(i) );
202                        aPointVariant.push_back( poLine->getZ(i) ); // this could be getM() if available
203        }
204        aPartSize.push_back( nCount );
205    }
206    else if ( EQUAL(poShape->getGeometryName(),"LINEARRING") )
207    {
208        OGRLinearRing *poRing = (OGRLinearRing *) poShape;
209        int nCount = poRing->getNumPoints();
210        int nNewCount = aPointX.size() + nCount;
211
212        aPointX.reserve( nNewCount );
213        aPointY.reserve( nNewCount );
214                aPointVariant.reserve( nNewCount );
215        for ( i = nCount - 1; i >= 0; i-- )
216        {
217            aPointX.push_back( poRing->getX(i) );
218            aPointY.push_back( poRing->getY(i) );
219                        aPointVariant.push_back( poRing->getY(i) ); // this could be getM() if available
220        }
221        aPartSize.push_back( nCount );
222    }
223    else if( eFlatType == wkbPolygon )
224    {
225        OGRPolygon *poPolygon = (OGRPolygon *) poShape;
226       
227        GDALCollectRingsFromGeometry( poPolygon->getExteriorRing(), 
228                                      aPointX, aPointY, aPointVariant, aPartSize );
229
230        for( i = 0; i < poPolygon->getNumInteriorRings(); i++ )
231            GDALCollectRingsFromGeometry( poPolygon->getInteriorRing(i), 
232                                          aPointX, aPointY, aPointVariant, aPartSize );
233    }
234   
235    else if( eFlatType == wkbMultiPoint
236             || eFlatType == wkbMultiLineString
237             || eFlatType == wkbMultiPolygon
238             || eFlatType == wkbGeometryCollection )
239    {
240        OGRGeometryCollection *poGC = (OGRGeometryCollection *) poShape;
241
242        for( i = 0; i < poGC->getNumGeometries(); i++ )
243            GDALCollectRingsFromGeometry( poGC->getGeometryRef(i),
244                                          aPointX, aPointY, aPointVariant, aPartSize );
245    }
246    else
247    {
248        CPLDebug( "GDAL", "Rasterizer ignoring non-polygonal geometry." );
249    }
250}
251
252/************************************************************************/
253/*                       gv_rasterize_one_shape()                       */
254/************************************************************************/
255static void 
256gv_rasterize_new_one_shape( unsigned char *pabyChunkBuf, int nYOff, int nYSize,
257                            int nBands, GDALDataType eType, GDALDataset *poDS,
258                            OGRGeometry *poShape, double *padfBurnValue, int bBurnValueSource,
259                            GDALTransformerFunc pfnTransformer, 
260                            void *pTransformArg )
261
262{
263    GDALRasterizeInfo sInfo;
264
265    sInfo.nXSize = poDS->GetRasterXSize();
266    sInfo.nYSize = nYSize;
267    sInfo.nBands = nBands;
268    sInfo.pabyChunkBuf = pabyChunkBuf;
269    sInfo.eType = eType;
270    sInfo.padfBurnValue = padfBurnValue;
271        sInfo.bBurnValueSource = bBurnValueSource;
272
273/* -------------------------------------------------------------------- */
274/*      Transform polygon geometries into a set of rings and a part     */
275/*      size list.                                                      */
276/* -------------------------------------------------------------------- */
277    std::vector<double> aPointX;
278    std::vector<double> aPointY;
279        std::vector<double> aPointVariant;
280    std::vector<int> aPartSize;
281
282    GDALCollectRingsFromGeometry( poShape, aPointX, aPointY, aPointVariant, aPartSize );
283
284/* -------------------------------------------------------------------- */
285/*      Transform points if needed.                                     */
286/* -------------------------------------------------------------------- */
287    if( pfnTransformer != NULL )
288    {
289        int *panSuccess = (int *) CPLCalloc(sizeof(int),aPointX.size());
290
291        // TODO: we need to add all appropriate error checking at some point.
292        pfnTransformer( pTransformArg, FALSE, aPointX.size(), 
293                        &(aPointX[0]), &(aPointY[0]), NULL, panSuccess );
294        CPLFree( panSuccess );
295    }
296
297/* -------------------------------------------------------------------- */
298/*      Shift to account for the buffer offset of this buffer.          */
299/* -------------------------------------------------------------------- */
300    unsigned int i;
301
302    for( i = 0; i < aPointY.size(); i++ )
303        aPointY[i] -= nYOff;
304
305/* -------------------------------------------------------------------- */
306/*      Perform the rasterization.                                      */
307/*      According to the C++ Standard/23.2.4, elements of a vector are  */
308/*      stored in continuous memory block.                              */
309/* -------------------------------------------------------------------- */
310
311    // TODO - mloskot: Check if vectors are empty, otherwise it may
312    // lead to undefined behavior by returning non-referencable pointer.
313    // if (!aPointX.empty())
314    //    /* fill polygon */
315    // else
316    //    /* How to report this problem? */   
317    switch ( wkbFlatten(poShape->getGeometryType()) )
318    {
319        case wkbPoint:
320        case wkbMultiPoint:
321            GDALdllImagePoint( sInfo.nXSize, nYSize, 
322                               aPartSize.size(), &(aPartSize[0]), 
323                               &(aPointX[0]), &(aPointY[0]), &(aPointVariant[0]),
324                               gvBurnPoint, &sInfo );
325            break;
326        case wkbLineString:
327        case wkbMultiLineString:
328            GDALdllImageLine( sInfo.nXSize, nYSize, 
329                              aPartSize.size(), &(aPartSize[0]), 
330                              &(aPointX[0]), &(aPointY[0]), &(aPointVariant[0]),
331                              gvBurnPoint, &sInfo );
332            break;
333        default:
334            GDALdllImageFilledPolygon( sInfo.nXSize, nYSize, 
335                                       aPartSize.size(), &(aPartSize[0]), 
336                                       &(aPointX[0]), &(aPointY[0]), &(aPointVariant[0]),
337                                       gvBurnScanline, &sInfo );
338            break;
339    }
340}
341
342/************************************************************************/
343/*                      GDALRasterizeGeometries()                       */
344/************************************************************************/
345
346/**
347 * Burn geometries into raster.
348 *
349 * Rasterize a list of geometric objects into a raster dataset.  The
350 * geometries are passed as an array of OGRGeometryH handlers. 
351 *
352 * If the geometries are in the georferenced coordinates of the raster
353 * dataset, then the pfnTransform may be passed in NULL and one will be
354 * derived internally from the geotransform of the dataset.  The transform
355 * needs to transform the geometry locations into pixel/line coordinates
356 * on the raster dataset.
357 *
358 * The output raster may be of any GDAL supported datatype, though currently
359 * internally the burning is done either as GDT_Byte or GDT_Float32.  This
360 * may be improved in the future.  An explicit list of burn values for
361 * each geometry for each band must be passed in.
362 *
363 * @param hDS output data, must be opened in update mode.
364 * @param nBandCount the number of bands to be updated.
365 * @param panBandList the list of bands to be updated.
366 * @param nGeomCount the number of geometries being passed in pahGeometries.
367 * @param pahGeometries the array of geometries to burn in.
368 * @param pfnTransformer transformation to apply to geometries to put into
369 * pixel/line coordinates on raster.  If NULL a geotransform based one will
370 * be created internally.
371 * @param pTransformerArg callback data for transformer.
372 * @param padfGeomBurnValue the array of values to burn into the raster. 
373 * There should be nBandCount values for each geometry.
374 * @param papszOption special options controlling rasterization.
375 * Currently one is defined, "BURN_VALUE_FROM". This option instructs
376 * GDALRasterizeGeometries(0 to derive the burn value from an alternate source.
377 * Currently the only recognized source is "Z", which is the interpolated Z
378 * coordinate of the source geometry, added to the explicit padfGeomBurnValue
379 * value. This is sensible for wkb[Multi]Point and wkb[Multi]LineString geometries
380 * but unimplemented for wkb[Multi]Polygon because interpolating Z within the
381 * filled polygon interior is extremely non-trivial.
382 * @param pfnProgress the progress function to report completion.
383 * @param pProgressArg callback data for progress function.
384 *
385 * @return CE_None on success or CE_Failure on error.
386 */
387
388CPLErr GDALRasterizeGeometries( GDALDatasetH hDS, 
389                                int nBandCount, int *panBandList,
390                                int nGeomCount, OGRGeometryH *pahGeometries,
391                                GDALTransformerFunc pfnTransformer, 
392                                void *pTransformArg, 
393                                double *padfGeomBurnValue,
394                                char **papszOptions,
395                                GDALProgressFunc pfnProgress, 
396                                void *pProgressArg )
397
398{
399    GDALDataType   eType;
400    int            nYChunkSize, nScanlineBytes;
401    unsigned char *pabyChunkBuf;
402    int            iY;
403    GDALDataset *poDS = (GDALDataset *) hDS;
404
405    if( pfnProgress == NULL )
406        pfnProgress = GDALDummyProgress;
407
408/* -------------------------------------------------------------------- */
409/*      Do some rudimentary arg checking.                               */
410/* -------------------------------------------------------------------- */
411    if( nBandCount == 0 || nGeomCount == 0 )
412        return CE_None;
413
414    // prototype band.
415    GDALRasterBand *poBand = poDS->GetRasterBand( panBandList[0] );
416
417        const char *pszOpt = CSLFetchNameValue( papszOptions, "BURN_VALUE_FROM" );
418        int bBurnValueSource = ( pszOpt && EQUAL(pszOpt,"Z")) ? 1 : 0; // Z = 1 other values undefined currently
419
420/* -------------------------------------------------------------------- */
421/*      If we have no transformer, assume the geometries are in file    */
422/*      georeferenced coordinates, and create a transformer to          */
423/*      convert that to pixel/line coordinates.                         */
424/*                                                                      */
425/*      We really just need to apply an affine transform, but for       */
426/*      simplicity we use the more general GenImgProjTransformer.       */
427/* -------------------------------------------------------------------- */
428    int bNeedToFreeTransformer = FALSE;
429
430    if( pfnTransformer == NULL )
431    {
432        bNeedToFreeTransformer = TRUE;
433
434        pTransformArg = 
435            GDALCreateGenImgProjTransformer( NULL, NULL, hDS, NULL, 
436                                             FALSE, 0.0, 0);
437        pfnTransformer = GDALGenImgProjTransform;
438    }
439
440/* -------------------------------------------------------------------- */
441/*      Establish a chunksize to operate on.  The larger the chunk      */
442/*      size the less times we need to make a pass through all the      */
443/*      shapes.                                                         */
444/* -------------------------------------------------------------------- */
445    if( poBand->GetRasterDataType() == GDT_Byte )
446        eType = GDT_Byte;
447    else
448        eType = GDT_Float32;
449
450    nScanlineBytes = nBandCount * poDS->GetRasterXSize()
451        * (GDALGetDataTypeSize(eType)/8);
452    nYChunkSize = 10000000 / nScanlineBytes;
453    if( nYChunkSize > poDS->GetRasterYSize() )
454        nYChunkSize = poDS->GetRasterYSize();
455
456    pabyChunkBuf = (unsigned char *) VSIMalloc(nYChunkSize * nScanlineBytes);
457    if( pabyChunkBuf == NULL )
458    {
459        CPLError( CE_Failure, CPLE_OutOfMemory, 
460                  "Unable to allocate rasterization buffer." );
461        return CE_Failure;
462    }
463
464/* ==================================================================== */
465/*      Loop over image in designated chunks.                           */
466/* ==================================================================== */
467    CPLErr  eErr = CE_None;
468
469    pfnProgress( 0.0, NULL, pProgressArg );
470
471    for( iY = 0; 
472         iY < poDS->GetRasterYSize() && eErr == CE_None; 
473         iY += nYChunkSize )
474    {
475        int     nThisYChunkSize;
476        int     iShape;
477
478        nThisYChunkSize = nYChunkSize;
479        if( nThisYChunkSize + iY > poDS->GetRasterYSize() )
480            nThisYChunkSize = poDS->GetRasterYSize() - iY;
481
482        eErr = 
483            poDS->RasterIO(GF_Read, 
484                           0, iY, poDS->GetRasterXSize(), nThisYChunkSize, 
485                           pabyChunkBuf,poDS->GetRasterXSize(),nThisYChunkSize,
486                           eType, nBandCount, panBandList,
487                           0, 0, 0 );
488        if( eErr != CE_None )
489            break;
490
491        for( iShape = 0; iShape < nGeomCount; iShape++ )
492        {
493            gv_rasterize_new_one_shape( pabyChunkBuf, iY, nThisYChunkSize,
494                                        nBandCount, eType, poDS,
495                                        (OGRGeometry *) pahGeometries[iShape],
496                                        padfGeomBurnValue + iShape*nBandCount, bBurnValueSource,
497                                        pfnTransformer, pTransformArg );
498        }
499
500        eErr = 
501            poDS->RasterIO( GF_Write, 0, iY,
502                            poDS->GetRasterXSize(), nThisYChunkSize, 
503                            pabyChunkBuf,
504                            poDS->GetRasterXSize(), nThisYChunkSize, 
505                            eType, nBandCount, panBandList, 0, 0, 0 );
506
507        if( !pfnProgress((iY+nThisYChunkSize)/((double)poDS->GetRasterYSize()),
508                         "", pProgressArg ) )
509        {
510            CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
511            eErr = CE_Failure;
512        }
513    }
514   
515/* -------------------------------------------------------------------- */
516/*      cleanup                                                         */
517/* -------------------------------------------------------------------- */
518    VSIFree( pabyChunkBuf );
519   
520    if( bNeedToFreeTransformer )
521        GDALDestroyTransformer( pTransformArg );
522
523    return eErr;
524}
525
526/************************************************************************/
527/*                        GDALRasterizeLayers()                         */
528/************************************************************************/
529
530/**
531 * Burn geometries from the specified layer into raster.
532 *
533 * Rasterize all the geometric objects from a list of layers into a raster
534 * dataset.  The layers are passed as an array of OGRLayerH handlers.
535 *
536 * If the geometries are in the georferenced coordinates of the raster
537 * dataset, then the pfnTransform may be passed in NULL and one will be
538 * derived internally from the geotransform of the dataset.  The transform
539 * needs to transform the geometry locations into pixel/line coordinates
540 * on the raster dataset.
541 *
542 * The output raster may be of any GDAL supported datatype, though currently
543 * internally the burning is done either as GDT_Byte or GDT_Float32.  This
544 * may be improved in the future.  An explicit list of burn values for
545 * each layer for each band must be passed in.
546 *
547 * @param hDS output data, must be opened in update mode.
548 * @param nBandCount the number of bands to be updated.
549 * @param panBandList the list of bands to be updated.
550 * @param nLayerCount the number of layers being passed in pahLayers array.
551 * @param pahLayers the array of layers to burn in.
552 * @param pfnTransformer transformation to apply to geometries to put into
553 * pixel/line coordinates on raster.  If NULL a geotransform based one will
554 * be created internally.
555 * @param pTransformerArg callback data for transformer.
556 * @param padfLayerBurnValues the array of values to burn into the raster. 
557 * There should be nBandCount values for each layer.
558 * @param papszOption special options controlling rasterization:
559 * <dl>
560 * <dt>"ATTRIBUTE":<dt> <dd>Identifies an attribute field on the features to be
561 * used for a burn in value. The value will be burned into all output
562 * bands. If specified, padfLayerBurnValues will not be used and can be a NULL
563 * pointer.</dd>
564 * <dt>"CHUNKYSIZE":<dt> <dd>The height in lines of the chunk to operate on.
565 * The larger the chunk size the less times we need to make a pass through all
566 * the shapes. If it is not set or set to zero the default chunk size will be
567 * used. Default size will be estimated based on the GDAL cache buffer size
568 * using formula: cache_size_bytes/scanline_size_bytes, so the chunk will
569 * not exceed the cache.</dd>
570 * </dl>
571 * @param pfnProgress the progress function to report completion.
572 * @param pProgressArg callback data for progress function.
573 *
574 * @return CE_None on success or CE_Failure on error.
575 */
576
577#ifdef OGR_ENABLED
578
579CPLErr GDALRasterizeLayers( GDALDatasetH hDS, 
580                            int nBandCount, int *panBandList,
581                            int nLayerCount, OGRLayerH *pahLayers,
582                            GDALTransformerFunc pfnTransformer, 
583                            void *pTransformArg, 
584                            double *padfLayerBurnValues,
585                            char **papszOptions,
586                            GDALProgressFunc pfnProgress, 
587                            void *pProgressArg )
588
589{
590    GDALDataType   eType;
591    unsigned char *pabyChunkBuf;
592    int            iY;
593    GDALDataset *poDS = (GDALDataset *) hDS;
594
595    if( pfnProgress == NULL )
596        pfnProgress = GDALDummyProgress;
597
598/* -------------------------------------------------------------------- */
599/*      Do some rudimentary arg checking.                               */
600/* -------------------------------------------------------------------- */
601    if( nBandCount == 0 || nLayerCount == 0 )
602        return CE_None;
603
604    // prototype band.
605    GDALRasterBand *poBand = poDS->GetRasterBand( panBandList[0] );
606
607        const char *pszOpt = CSLFetchNameValue( papszOptions, "BURN_VALUE_FROM" );
608        int bBurnValueSource = ( pszOpt && EQUAL(pszOpt,"Z")) ? 1 : 0; // Z = 1 other values undefined currently
609
610/* -------------------------------------------------------------------- */
611/*      Establish a chunksize to operate on.  The larger the chunk      */
612/*      size the less times we need to make a pass through all the      */
613/*      shapes.                                                         */
614/* -------------------------------------------------------------------- */
615    int         nYChunkSize, nScanlineBytes;
616    const char  *pszYChunkSize =
617        CSLFetchNameValue( papszOptions, "CHUNKYSIZE" );
618
619    if( poBand->GetRasterDataType() == GDT_Byte )
620        eType = GDT_Byte;
621    else
622        eType = GDT_Float32;
623
624    nScanlineBytes = nBandCount * poDS->GetRasterXSize()
625        * (GDALGetDataTypeSize(eType)/8);
626
627    if ( pszYChunkSize && (nYChunkSize = atoi(pszYChunkSize)) )
628        ;
629    else
630        nYChunkSize = GDALGetCacheMax() / nScanlineBytes;
631
632    if( nYChunkSize > poDS->GetRasterYSize() )
633        nYChunkSize = poDS->GetRasterYSize();
634
635    pabyChunkBuf = (unsigned char *) VSIMalloc(nYChunkSize * nScanlineBytes);
636    if( pabyChunkBuf == NULL )
637    {
638        CPLError( CE_Failure, CPLE_OutOfMemory, 
639                  "Unable to allocate rasterization buffer." );
640        return CE_Failure;
641    }
642
643/* ==================================================================== */
644/*      Read thespecified layers transfoming and rasterizing            */
645/*      geometries.                                                     */
646/* ==================================================================== */
647    CPLErr      eErr = CE_None;
648    int         iLayer;
649    const char  *pszBurnAttribute =
650        CSLFetchNameValue( papszOptions, "ATTRIBUTE" );
651
652    pfnProgress( 0.0, NULL, pProgressArg );
653
654    for( iLayer = 0; iLayer < nLayerCount; iLayer++ )
655    {
656        int         iBurnField = -1;
657        double      *padfBurnValues = NULL;
658        OGRLayer    *poLayer = (OGRLayer *) pahLayers[iLayer];
659
660        if ( !poLayer )
661        {
662            CPLError( CE_Warning, CPLE_AppDefined, 
663                      "Layer element number %d is NULL, skipping.\n", iLayer );
664            continue;
665        }
666
667/* -------------------------------------------------------------------- */
668/*      If the layer does not contain any features just skip it.        */
669/*      Do not force the feature count, so if driver doesn't know       */
670/*      exact number of features, go down the normal way.               */
671/* -------------------------------------------------------------------- */
672        if ( poLayer->GetFeatureCount(FALSE) == 0 )
673            continue;
674
675        if ( pszBurnAttribute )
676        {
677            iBurnField =
678                poLayer->GetLayerDefn()->GetFieldIndex( pszBurnAttribute );
679            if ( iBurnField == -1 )
680            {
681                CPLError( CE_Warning, CPLE_AppDefined, 
682                          "Failed to find field %s on layer %s, skipping.\n",
683                          pszBurnAttribute, 
684                          poLayer->GetLayerDefn()->GetName() );
685                continue;
686            }
687        }
688        else
689            padfBurnValues = padfLayerBurnValues + iLayer * nBandCount;
690
691/* -------------------------------------------------------------------- */
692/*      If we have no transformer, create the one from input file       */
693/*      projection. Note that each layer can be georefernced            */
694/*      separately.                                                     */
695/* -------------------------------------------------------------------- */
696        int bNeedToFreeTransformer = FALSE;
697
698        if( pfnTransformer == NULL )
699        {
700            char    *pszProjection;
701            bNeedToFreeTransformer = TRUE;
702
703            OGRSpatialReference *poSRS = poLayer->GetSpatialRef();
704            if ( !poSRS )
705            {
706                CPLError( CE_Warning, CPLE_AppDefined, 
707                          "Failed to fetch spatial reference on layer %s "
708                          "to build transformer, skipping.\n",
709                          poLayer->GetLayerDefn()->GetName() );
710                continue;
711            }
712
713            poSRS->exportToWkt( &pszProjection );
714
715            pTransformArg = 
716                GDALCreateGenImgProjTransformer( NULL, pszProjection,
717                                                 hDS, NULL, FALSE, 0.0, 0 );
718            pfnTransformer = GDALGenImgProjTransform;
719
720            CPLFree( pszProjection );
721        }
722
723        OGRFeature *poFeat;
724
725        poLayer->ResetReading();
726
727/* -------------------------------------------------------------------- */
728/*      Loop over image in designated chunks.                           */
729/* -------------------------------------------------------------------- */
730       for( iY = 0; 
731            iY < poDS->GetRasterYSize() && eErr == CE_None; 
732            iY += nYChunkSize )
733        {
734            int nThisYChunkSize;
735
736            nThisYChunkSize = nYChunkSize;
737            if( nThisYChunkSize + iY > poDS->GetRasterYSize() )
738                nThisYChunkSize = poDS->GetRasterYSize() - iY;
739
740            eErr = 
741                poDS->RasterIO( GF_Read, 0, iY,
742                                poDS->GetRasterXSize(), nThisYChunkSize, 
743                                pabyChunkBuf,
744                                poDS->GetRasterXSize(), nThisYChunkSize,
745                                eType, nBandCount, panBandList, 0, 0, 0 );
746            if( eErr != CE_None )
747                break;
748
749            while( (poFeat = poLayer->GetNextFeature()) != NULL )
750            {
751                OGRGeometry *poGeom = poFeat->GetGeometryRef();
752                double      dfBurnValue;
753
754                if ( pszBurnAttribute )
755                {
756                    dfBurnValue = poFeat->GetFieldAsDouble( iBurnField );
757                    padfBurnValues = &dfBurnValue;
758                }
759               
760                gv_rasterize_new_one_shape( pabyChunkBuf, iY, nThisYChunkSize,
761                                            nBandCount, eType, poDS, poGeom,
762                                            padfBurnValues, bBurnValueSource,
763                                            pfnTransformer, pTransformArg );
764
765                delete poFeat;
766            }
767
768            eErr = 
769                poDS->RasterIO( GF_Write, 0, iY,
770                                poDS->GetRasterXSize(), nThisYChunkSize, 
771                                pabyChunkBuf,
772                                poDS->GetRasterXSize(), nThisYChunkSize, 
773                                eType, nBandCount, panBandList, 0, 0, 0 );
774
775            poLayer->ResetReading();
776
777            if( !pfnProgress((iY+nThisYChunkSize)/((double)poDS->GetRasterYSize()),
778                             "", pProgressArg) )
779            {
780                CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
781                eErr = CE_Failure;
782            }
783        }
784
785        if ( bNeedToFreeTransformer )
786        {
787            GDALDestroyTransformer( pTransformArg );
788            pTransformArg = NULL;
789            pfnTransformer = NULL;
790        }
791    }
792   
793/* -------------------------------------------------------------------- */
794/*      cleanup                                                         */
795/* -------------------------------------------------------------------- */
796    VSIFree( pabyChunkBuf );
797   
798    return eErr;
799}
800
801#endif /* def OGR_ENABLED */