EVOLUTION-MANAGER

Edit File: gdalrasterize.cpp

/******************************************************************************
 *
 * Project:  GDAL
 * Purpose:  Vector rasterization.
 * Author:   Frank Warmerdam, warmerdam@pobox.com
 *
 ******************************************************************************
 * Copyright (c) 2005, Frank Warmerdam <warmerdam@pobox.com>
 * Copyright (c) 2008-2013, Even Rouault <even dot rouault at mines-paris dot org>
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 ****************************************************************************/

#include "cpl_port.h"
#include "gdal_alg.h"
#include "gdal_alg_priv.h"

#include <climits>
#include <cstddef>
#include <cstdlib>
#include <cstring>
#include <vector>
#include <algorithm>

#include "cpl_conv.h"
#include "cpl_error.h"
#include "cpl_progress.h"
#include "cpl_string.h"
#include "cpl_vsi.h"
#include "gdal.h"
#include "gdal_priv.h"
#include "ogr_api.h"
#include "ogr_core.h"
#include "ogr_feature.h"
#include "ogr_geometry.h"
#include "ogr_spatialref.h"
#include "ogrsf_frmts.h"

CPL_CVSID("$Id: gdalrasterize.cpp 9ff327806cd64df6d73a6c91f92d12ca0c5e07df 2018-04-07 20:25:06 +0200 Even Rouault $")

/************************************************************************/
/*                           gvBurnScanline()                           */
/************************************************************************/
static
void gvBurnScanline( void *pCBData, int nY, int nXStart, int nXEnd,
                     double dfVariant )

{
    if( nXStart > nXEnd )
        return;

GDALRasterizeInfo *psInfo = static_cast<GDALRasterizeInfo *>(pCBData);

CPLAssert( nY >= 0 && nY < psInfo->nYSize );
    CPLAssert( nXStart <= nXEnd );
    CPLAssert( nXStart < psInfo->nXSize );
    CPLAssert( nXEnd >= 0 );

if( nXStart < 0 )
        nXStart = 0;
    if( nXEnd >= psInfo->nXSize )
        nXEnd = psInfo->nXSize - 1;

if( psInfo->eType == GDT_Byte )
    {
        for( int iBand = 0; iBand < psInfo->nBands; iBand++ )
        {
            unsigned char nBurnValue = static_cast<unsigned char>
                ( psInfo->padfBurnValue[iBand] +
                  ( (psInfo->eBurnValueSource == GBV_UserBurnValue)?
                             0 : dfVariant ) );

unsigned char *pabyInsert =
                psInfo->pabyChunkBuf
                + iBand * psInfo->nXSize * psInfo->nYSize
                + nY * psInfo->nXSize + nXStart;

if( psInfo->eMergeAlg == GRMA_Add ) {
                int nPixels = nXEnd - nXStart + 1;
                while( nPixels-- > 0 )
                    *(pabyInsert++) += nBurnValue;
            } else {
                memset( pabyInsert, nBurnValue, nXEnd - nXStart + 1 );
            }
        }
    }
    else if( psInfo->eType == GDT_Float64 )
    {
        for( int iBand = 0; iBand < psInfo->nBands; iBand++ )
        {
            int nPixels = nXEnd - nXStart + 1;
            const double dfBurnValue =
                ( psInfo->padfBurnValue[iBand] +
                  ( (psInfo->eBurnValueSource == GBV_UserBurnValue)?
                             0 : dfVariant ) );

double *padfInsert =
                (reinterpret_cast<double *>(psInfo->pabyChunkBuf))
                + iBand * psInfo->nXSize * psInfo->nYSize
                + nY * psInfo->nXSize + nXStart;

if( psInfo->eMergeAlg == GRMA_Add ) {
                while( nPixels-- > 0 )
                    *(padfInsert++) += dfBurnValue;
            } else {
                while( nPixels-- > 0 )
                    *(padfInsert++) = dfBurnValue;
            }
        }
    }
    else {
        CPLAssert(false);
    }
}

/************************************************************************/
/*                            gvBurnPoint()                             */
/************************************************************************/
static
void gvBurnPoint( void *pCBData, int nY, int nX, double dfVariant )

{
    GDALRasterizeInfo *psInfo = static_cast<GDALRasterizeInfo *>(pCBData);

CPLAssert( nY >= 0 && nY < psInfo->nYSize );
    CPLAssert( nX >= 0 && nX < psInfo->nXSize );

if( psInfo->eType == GDT_Byte )
    {
        for( int iBand = 0; iBand < psInfo->nBands; iBand++ )
        {
            unsigned char *pbyInsert = psInfo->pabyChunkBuf
                                      + iBand * psInfo->nXSize * psInfo->nYSize
                                      + nY * psInfo->nXSize + nX;
            double dfVal;
            if( psInfo->eMergeAlg == GRMA_Add ) {
                dfVal = *pbyInsert + ( psInfo->padfBurnValue[iBand] +
                          ( (psInfo->eBurnValueSource == GBV_UserBurnValue)?
                             0 : dfVariant ) );
            } else {
                dfVal = psInfo->padfBurnValue[iBand] +
                          ( (psInfo->eBurnValueSource == GBV_UserBurnValue)?
                             0 : dfVariant );
            }
            if( dfVal > 255.0 )
                *pbyInsert = 255;
            else if( dfVal < 0.0 )
                *pbyInsert = 0;
            else
                *pbyInsert = static_cast<unsigned char>( dfVal );
        }
    }
    else if( psInfo->eType == GDT_Float64 )
    {
        for( int iBand = 0; iBand < psInfo->nBands; iBand++ )
        {
            double *pdfInsert = reinterpret_cast<double *>(psInfo->pabyChunkBuf)
                                + iBand * psInfo->nXSize * psInfo->nYSize
                                + nY * psInfo->nXSize + nX;

if( psInfo->eMergeAlg == GRMA_Add ) {
                *pdfInsert += ( psInfo->padfBurnValue[iBand] +
                         ( (psInfo->eBurnValueSource == GBV_UserBurnValue)?
                            0 : dfVariant ) );
            } else {
                *pdfInsert = ( psInfo->padfBurnValue[iBand] +
                         ( (psInfo->eBurnValueSource == GBV_UserBurnValue)?
                            0 : dfVariant ) );
            }
        }
    }
    else {
        CPLAssert(false);
    }
}

/************************************************************************/
/*                    GDALCollectRingsFromGeometry()                    */
/************************************************************************/

static void GDALCollectRingsFromGeometry(
    OGRGeometry *poShape,
    std::vector<double> &aPointX, std::vector<double> &aPointY,
    std::vector<double> &aPointVariant,
    std::vector<int> &aPartSize, GDALBurnValueSrc eBurnValueSrc)

{
    if( poShape == nullptr || poShape->IsEmpty() )
        return;

const OGRwkbGeometryType eFlatType = wkbFlatten(poShape->getGeometryType());

if( eFlatType == wkbPoint )
    {
        OGRPoint *poPoint = dynamic_cast<OGRPoint *>(poShape);
        CPLAssert(poPoint != nullptr);
        const size_t nNewCount = aPointX.size() + 1;

aPointX.reserve( nNewCount );
        aPointY.reserve( nNewCount );
        aPointX.push_back( poPoint->getX() );
        aPointY.push_back( poPoint->getY() );
        aPartSize.push_back( 1 );
        if( eBurnValueSrc != GBV_UserBurnValue )
        {
            // TODO(schwehr): Why not have the option for M r18164?
            // switch( eBurnValueSrc )
            // {
            // case GBV_Z:*/
            aPointVariant.reserve( nNewCount );
            aPointVariant.push_back( poPoint->getZ() );
            // break;
            // case GBV_M:
            //    aPointVariant.reserve( nNewCount );
            //    aPointVariant.push_back( poPoint->getM() );
        }
    }
    else if( eFlatType == wkbLineString )
    {
        OGRLineString *poLine = dynamic_cast<OGRLineString *>(poShape);
        CPLAssert(poLine != nullptr);
        const int nCount = poLine->getNumPoints();
        const size_t nNewCount = aPointX.size() + static_cast<size_t>(nCount);

aPointX.reserve( nNewCount );
        aPointY.reserve( nNewCount );
        if( eBurnValueSrc != GBV_UserBurnValue )
            aPointVariant.reserve( nNewCount );
        for( int i = nCount - 1; i >= 0; i-- )
        {
            aPointX.push_back( poLine->getX(i) );
            aPointY.push_back( poLine->getY(i) );
            if( eBurnValueSrc != GBV_UserBurnValue )
            {
                /*switch( eBurnValueSrc )
                {
                    case GBV_Z:*/
                        aPointVariant.push_back( poLine->getZ(i) );
                        /*break;
                    case GBV_M:
                        aPointVariant.push_back( poLine->getM(i) );
                }*/
            }
        }
        aPartSize.push_back( nCount );
    }
    else if( EQUAL(poShape->getGeometryName(), "LINEARRING") )
    {
        OGRLinearRing *poRing = dynamic_cast<OGRLinearRing *>(poShape);
        CPLAssert(poRing != nullptr);
        const int nCount = poRing->getNumPoints();
        const size_t nNewCount = aPointX.size() + static_cast<size_t>(nCount);

aPointX.reserve( nNewCount );
        aPointY.reserve( nNewCount );
        if( eBurnValueSrc != GBV_UserBurnValue )
            aPointVariant.reserve( nNewCount );
        int i = nCount - 1;  // Used after for.
        for( ; i >= 0; i-- )
        {
            aPointX.push_back( poRing->getX(i) );
            aPointY.push_back( poRing->getY(i) );
        }
        if( eBurnValueSrc != GBV_UserBurnValue )
        {
            /*switch( eBurnValueSrc )
            {
            case GBV_Z:*/
                aPointVariant.push_back( poRing->getZ(i) );
                /*break;
            case GBV_M:
                aPointVariant.push_back( poRing->getM(i) );
            }*/
        }
        aPartSize.push_back( nCount );
    }
    else if( eFlatType == wkbPolygon )
    {
        OGRPolygon *poPolygon = dynamic_cast<OGRPolygon *>(poShape);
        CPLAssert(poPolygon != nullptr);

GDALCollectRingsFromGeometry( poPolygon->getExteriorRing(),
                                      aPointX, aPointY, aPointVariant,
                                      aPartSize, eBurnValueSrc );

for( int i = 0; i < poPolygon->getNumInteriorRings(); i++ )
            GDALCollectRingsFromGeometry( poPolygon->getInteriorRing(i),
                                          aPointX, aPointY, aPointVariant,
                                          aPartSize, eBurnValueSrc );
    }
    else if( eFlatType == wkbMultiPoint
             || eFlatType == wkbMultiLineString
             || eFlatType == wkbMultiPolygon
             || eFlatType == wkbGeometryCollection )
    {
        OGRGeometryCollection *poGC = dynamic_cast<OGRGeometryCollection *>(poShape);
        CPLAssert(poGC != nullptr);

for( int i = 0; i < poGC->getNumGeometries(); i++ )
            GDALCollectRingsFromGeometry( poGC->getGeometryRef(i),
                                          aPointX, aPointY, aPointVariant,
                                          aPartSize, eBurnValueSrc );
    }
    else
    {
        CPLDebug( "GDAL", "Rasterizer ignoring non-polygonal geometry." );
    }
}

/************************************************************************/
/*                       gv_rasterize_one_shape()                       */
/************************************************************************/
static void
gv_rasterize_one_shape( unsigned char *pabyChunkBuf, int nXOff, int nYOff,
                        int nXSize, int nYSize,
                        int nBands, GDALDataType eType, int bAllTouched,
                        OGRGeometry *poShape, double *padfBurnValue,
                        GDALBurnValueSrc eBurnValueSrc,
                        GDALRasterMergeAlg eMergeAlg,
                        GDALTransformerFunc pfnTransformer,
                        void *pTransformArg )

{
    if( poShape == nullptr || poShape->IsEmpty() )
        return;

GDALRasterizeInfo sInfo;
    sInfo.nXSize = nXSize;
    sInfo.nYSize = nYSize;
    sInfo.nBands = nBands;
    sInfo.pabyChunkBuf = pabyChunkBuf;
    sInfo.eType = eType;
    sInfo.padfBurnValue = padfBurnValue;
    sInfo.eBurnValueSource = eBurnValueSrc;
    sInfo.eMergeAlg = eMergeAlg;

/* -------------------------------------------------------------------- */
/*      Transform polygon geometries into a set of rings and a part     */
/*      size list.                                                      */
/* -------------------------------------------------------------------- */
    std::vector<double> aPointX;
    std::vector<double> aPointY;
    std::vector<double> aPointVariant;
    std::vector<int> aPartSize;

GDALCollectRingsFromGeometry( poShape, aPointX, aPointY, aPointVariant,
                                  aPartSize, eBurnValueSrc );

/* -------------------------------------------------------------------- */
/*      Transform points if needed.                                     */
/* -------------------------------------------------------------------- */
    if( pfnTransformer != nullptr )
    {
        int *panSuccess =
            static_cast<int *>(CPLCalloc(sizeof(int), aPointX.size()));

// TODO: We need to add all appropriate error checking at some point.
        pfnTransformer( pTransformArg, FALSE, static_cast<int>(aPointX.size()),
                        &(aPointX[0]), &(aPointY[0]), nullptr, panSuccess );
        CPLFree( panSuccess );
    }

/* -------------------------------------------------------------------- */
/*      Shift to account for the buffer offset of this buffer.          */
/* -------------------------------------------------------------------- */
    for( unsigned int i = 0; i < aPointX.size(); i++ )
        aPointX[i] -= nXOff;
    for( unsigned int i = 0; i < aPointY.size(); i++ )
        aPointY[i] -= nYOff;

/* -------------------------------------------------------------------- */
/*      Perform the rasterization.                                      */
/*      According to the C++ Standard/23.2.4, elements of a vector are  */
/*      stored in continuous memory block.                              */
/* -------------------------------------------------------------------- */

// TODO - mloskot: Check if vectors are empty, otherwise it may
    // lead to undefined behavior by returning non-referencable pointer.
    // if( !aPointX.empty() )
    //    // Fill polygon.
    // else
    //    // How to report this problem?
    switch( wkbFlatten(poShape->getGeometryType()) )
    {
      case wkbPoint:
      case wkbMultiPoint:
        GDALdllImagePoint( sInfo.nXSize, nYSize,
                           static_cast<int>(aPartSize.size()), &(aPartSize[0]),
                           &(aPointX[0]), &(aPointY[0]),
                           (eBurnValueSrc == GBV_UserBurnValue)?
                           nullptr : &(aPointVariant[0]),
                           gvBurnPoint, &sInfo );
        break;
      case wkbLineString:
      case wkbMultiLineString:
      {
          if( bAllTouched )
              GDALdllImageLineAllTouched( sInfo.nXSize, nYSize,
                                          static_cast<int>(aPartSize.size()),
                                          &(aPartSize[0]),
                                          &(aPointX[0]), &(aPointY[0]),
                                          (eBurnValueSrc == GBV_UserBurnValue)?
                                          nullptr : &(aPointVariant[0]),
                                          gvBurnPoint, &sInfo );
          else
              GDALdllImageLine( sInfo.nXSize, nYSize,
                                static_cast<int>(aPartSize.size()),
                                &(aPartSize[0]),
                                &(aPointX[0]), &(aPointY[0]),
                                (eBurnValueSrc == GBV_UserBurnValue)?
                                nullptr : &(aPointVariant[0]),
                                gvBurnPoint, &sInfo );
      }
      break;

default:
      {
          GDALdllImageFilledPolygon(
              sInfo.nXSize, nYSize,
              static_cast<int>(aPartSize.size()), &(aPartSize[0]),
              &(aPointX[0]), &(aPointY[0]),
              (eBurnValueSrc == GBV_UserBurnValue)?
              nullptr : &(aPointVariant[0]),
              gvBurnScanline, &sInfo );
          if( bAllTouched )
          {
              // Reverting the variants to the first value because the
              // polygon is filled using the variant from the first point of
              // the first segment. Should be removed when the code to full
              // polygons more appropriately is added.
              if( eBurnValueSrc == GBV_UserBurnValue )
              {
                  GDALdllImageLineAllTouched(
                      sInfo.nXSize, nYSize,
                      static_cast<int>(aPartSize.size()), &(aPartSize[0]),
                      &(aPointX[0]), &(aPointY[0]),
                      nullptr,
                      gvBurnPoint, &sInfo );
              }
              else
              {
                  for( unsigned int i = 0, n = 0;
                       i < static_cast<unsigned int>(aPartSize.size());
                       i++ )
                  {
                      for( int j = 0; j < aPartSize[i]; j++ )
                          aPointVariant[n++] = aPointVariant[0];
                  }

GDALdllImageLineAllTouched(
                      sInfo.nXSize, nYSize,
                      static_cast<int>(aPartSize.size()), &(aPartSize[0]),
                      &(aPointX[0]), &(aPointY[0]),
                      &(aPointVariant[0]),
                      gvBurnPoint, &sInfo );
              }
          }
      }
      break;
    }
}

/************************************************************************/
/*                        GDALRasterizeOptions()                        */
/*                                                                      */
/*      Recognise a few rasterize options used by all three entry       */
/*      points.                                                         */
/************************************************************************/

static CPLErr GDALRasterizeOptions( char **papszOptions,
                                    int *pbAllTouched,
                                    GDALBurnValueSrc *peBurnValueSource,
                                    GDALRasterMergeAlg *peMergeAlg,
                                    GDALRasterizeOptim *peOptim)
{
    *pbAllTouched = CPLFetchBool( papszOptions, "ALL_TOUCHED", false );

const char *pszOpt = CSLFetchNameValue( papszOptions, "BURN_VALUE_FROM" );
    *peBurnValueSource = GBV_UserBurnValue;
    if( pszOpt )
    {
        if( EQUAL(pszOpt, "Z"))
        {
            *peBurnValueSource = GBV_Z;
        }
        // else if( EQUAL(pszOpt, "M"))
        //     eBurnValueSource = GBV_M;
        else
        {
            CPLError( CE_Failure, CPLE_AppDefined,
                      "Unrecognized value '%s' for BURN_VALUE_FROM.",
                      pszOpt );
            return CE_Failure;
        }
    }

/* -------------------------------------------------------------------- */
/*      MERGE_ALG=[REPLACE]/ADD                                         */
/* -------------------------------------------------------------------- */
    *peMergeAlg = GRMA_Replace;
    pszOpt = CSLFetchNameValue( papszOptions, "MERGE_ALG" );
    if( pszOpt )
    {
        if( EQUAL(pszOpt, "ADD"))
        {
            *peMergeAlg = GRMA_Add;
        }
        else if( EQUAL(pszOpt, "REPLACE"))
        {
            *peMergeAlg = GRMA_Replace;
        }
        else
        {
            CPLError( CE_Failure, CPLE_AppDefined,
                      "Unrecognized value '%s' for MERGE_ALG.",
                      pszOpt );
            return CE_Failure;
        }
    }

/* -------------------------------------------------------------------- */
/*      OPTIM=[AUTO]/RASTER/VECTOR                               */
/* -------------------------------------------------------------------- */
    *peOptim = GRO_Auto;
    pszOpt = CSLFetchNameValue( papszOptions, "OPTIM" );
    if( pszOpt )
    {
        if( EQUAL(pszOpt, "RASTER"))
        {
            *peOptim = GRO_Raster;
        }
        else if( EQUAL(pszOpt, "VECTOR"))
        {
            *peOptim = GRO_Vector;
        }
        else if( EQUAL(pszOpt, "AUTO"))
        {
            *peOptim = GRO_Auto;
        }
        else
        {
            CPLError( CE_Failure, CPLE_AppDefined,
                      "Unrecognized value '%s' for OPTIM.",
                      pszOpt );
            return CE_Failure;
        }
    }

return CE_None;
}

/************************************************************************/
/*                      GDALRasterizeGeometries()                       */
/************************************************************************/

/**
 * Burn geometries into raster.
 *
 * Rasterize a list of geometric objects into a raster dataset.  The
 * geometries are passed as an array of OGRGeometryH handlers.
 *
 * If the geometries are in the georeferenced coordinates of the raster
 * dataset, then the pfnTransform may be passed in NULL and one will be
 * derived internally from the geotransform of the dataset.  The transform
 * needs to transform the geometry locations into pixel/line coordinates
 * on the raster dataset.
 *
 * The output raster may be of any GDAL supported datatype, though currently
 * internally the burning is done either as GDT_Byte or GDT_Float32.  This
 * may be improved in the future.  An explicit list of burn values for
 * each geometry for each band must be passed in.
 *
 * The papszOption list of options currently only supports one option. The
 * "ALL_TOUCHED" option may be enabled by setting it to "TRUE".
 *
 * @param hDS output data, must be opened in update mode.
 * @param nBandCount the number of bands to be updated.
 * @param panBandList the list of bands to be updated.
 * @param nGeomCount the number of geometries being passed in pahGeometries.
 * @param pahGeometries the array of geometries to burn in.
 * @param pfnTransformer transformation to apply to geometries to put into
 * pixel/line coordinates on raster.  If NULL a geotransform based one will
 * be created internally.
 * @param pTransformArg callback data for transformer.
 * @param padfGeomBurnValue the array of values to burn into the raster.
 * There should be nBandCount values for each geometry.
 * @param papszOptions special options controlling rasterization
 * <ul>
 * <li>"ALL_TOUCHED": May be set to TRUE to set all pixels touched
 * by the line or polygons, not just those whose center is within the polygon
 * or that are selected by brezenhams line algorithm.  Defaults to FALSE.</li>
 * <li>"BURN_VALUE_FROM": May be set to "Z" to use the Z values of the
 * geometries. dfBurnValue is added to this before burning.
 * Defaults to GDALBurnValueSrc.GBV_UserBurnValue in which case just the
 * dfBurnValue is burned. This is implemented only for points and lines for
 * now. The M value may be supported in the future.</li>
 * <li>"MERGE_ALG": May be REPLACE (the default) or ADD.  REPLACE results in
 * overwriting of value, while ADD adds the new value to the existing raster,
 * suitable for heatmaps for instance.</li>
 * <li>"CHUNKYSIZE": The height in lines of the chunk to operate on.
 * The larger the chunk size the less times we need to make a pass through all
 * the shapes. If it is not set or set to zero the default chunk size will be
 * used. Default size will be estimated based on the GDAL cache buffer size
 * using formula: cache_size_bytes/scanline_size_bytes, so the chunk will
 * not exceed the cache. Not used in OPTIM=RASTER mode.</li>
 * </ul>
 * @param pfnProgress the progress function to report completion.
 * @param pProgressArg callback data for progress function.
 *
 * @return CE_None on success or CE_Failure on error.
 */

CPLErr GDALRasterizeGeometries( GDALDatasetH hDS,
                                int nBandCount, int *panBandList,
                                int nGeomCount, OGRGeometryH *pahGeometries,
                                GDALTransformerFunc pfnTransformer,
                                void *pTransformArg,
                                double *padfGeomBurnValue,
                                char **papszOptions,
                                GDALProgressFunc pfnProgress,
                                void *pProgressArg )

{
    VALIDATE_POINTER1( hDS, "GDALRasterizeGeometries", CE_Failure);

if( pfnProgress == nullptr )
        pfnProgress = GDALDummyProgress;

GDALDataset *poDS = reinterpret_cast<GDALDataset *>(hDS);

/* -------------------------------------------------------------------- */
/*      Do some rudimentary arg checking.                               */
/* -------------------------------------------------------------------- */
    if( nBandCount == 0 || nGeomCount == 0 )
    {
        pfnProgress(1.0, "", pProgressArg );
        return CE_None;
    }

// Prototype band.
    GDALRasterBand *poBand = poDS->GetRasterBand( panBandList[0] );
    if( poBand == nullptr )
        return CE_Failure;

/* -------------------------------------------------------------------- */
/*      Options                                                         */
/* -------------------------------------------------------------------- */
    int bAllTouched = FALSE;
    GDALBurnValueSrc eBurnValueSource = GBV_UserBurnValue;
    GDALRasterMergeAlg eMergeAlg = GRMA_Replace;
    GDALRasterizeOptim eOptim = GRO_Auto;
    if( GDALRasterizeOptions(papszOptions, &bAllTouched,
                             &eBurnValueSource, &eMergeAlg,
                             &eOptim) == CE_Failure )
    {
        return CE_Failure;
    }

/* -------------------------------------------------------------------- */
/*      If we have no transformer, assume the geometries are in file    */
/*      georeferenced coordinates, and create a transformer to          */
/*      convert that to pixel/line coordinates.                         */
/*                                                                      */
/*      We really just need to apply an affine transform, but for       */
/*      simplicity we use the more general GenImgProjTransformer.       */
/* -------------------------------------------------------------------- */
    bool bNeedToFreeTransformer = false;

if( pfnTransformer == nullptr )
    {
        bNeedToFreeTransformer = true;

char** papszTransformerOptions = nullptr;
        double adfGeoTransform[6] = { 0.0 };
        if( poDS->GetGeoTransform( adfGeoTransform ) != CE_None &&
            poDS->GetGCPCount() == 0 &&
            poDS->GetMetadata("RPC") == nullptr )
        {
            papszTransformerOptions = CSLSetNameValue(
                papszTransformerOptions, "DST_METHOD", "NO_GEOTRANSFORM");
        }

pTransformArg =
            GDALCreateGenImgProjTransformer2( nullptr, hDS,
                                                papszTransformerOptions );
        CSLDestroy( papszTransformerOptions );

pfnTransformer = GDALGenImgProjTransform;
        if( pTransformArg == nullptr )
        {
            return CE_Failure;
        }
    }

/* -------------------------------------------------------------------- */
/*      Choice of optimisation in auto mode. Use vector optim :         */
/*      1) if output is tiled                                           */
/*      2) if large number of features is present (>10000)              */
/*      3) if the nb of pixels > 50 * nb of features (not-too-small ft) */
/* -------------------------------------------------------------------- */
    int nXBlockSize, nYBlockSize;
    poBand->GetBlockSize(&nXBlockSize, &nYBlockSize);

if( eOptim == GRO_Auto )
    {
        eOptim = GRO_Raster;
        // TODO make more tests with various inputs/outputs to ajust the parameters
        if( nYBlockSize > 1 && nGeomCount > 10000 && (poBand->GetXSize() * static_cast<long long>(poBand->GetYSize()) / nGeomCount > 50) )
        {
            eOptim = GRO_Vector;
            CPLDebug("GDAL", "The vector optim has been chosen automatically");
        }
    }

/* -------------------------------------------------------------------- */
/*      The original algorithm                                          */
/*      Optimized for raster writing                                    */
/*      (optimal on a small number of large vectors)                    */
/* -------------------------------------------------------------------- */
    unsigned char *pabyChunkBuf;
    CPLErr eErr = CE_None;
    if( eOptim == GRO_Raster )
    {
/* -------------------------------------------------------------------- */
/*      Establish a chunksize to operate on.  The larger the chunk      */
/*      size the less times we need to make a pass through all the      */
/*      shapes.                                                         */
/* -------------------------------------------------------------------- */
        const GDALDataType eType =
            poBand->GetRasterDataType() == GDT_Byte ? GDT_Byte : GDT_Float64;

const int nScanlineBytes =
            nBandCount * poDS->GetRasterXSize() * GDALGetDataTypeSizeBytes(eType);

int nYChunkSize = 0;
        const char *pszYChunkSize = CSLFetchNameValue(papszOptions, "CHUNKYSIZE");
        if( pszYChunkSize == nullptr || ((nYChunkSize = atoi(pszYChunkSize))) == 0)
        {
            const GIntBig nYChunkSize64 = GDALGetCacheMax64() / nScanlineBytes;
            const int knIntMax = std::numeric_limits<int>::max();
            nYChunkSize = nYChunkSize64 > knIntMax ? knIntMax
                          : static_cast<int>(nYChunkSize64);
        }

if( nYChunkSize < 1 )
            nYChunkSize = 1;
        if( nYChunkSize > poDS->GetRasterYSize() )
            nYChunkSize = poDS->GetRasterYSize();

CPLDebug( "GDAL", "Rasterizer operating on %d swaths of %d scanlines.",
                  (poDS->GetRasterYSize() + nYChunkSize - 1) / nYChunkSize,
                  nYChunkSize );

pabyChunkBuf = static_cast<unsigned char *>(
            VSI_MALLOC2_VERBOSE(nYChunkSize, nScanlineBytes));
        if( pabyChunkBuf == nullptr )
        {
            if( bNeedToFreeTransformer )
                GDALDestroyTransformer( pTransformArg );
            return CE_Failure;
        }

/* ==================================================================== */
/*      Loop over image in designated chunks.                           */
/* ==================================================================== */
        pfnProgress( 0.0, nullptr, pProgressArg );

for( int iY = 0;
             iY < poDS->GetRasterYSize() && eErr == CE_None;
             iY += nYChunkSize )
        {
            int nThisYChunkSize = nYChunkSize;
            if( nThisYChunkSize + iY > poDS->GetRasterYSize() )
                nThisYChunkSize = poDS->GetRasterYSize() - iY;

eErr =
                poDS->RasterIO(GF_Read,
                               0, iY, poDS->GetRasterXSize(), nThisYChunkSize,
                               pabyChunkBuf,
                               poDS->GetRasterXSize(), nThisYChunkSize,
                               eType, nBandCount, panBandList,
                               0, 0, 0, nullptr);
            if( eErr != CE_None )
                break;

for( int iShape = 0; iShape < nGeomCount; iShape++ )
            {
                gv_rasterize_one_shape( pabyChunkBuf, 0, iY,
                                        poDS->GetRasterXSize(), nThisYChunkSize,
                                        nBandCount, eType, bAllTouched,
                                        reinterpret_cast<OGRGeometry *>(
                                                            pahGeometries[iShape]),
                                        padfGeomBurnValue + iShape*nBandCount,
                                        eBurnValueSource, eMergeAlg,
                                        pfnTransformer, pTransformArg );
            }

eErr =
                poDS->RasterIO( GF_Write, 0, iY,
                                poDS->GetRasterXSize(), nThisYChunkSize,
                                pabyChunkBuf,
                                poDS->GetRasterXSize(), nThisYChunkSize,
                                eType, nBandCount, panBandList, 0, 0, 0, nullptr);

if( !pfnProgress((iY + nThisYChunkSize) /
                             static_cast<double>(poDS->GetRasterYSize()),
                             "", pProgressArg ) )
            {
                CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
                eErr = CE_Failure;
            }
        }
    }
/* -------------------------------------------------------------------- */
/*      The new algorithm                                               */
/*      Optimized to minimize the vector computation                    */
/*      (optimal on a large number of vectors & tiled raster)           */
/* -------------------------------------------------------------------- */
    else
    {
/* -------------------------------------------------------------------- */
/*      Establish a chunksize to operate on.  Its size is defined by    */
/*      the block size of the output file.                              */
/* -------------------------------------------------------------------- */
        const int nXBlocks = (poBand->GetXSize() + nXBlockSize - 1) / nXBlockSize;
        const int nYBlocks = (poBand->GetYSize() + nYBlockSize - 1) / nYBlockSize;

const GDALDataType eType =
            poBand->GetRasterDataType() == GDT_Byte ? GDT_Byte : GDT_Float64;

const int nPixelSize = nBandCount * GDALGetDataTypeSizeBytes(eType);

// rem: optimized for square blocks
        const GIntBig nbMaxBlocks64 = GDALGetCacheMax64() / nPixelSize / nYBlockSize / nXBlockSize;
        const int knIntMax = std::numeric_limits<int>::max();
        const int nbMaxBlocks =
            nbMaxBlocks64 > knIntMax
            ? knIntMax : static_cast<int>(nbMaxBlocks64);
        const int nbBlocsX = std::max(1, std::min(static_cast<int>(sqrt(static_cast<double>(nbMaxBlocks))), nXBlocks));
        const int nbBlocsY = std::max(1, std::min(nbMaxBlocks / nbBlocsX, nYBlocks));

const int nScanblocks =
            nXBlockSize * nbBlocsX * nYBlockSize * nbBlocsY;

pabyChunkBuf = static_cast<unsigned char *>( VSI_MALLOC2_VERBOSE(nPixelSize, nScanblocks) );
        if( pabyChunkBuf == nullptr )
        {
            if( bNeedToFreeTransformer )
                GDALDestroyTransformer( pTransformArg );
            return CE_Failure;
        }

int * panSuccessTransform = static_cast<int *>(CPLCalloc(sizeof(int), 2));

/* -------------------------------------------------------------------- */
/*      loop over the vectorial geometries                              */
/* -------------------------------------------------------------------- */
        pfnProgress( 0.0, nullptr, pProgressArg );
        for( int iShape = 0; iShape < nGeomCount; iShape++ )
        {

OGRGeometry * poGeometry = reinterpret_cast<OGRGeometry *>(pahGeometries[iShape]);
            if ( poGeometry == nullptr || poGeometry->IsEmpty() )
              continue;
/* -------------------------------------------------------------------- */
/*      get the envelope of the geometry and transform it to pixels coo */
/* -------------------------------------------------------------------- */
            OGREnvelope psGeomEnvelope;
            poGeometry->getEnvelope(&psGeomEnvelope);
            if( pfnTransformer != nullptr )
            {
                double apCorners[4];
                apCorners[0] = psGeomEnvelope.MinX;
                apCorners[1] = psGeomEnvelope.MaxX;
                apCorners[2] = psGeomEnvelope.MinY;
                apCorners[3] = psGeomEnvelope.MaxY;
                // TODO: need to add all appropriate error checking
                pfnTransformer( pTransformArg, FALSE, 2, &(apCorners[0]),
                                &(apCorners[2]), nullptr, panSuccessTransform );
                psGeomEnvelope.MinX = std::min(apCorners[0], apCorners[1]);
                psGeomEnvelope.MaxX = std::max(apCorners[0], apCorners[1]);
                psGeomEnvelope.MinY = std::min(apCorners[2], apCorners[3]);
                psGeomEnvelope.MaxY = std::max(apCorners[2], apCorners[3]);
            }

int minBlockX = std::max(0, int(psGeomEnvelope.MinX) / nXBlockSize );
            int minBlockY = std::max(0, int(psGeomEnvelope.MinY) / nYBlockSize );
            int maxBlockX = std::min(nXBlocks-1, int(psGeomEnvelope.MaxX+1) / nXBlockSize );
            int maxBlockY = std::min(nYBlocks-1, int(psGeomEnvelope.MaxY+1) / nYBlockSize );

/* -------------------------------------------------------------------- */
/*      loop over the blocks concerned by the geometry                  */
/*      (by packs of nbBlocsX x nbBlocsY)                                 */
/* -------------------------------------------------------------------- */

for(int xB = minBlockX; xB <= maxBlockX; xB += nbBlocsX)
            {
                for(int yB = minBlockY; yB <= maxBlockY; yB += nbBlocsY)
                {

/* -------------------------------------------------------------------- */
/*      ensure to stay in the image                                     */
/* -------------------------------------------------------------------- */
                    int remSBX = std::min(maxBlockX - xB + 1, nbBlocsX);
                    int remSBY = std::min(maxBlockY - yB + 1, nbBlocsY);
                    int nThisXChunkSize = nXBlockSize * remSBX;
                    int nThisYChunkSize = nYBlockSize * remSBY;
                    if( xB * nXBlockSize + nThisXChunkSize > poDS->GetRasterXSize() )
                        nThisXChunkSize = poDS->GetRasterXSize() - xB * nXBlockSize;
                    if( yB * nYBlockSize + nThisYChunkSize > poDS->GetRasterYSize() )
                        nThisYChunkSize = poDS->GetRasterYSize() - yB * nYBlockSize;

/* -------------------------------------------------------------------- */
/*      read image / process buffer / write buffer                      */
/* -------------------------------------------------------------------- */
                    eErr = poDS->RasterIO(GF_Read, xB * nXBlockSize, yB * nYBlockSize, nThisXChunkSize, nThisYChunkSize,
                                       pabyChunkBuf, nThisXChunkSize, nThisYChunkSize, eType, nBandCount, panBandList,
                                       0, 0, 0, nullptr);
                    if( eErr != CE_None )
                        break;

gv_rasterize_one_shape( pabyChunkBuf, xB * nXBlockSize, yB * nYBlockSize,
                                            nThisXChunkSize, nThisYChunkSize,
                                            nBandCount, eType, bAllTouched,
                                            reinterpret_cast<OGRGeometry *>(pahGeometries[iShape]),
                                            padfGeomBurnValue + iShape*nBandCount,
                                            eBurnValueSource, eMergeAlg,
                                            pfnTransformer, pTransformArg );

eErr = poDS->RasterIO(GF_Write, xB * nXBlockSize, yB * nYBlockSize, nThisXChunkSize, nThisYChunkSize,
                                       pabyChunkBuf, nThisXChunkSize, nThisYChunkSize, eType, nBandCount, panBandList,
                                       0, 0, 0, nullptr);
                    if( eErr != CE_None )
                        break;
                }
            }

if( !pfnProgress(iShape / static_cast<double>(nGeomCount), "", pProgressArg ) )
            {
                CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
                eErr = CE_Failure;
            }

}

CPLFree( panSuccessTransform );

if( !pfnProgress(1., "", pProgressArg ) )
        {
            CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
            eErr = CE_Failure;
        }

}

/* -------------------------------------------------------------------- */
/*      cleanup                                                         */
/* -------------------------------------------------------------------- */
    VSIFree( pabyChunkBuf );

if( bNeedToFreeTransformer )
        GDALDestroyTransformer( pTransformArg );

return eErr;
}

/************************************************************************/
/*                        GDALRasterizeLayers()                         */
/************************************************************************/

/**
 * Burn geometries from the specified list of layers into raster.
 *
 * Rasterize all the geometric objects from a list of layers into a raster
 * dataset.  The layers are passed as an array of OGRLayerH handlers.
 *
 * If the geometries are in the georeferenced coordinates of the raster
 * dataset, then the pfnTransform may be passed in NULL and one will be
 * derived internally from the geotransform of the dataset.  The transform
 * needs to transform the geometry locations into pixel/line coordinates
 * on the raster dataset.
 *
 * The output raster may be of any GDAL supported datatype, though currently
 * internally the burning is done either as GDT_Byte or GDT_Float32.  This
 * may be improved in the future.  An explicit list of burn values for
 * each layer for each band must be passed in.
 *
 * @param hDS output data, must be opened in update mode.
 * @param nBandCount the number of bands to be updated.
 * @param panBandList the list of bands to be updated.
 * @param nLayerCount the number of layers being passed in pahLayers array.
 * @param pahLayers the array of layers to burn in.
 * @param pfnTransformer transformation to apply to geometries to put into
 * pixel/line coordinates on raster.  If NULL a geotransform based one will
 * be created internally.
 * @param pTransformArg callback data for transformer.
 * @param padfLayerBurnValues the array of values to burn into the raster.
 * There should be nBandCount values for each layer.
 * @param papszOptions special options controlling rasterization:
 * <ul>
 * <li>"ATTRIBUTE": Identifies an attribute field on the features to be
 * used for a burn in value. The value will be burned into all output
 * bands. If specified, padfLayerBurnValues will not be used and can be a NULL
 * pointer.</li>
 * <li>"CHUNKYSIZE": The height in lines of the chunk to operate on.
 * The larger the chunk size the less times we need to make a pass through all
 * the shapes. If it is not set or set to zero the default chunk size will be
 * used. Default size will be estimated based on the GDAL cache buffer size
 * using formula: cache_size_bytes/scanline_size_bytes, so the chunk will
 * not exceed the cache.</li>
 * <li>"ALL_TOUCHED": May be set to TRUE to set all pixels touched
 * by the line or polygons, not just those whose center is within the polygon
 * or that are selected by brezenhams line algorithm.  Defaults to FALSE.
 * <li>"BURN_VALUE_FROM": May be set to "Z" to use the Z values of the</li>
 * geometries. The value from padfLayerBurnValues or the attribute field value
 * is added to this before burning. In default case dfBurnValue is burned as it
 * is. This is implemented properly only for points and lines for now. Polygons
 * will be burned using the Z value from the first point. The M value may be
 * supported in the future.</li>
 * <li>"MERGE_ALG": May be REPLACE (the default) or ADD.  REPLACE results in
 * overwriting of value, while ADD adds the new value to the existing raster,
 * suitable for heatmaps for instance.</li>
 * </ul>
 * @param pfnProgress the progress function to report completion.
 * @param pProgressArg callback data for progress function.
 *
 * @return CE_None on success or CE_Failure on error.
 */

CPLErr GDALRasterizeLayers( GDALDatasetH hDS,
                            int nBandCount, int *panBandList,
                            int nLayerCount, OGRLayerH *pahLayers,
                            GDALTransformerFunc pfnTransformer,
                            void *pTransformArg,
                            double *padfLayerBurnValues,
                            char **papszOptions,
                            GDALProgressFunc pfnProgress,
                            void *pProgressArg )

{
    VALIDATE_POINTER1( hDS, "GDALRasterizeLayers", CE_Failure);

if( pfnProgress == nullptr )
        pfnProgress = GDALDummyProgress;

GDALDataset *poDS = reinterpret_cast<GDALDataset *>(hDS);

// Prototype band.
    GDALRasterBand *poBand = poDS->GetRasterBand( panBandList[0] );
    if( poBand == nullptr )
        return CE_Failure;

/* -------------------------------------------------------------------- */
/*      Establish a chunksize to operate on.  The larger the chunk      */
/*      size the less times we need to make a pass through all the      */
/*      shapes.                                                         */
/* -------------------------------------------------------------------- */
    const char  *pszYChunkSize =
        CSLFetchNameValue( papszOptions, "CHUNKYSIZE" );

const GDALDataType eType =
        poBand->GetRasterDataType() == GDT_Byte ? GDT_Byte : GDT_Float64;

const int nScanlineBytes =
        nBandCount * poDS->GetRasterXSize() * GDALGetDataTypeSizeBytes(eType);

int nYChunkSize = 0;
    if( !(pszYChunkSize && ((nYChunkSize = atoi(pszYChunkSize))) != 0) )
    {
        const GIntBig nYChunkSize64 = GDALGetCacheMax64() / nScanlineBytes;
        const int knIntMax = std::numeric_limits<int>::max();
        if( nYChunkSize64 > knIntMax )
            nYChunkSize = knIntMax;
        else
            nYChunkSize = static_cast<int>(nYChunkSize64);
    }

if( nYChunkSize < 1 )
        nYChunkSize = 1;
    if( nYChunkSize > poDS->GetRasterYSize() )
        nYChunkSize = poDS->GetRasterYSize();

CPLDebug( "GDAL", "Rasterizer operating on %d swaths of %d scanlines.",
              (poDS->GetRasterYSize() + nYChunkSize - 1) / nYChunkSize,
              nYChunkSize );
    unsigned char *pabyChunkBuf = static_cast<unsigned char *>(
        VSI_MALLOC2_VERBOSE(nYChunkSize, nScanlineBytes));
    if( pabyChunkBuf == nullptr )
    {
        return CE_Failure;
    }

/* -------------------------------------------------------------------- */
/*      Read the image once for all layers if user requested to render  */
/*      the whole raster in single chunk.                               */
/* -------------------------------------------------------------------- */
    if( nYChunkSize == poDS->GetRasterYSize() )
    {
        if( poDS->RasterIO( GF_Read, 0, 0, poDS->GetRasterXSize(),
                            nYChunkSize, pabyChunkBuf,
                            poDS->GetRasterXSize(), nYChunkSize,
                            eType, nBandCount, panBandList, 0, 0, 0, nullptr )
             != CE_None )
        {
            CPLFree( pabyChunkBuf );
            return CE_Failure;
        }
    }

/* ==================================================================== */
/*      Read the specified layers transforming and rasterizing          */
/*      geometries.                                                     */
/* ==================================================================== */
    CPLErr eErr = CE_None;
    const char *pszBurnAttribute = CSLFetchNameValue(papszOptions, "ATTRIBUTE");

pfnProgress( 0.0, nullptr, pProgressArg );

for( int iLayer = 0; iLayer < nLayerCount; iLayer++ )
    {
        OGRLayer *poLayer = reinterpret_cast<OGRLayer *>(pahLayers[iLayer]);

if( !poLayer )
        {
            CPLError( CE_Warning, CPLE_AppDefined,
                      "Layer element number %d is NULL, skipping.", iLayer );
            continue;
        }

/* -------------------------------------------------------------------- */
/*      If the layer does not contain any features just skip it.        */
/*      Do not force the feature count, so if driver doesn't know       */
/*      exact number of features, go down the normal way.               */
/* -------------------------------------------------------------------- */
        if( poLayer->GetFeatureCount(FALSE) == 0 )
            continue;

int iBurnField = -1;
        double *padfBurnValues = nullptr;

if( pszBurnAttribute )
        {
            iBurnField =
                poLayer->GetLayerDefn()->GetFieldIndex( pszBurnAttribute );
            if( iBurnField == -1 )
            {
                CPLError( CE_Warning, CPLE_AppDefined,
                          "Failed to find field %s on layer %s, skipping.",
                          pszBurnAttribute,
                          poLayer->GetLayerDefn()->GetName() );
                continue;
            }
        }
        else
        {
            padfBurnValues = padfLayerBurnValues + iLayer * nBandCount;
        }

/* -------------------------------------------------------------------- */
/*      If we have no transformer, create the one from input file       */
/*      projection. Note that each layer can be georefernced            */
/*      separately.                                                     */
/* -------------------------------------------------------------------- */
        bool bNeedToFreeTransformer = false;

if( pfnTransformer == nullptr )
        {
            char *pszProjection = nullptr;
            bNeedToFreeTransformer = true;

OGRSpatialReference *poSRS = poLayer->GetSpatialRef();
            if( !poSRS )
            {
                CPLError( CE_Warning, CPLE_AppDefined,
                          "Failed to fetch spatial reference on layer %s "
                          "to build transformer, assuming matching coordinate "
                          "systems.",
                          poLayer->GetLayerDefn()->GetName() );
            }
            else
            {
                poSRS->exportToWkt( &pszProjection );
            }

char** papszTransformerOptions = nullptr;
            if( pszProjection != nullptr )
                papszTransformerOptions = CSLSetNameValue(
                        papszTransformerOptions, "SRC_SRS", pszProjection );
            double adfGeoTransform[6] = {};
            if( poDS->GetGeoTransform( adfGeoTransform ) != CE_None &&
                poDS->GetGCPCount() == 0 &&
                poDS->GetMetadata("RPC") == nullptr )
            {
                papszTransformerOptions = CSLSetNameValue(
                    papszTransformerOptions, "DST_METHOD", "NO_GEOTRANSFORM");
            }

pTransformArg =
                GDALCreateGenImgProjTransformer2( nullptr, hDS,
                                                  papszTransformerOptions );
            pfnTransformer = GDALGenImgProjTransform;

CPLFree( pszProjection );
            CSLDestroy( papszTransformerOptions );
            if( pTransformArg == nullptr )
            {
                CPLFree( pabyChunkBuf );
                return CE_Failure;
            }
        }

poLayer->ResetReading();

/* -------------------------------------------------------------------- */
/*      Loop over image in designated chunks.                           */
/* -------------------------------------------------------------------- */

double *padfAttrValues = static_cast<double *>(
            VSI_MALLOC_VERBOSE(sizeof(double) * nBandCount));
        if( padfAttrValues == nullptr )
            eErr = CE_Failure;

// Only re-read image if not a single chunk is being rendered.
            if( nYChunkSize < poDS->GetRasterYSize() )
            {
                eErr =
                    poDS->RasterIO( GF_Read, 0, iY,
                                    poDS->GetRasterXSize(), nThisYChunkSize,
                                    pabyChunkBuf,
                                    poDS->GetRasterXSize(), nThisYChunkSize,
                                    eType, nBandCount, panBandList,
                                    0, 0, 0, nullptr );
                if( eErr != CE_None )
                    break;
            }

OGRFeature *poFeat = nullptr;
            while( (poFeat = poLayer->GetNextFeature()) != nullptr )
            {
                OGRGeometry *poGeom = poFeat->GetGeometryRef();

if( pszBurnAttribute )
                {
                    const double dfAttrValue =
                        poFeat->GetFieldAsDouble( iBurnField );
                    for( int iBand = 0 ; iBand < nBandCount ; iBand++)
                        padfAttrValues[iBand] = dfAttrValue;

padfBurnValues = padfAttrValues;
                }

gv_rasterize_one_shape( pabyChunkBuf, 0, iY,
                                        poDS->GetRasterXSize(),
                                        nThisYChunkSize,
                                        nBandCount, eType, bAllTouched, poGeom,
                                        padfBurnValues, eBurnValueSource,
                                        eMergeAlg,
                                        pfnTransformer, pTransformArg );

delete poFeat;
            }

// Only write image if not a single chunk is being rendered.
            if( nYChunkSize < poDS->GetRasterYSize() )
            {
                eErr =
                    poDS->RasterIO( GF_Write, 0, iY,
                                    poDS->GetRasterXSize(), nThisYChunkSize,
                                    pabyChunkBuf,
                                    poDS->GetRasterXSize(), nThisYChunkSize,
                                    eType, nBandCount, panBandList,
                                    0, 0, 0, nullptr );
            }

poLayer->ResetReading();

VSIFree( padfAttrValues );

if( bNeedToFreeTransformer )
        {
            GDALDestroyTransformer( pTransformArg );
            pTransformArg = nullptr;
            pfnTransformer = nullptr;
        }
    }

/* -------------------------------------------------------------------- */
/*      Write out the image once for all layers if user requested       */
/*      to render the whole raster in single chunk.                     */
/* -------------------------------------------------------------------- */
    if( eErr == CE_None && nYChunkSize == poDS->GetRasterYSize() )
    {
        eErr = poDS->RasterIO( GF_Write, 0, 0,
                                poDS->GetRasterXSize(), nYChunkSize,
                                pabyChunkBuf,
                                poDS->GetRasterXSize(), nYChunkSize,
                                eType, nBandCount, panBandList, 0, 0, 0, nullptr );
    }

return eErr;
}

/************************************************************************/
/*                        GDALRasterizeLayersBuf()                      */
/************************************************************************/

/**
 * Burn geometries from the specified list of layer into raster.
 *
 * Rasterize all the geometric objects from a list of layers into supplied
 * raster buffer.  The layers are passed as an array of OGRLayerH handlers.
 *
 * If the geometries are in the georeferenced coordinates of the raster
 * dataset, then the pfnTransform may be passed in NULL and one will be
 * derived internally from the geotransform of the dataset.  The transform
 * needs to transform the geometry locations into pixel/line coordinates
 * of the target raster.
 *
 * The output raster may be of any GDAL supported datatype, though currently
 * internally the burning is done either as GDT_Byte or GDT_Float32.  This
 * may be improved in the future.
 *
 * @param pData pointer to the output data array.
 *
 * @param nBufXSize width of the output data array in pixels.
 *
 * @param nBufYSize height of the output data array in pixels.
 *
 * @param eBufType data type of the output data array.
 *
 * @param nPixelSpace The byte offset from the start of one pixel value in
 * pData to the start of the next pixel value within a scanline.  If defaulted
 * (0) the size of the datatype eBufType is used.
 *
 * @param nLineSpace The byte offset from the start of one scanline in
 * pData to the start of the next.  If defaulted the size of the datatype
 * eBufType * nBufXSize is used.
 *
 * @param nLayerCount the number of layers being passed in pahLayers array.
 *
 * @param pahLayers the array of layers to burn in.
 *
 * @param pszDstProjection WKT defining the coordinate system of the target
 * raster.
 *
 * @param padfDstGeoTransform geotransformation matrix of the target raster.
 *
 * @param pfnTransformer transformation to apply to geometries to put into
 * pixel/line coordinates on raster.  If NULL a geotransform based one will
 * be created internally.
 *
 * @param pTransformArg callback data for transformer.
 *
 * @param dfBurnValue the value to burn into the raster.
 *
 * @param papszOptions special options controlling rasterization:
 * <ul>
 * <li>"ATTRIBUTE": Identifies an attribute field on the features to be
 * used for a burn in value. The value will be burned into all output
 * bands. If specified, padfLayerBurnValues will not be used and can be a NULL
 * pointer.</li>
 * <li>"ALL_TOUCHED": May be set to TRUE to set all pixels touched
 * by the line or polygons, not just those whose center is within the polygon
 * or that are selected by brezenhams line algorithm.  Defaults to FALSE.</li>
 * <li>"BURN_VALUE_FROM": May be set to "Z" to use
 * the Z values of the geometries. dfBurnValue or the attribute field value is
 * added to this before burning. In default case dfBurnValue is burned as it
 * is. This is implemented properly only for points and lines for now. Polygons
 * will be burned using the Z value from the first point. The M value may
 * be supported in the future.</li>
 * <li>"MERGE_ALG": May be REPLACE (the default) or ADD.  REPLACE
 * results in overwriting of value, while ADD adds the new value to the
 * existing raster, suitable for heatmaps for instance.</li>
 * </ul>
 *
 * @param pfnProgress the progress function to report completion.
 *
 * @param pProgressArg callback data for progress function.
 *
 *
 * @return CE_None on success or CE_Failure on error.
 */

CPLErr GDALRasterizeLayersBuf( void *pData, int nBufXSize, int nBufYSize,
                               GDALDataType eBufType,
                               int nPixelSpace, int nLineSpace,
                               int nLayerCount, OGRLayerH *pahLayers,
                               const char *pszDstProjection,
                               double *padfDstGeoTransform,
                               GDALTransformerFunc pfnTransformer,
                               void *pTransformArg, double dfBurnValue,
                               char **papszOptions,
                               GDALProgressFunc pfnProgress,
                               void *pProgressArg )

{
/* -------------------------------------------------------------------- */
/*      If pixel and line spaceing are defaulted assign reasonable      */
/*      value assuming a packed buffer.                                 */
/* -------------------------------------------------------------------- */
    if( nPixelSpace != 0 )
    {
        nPixelSpace = GDALGetDataTypeSizeBytes( eBufType );
    }
    if( nPixelSpace != GDALGetDataTypeSizeBytes( eBufType ) )
    {
        CPLError(CE_Failure, CPLE_NotSupported,
                 "GDALRasterizeLayersBuf(): unsupported value of nPixelSpace");
        return CE_Failure;
    }

if( nLineSpace == 0 )
    {
        nLineSpace = nPixelSpace * nBufXSize;
    }
    if( nLineSpace != nPixelSpace * nBufXSize )
    {
        CPLError(CE_Failure, CPLE_NotSupported,
                 "GDALRasterizeLayersBuf(): unsupported value of nLineSpace");
        return CE_Failure;
    }

if( pfnProgress == nullptr )
        pfnProgress = GDALDummyProgress;

/* ==================================================================== */
/*      Read the specified layers transforming and rasterizing          */
/*      geometries.                                                     */
/* ==================================================================== */
    CPLErr eErr = CE_None;
    const char  *pszBurnAttribute =
        CSLFetchNameValue( papszOptions, "ATTRIBUTE" );

pfnProgress( 0.0, nullptr, pProgressArg );

for( int iLayer = 0; iLayer < nLayerCount; iLayer++ )
    {
        OGRLayer *poLayer = reinterpret_cast<OGRLayer *>(pahLayers[iLayer]);

if( !poLayer )
        {
            CPLError( CE_Warning, CPLE_AppDefined,
                      "Layer element number %d is NULL, skipping.", iLayer );
            continue;
        }

int iBurnField = -1;
        if( pszBurnAttribute )
        {
            iBurnField =
                poLayer->GetLayerDefn()->GetFieldIndex( pszBurnAttribute );
            if( iBurnField == -1 )
            {
                CPLError( CE_Warning, CPLE_AppDefined,
                          "Failed to find field %s on layer %s, skipping.",
                          pszBurnAttribute,
                          poLayer->GetLayerDefn()->GetName() );
                continue;
            }
        }

if( pfnTransformer == nullptr )
        {
            char *pszProjection = nullptr;
            bNeedToFreeTransformer = true;

pTransformArg =
                GDALCreateGenImgProjTransformer3( pszProjection, nullptr,
                                                  pszDstProjection,
                                                  padfDstGeoTransform );
            pfnTransformer = GDALGenImgProjTransform;

CPLFree( pszProjection );
        }

poLayer->ResetReading();

{
            OGRFeature *poFeat = nullptr;
            while( (poFeat = poLayer->GetNextFeature()) != nullptr )
            {
                OGRGeometry *poGeom = poFeat->GetGeometryRef();

if( pszBurnAttribute )
                    dfBurnValue = poFeat->GetFieldAsDouble( iBurnField );

gv_rasterize_one_shape( static_cast<unsigned char *>(pData), 0, 0,
                                        nBufXSize, nBufYSize,
                                        1, eBufType, bAllTouched, poGeom,
                                        &dfBurnValue, eBurnValueSource,
                                        eMergeAlg,
                                        pfnTransformer, pTransformArg );

delete poFeat;
            }
        }

poLayer->ResetReading();

if( !pfnProgress(1, "", pProgressArg) )
        {
            CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
            eErr = CE_Failure;
        }

if( bNeedToFreeTransformer )
        {
            GDALDestroyTransformer( pTransformArg );
            pTransformArg = nullptr;
            pfnTransformer = nullptr;
        }
    }

return eErr;
}