/****************************************************************************** * $Id: gdalrasterize.cpp 15665 2008-10-31 18:18:35Z warmerdam $ * * Project: GDAL * Purpose: Vector rasterization. * Author: Frank Warmerdam, warmerdam@pobox.com * ****************************************************************************** * Copyright (c) 2005, Frank Warmerdam * * 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 #include "gdal_alg.h" #include "gdal_alg_priv.h" #include "gdal_priv.h" #include "ogr_api.h" #include "ogr_geometry.h" #include "ogr_spatialref.h" #ifdef OGR_ENABLED #include "ogrsf_frmts.h" #endif typedef struct { unsigned char * pabyChunkBuf; int nXSize; int nYSize; int nBands; GDALDataType eType; double *padfBurnValue; int bBurnValueSource; // 0 indicates burn value comes from padfBurnValue, 1 indicates value from Variant (Z, or possibly M in the future) } GDALRasterizeInfo; /************************************************************************/ /* gvBurnScanline() */ /************************************************************************/ void gvBurnScanline( void *pCBData, int nY, int nXStart, int nXEnd ) { GDALRasterizeInfo *psInfo = (GDALRasterizeInfo *) pCBData; int iBand; if( nXStart > nXEnd ) return; 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( iBand = 0; iBand < psInfo->nBands; iBand++ ) { unsigned char *pabyInsert; unsigned char nBurnValue = (unsigned char) psInfo->padfBurnValue[iBand]; pabyInsert = psInfo->pabyChunkBuf + iBand * psInfo->nXSize * psInfo->nYSize + nY * psInfo->nXSize + nXStart; memset( pabyInsert, nBurnValue, nXEnd - nXStart + 1 ); } } else { for( iBand = 0; iBand < psInfo->nBands; iBand++ ) { int nPixels = nXEnd - nXStart + 1; float *pafInsert; float fBurnValue = (float) psInfo->padfBurnValue[iBand]; pafInsert = ((float *) psInfo->pabyChunkBuf) + iBand * psInfo->nXSize * psInfo->nYSize + nY * psInfo->nXSize + nXStart; while( nPixels-- > 0 ) *(pafInsert++) = fBurnValue; } } } /************************************************************************/ /* gvBurnPoint() */ /************************************************************************/ void gvBurnPoint( void *pCBData, int nY, int nX, double dfVariant ) { GDALRasterizeInfo *psInfo = (GDALRasterizeInfo *) pCBData; int iBand; CPLAssert( nY >= 0 && nY < psInfo->nYSize ); CPLAssert( nX >= 0 && nX < psInfo->nXSize ); if( psInfo->eType == GDT_Byte ) { for( iBand = 0; iBand < psInfo->nBands; iBand++ ) { unsigned char *pbyInsert = psInfo->pabyChunkBuf + iBand * psInfo->nXSize * psInfo->nYSize + nY * psInfo->nXSize + nX; if(psInfo->bBurnValueSource) // should we include the Variant as our Burn Value source? { *pbyInsert = (unsigned char) dfVariant + (unsigned char)psInfo->padfBurnValue[iBand]; } // if else { *pbyInsert = (unsigned char) psInfo->padfBurnValue[iBand]; } // else } } else { for( iBand = 0; iBand < psInfo->nBands; iBand++ ) { float *pfInsert = ((float *) psInfo->pabyChunkBuf) + iBand * psInfo->nXSize * psInfo->nYSize + nY * psInfo->nXSize + nX; if(psInfo->bBurnValueSource) // should we use the Variant as our Burn Value source? { *pfInsert = (float) dfVariant + (float) psInfo->padfBurnValue[iBand]; } // if else { *pfInsert = (float) psInfo->padfBurnValue[iBand]; } // else } } } /************************************************************************/ /* GDALCollectRingsFromGeometry() */ /************************************************************************/ static void GDALCollectRingsFromGeometry( OGRGeometry *poShape, std::vector &aPointX, std::vector &aPointY, std::vector &aPointVariant, std::vector &aPartSize ) { if( poShape == NULL ) return; OGRwkbGeometryType eFlatType = wkbFlatten(poShape->getGeometryType()); int i; if ( eFlatType == wkbPoint ) { OGRPoint *poPoint = (OGRPoint *) poShape; int nNewCount = aPointX.size() + 1; aPointX.reserve( nNewCount ); aPointY.reserve( nNewCount ); aPointVariant.reserve( nNewCount ); aPointX.push_back( poPoint->getX() ); aPointY.push_back( poPoint->getY() ); aPointVariant.push_back( poPoint->getZ() ); // this could be getM() if available aPartSize.push_back( 1 ); } else if ( eFlatType == wkbLineString ) { OGRLineString *poLine = (OGRLineString *) poShape; int nCount = poLine->getNumPoints(); int nNewCount = aPointX.size() + nCount; aPointX.reserve( nNewCount ); aPointY.reserve( nNewCount ); aPointVariant.reserve( nNewCount ); for ( i = nCount - 1; i >= 0; i-- ) { aPointX.push_back( poLine->getX(i) ); aPointY.push_back( poLine->getY(i) ); aPointVariant.push_back( poLine->getZ(i) ); // this could be getM() if available } aPartSize.push_back( nCount ); } else if ( EQUAL(poShape->getGeometryName(),"LINEARRING") ) { OGRLinearRing *poRing = (OGRLinearRing *) poShape; int nCount = poRing->getNumPoints(); int nNewCount = aPointX.size() + nCount; aPointX.reserve( nNewCount ); aPointY.reserve( nNewCount ); aPointVariant.reserve( nNewCount ); for ( i = nCount - 1; i >= 0; i-- ) { aPointX.push_back( poRing->getX(i) ); aPointY.push_back( poRing->getY(i) ); aPointVariant.push_back( poRing->getY(i) ); // this could be getM() if available } aPartSize.push_back( nCount ); } else if( eFlatType == wkbPolygon ) { OGRPolygon *poPolygon = (OGRPolygon *) poShape; GDALCollectRingsFromGeometry( poPolygon->getExteriorRing(), aPointX, aPointY, aPointVariant, aPartSize ); for( i = 0; i < poPolygon->getNumInteriorRings(); i++ ) GDALCollectRingsFromGeometry( poPolygon->getInteriorRing(i), aPointX, aPointY, aPointVariant, aPartSize ); } else if( eFlatType == wkbMultiPoint || eFlatType == wkbMultiLineString || eFlatType == wkbMultiPolygon || eFlatType == wkbGeometryCollection ) { OGRGeometryCollection *poGC = (OGRGeometryCollection *) poShape; for( i = 0; i < poGC->getNumGeometries(); i++ ) GDALCollectRingsFromGeometry( poGC->getGeometryRef(i), aPointX, aPointY, aPointVariant, aPartSize ); } else { CPLDebug( "GDAL", "Rasterizer ignoring non-polygonal geometry." ); } } /************************************************************************/ /* gv_rasterize_one_shape() */ /************************************************************************/ static void gv_rasterize_new_one_shape( unsigned char *pabyChunkBuf, int nYOff, int nYSize, int nBands, GDALDataType eType, GDALDataset *poDS, OGRGeometry *poShape, double *padfBurnValue, int bBurnValueSource, GDALTransformerFunc pfnTransformer, void *pTransformArg ) { GDALRasterizeInfo sInfo; sInfo.nXSize = poDS->GetRasterXSize(); sInfo.nYSize = nYSize; sInfo.nBands = nBands; sInfo.pabyChunkBuf = pabyChunkBuf; sInfo.eType = eType; sInfo.padfBurnValue = padfBurnValue; sInfo.bBurnValueSource = bBurnValueSource; /* -------------------------------------------------------------------- */ /* Transform polygon geometries into a set of rings and a part */ /* size list. */ /* -------------------------------------------------------------------- */ std::vector aPointX; std::vector aPointY; std::vector aPointVariant; std::vector aPartSize; GDALCollectRingsFromGeometry( poShape, aPointX, aPointY, aPointVariant, aPartSize ); /* -------------------------------------------------------------------- */ /* Transform points if needed. */ /* -------------------------------------------------------------------- */ if( pfnTransformer != NULL ) { int *panSuccess = (int *) CPLCalloc(sizeof(int),aPointX.size()); // TODO: we need to add all appropriate error checking at some point. pfnTransformer( pTransformArg, FALSE, aPointX.size(), &(aPointX[0]), &(aPointY[0]), NULL, panSuccess ); CPLFree( panSuccess ); } /* -------------------------------------------------------------------- */ /* Shift to account for the buffer offset of this buffer. */ /* -------------------------------------------------------------------- */ unsigned int i; for( 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, aPartSize.size(), &(aPartSize[0]), &(aPointX[0]), &(aPointY[0]), &(aPointVariant[0]), gvBurnPoint, &sInfo ); break; case wkbLineString: case wkbMultiLineString: GDALdllImageLine( sInfo.nXSize, nYSize, aPartSize.size(), &(aPartSize[0]), &(aPointX[0]), &(aPointY[0]), &(aPointVariant[0]), gvBurnPoint, &sInfo ); break; default: GDALdllImageFilledPolygon( sInfo.nXSize, nYSize, aPartSize.size(), &(aPartSize[0]), &(aPointX[0]), &(aPointY[0]), &(aPointVariant[0]), gvBurnScanline, &sInfo ); break; } } /************************************************************************/ /* 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 georferenced 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. * * @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 pTransformerArg callback data for transformer. * @param padfGeomBurnValue the array of values to burn into the raster. * There should be nBandCount values for each geometry. * @param papszOption special options controlling rasterization. * Currently one is defined, "BURN_VALUE_FROM". This option instructs * GDALRasterizeGeometries(0 to derive the burn value from an alternate source. * Currently the only recognized source is "Z", which is the interpolated Z * coordinate of the source geometry, added to the explicit padfGeomBurnValue * value. This is sensible for wkb[Multi]Point and wkb[Multi]LineString geometries * but unimplemented for wkb[Multi]Polygon because interpolating Z within the * filled polygon interior is extremely non-trivial. * @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 ) { GDALDataType eType; int nYChunkSize, nScanlineBytes; unsigned char *pabyChunkBuf; int iY; GDALDataset *poDS = (GDALDataset *) hDS; if( pfnProgress == NULL ) pfnProgress = GDALDummyProgress; /* -------------------------------------------------------------------- */ /* Do some rudimentary arg checking. */ /* -------------------------------------------------------------------- */ if( nBandCount == 0 || nGeomCount == 0 ) return CE_None; // prototype band. GDALRasterBand *poBand = poDS->GetRasterBand( panBandList[0] ); const char *pszOpt = CSLFetchNameValue( papszOptions, "BURN_VALUE_FROM" ); int bBurnValueSource = ( pszOpt && EQUAL(pszOpt,"Z")) ? 1 : 0; // Z = 1 other values undefined currently /* -------------------------------------------------------------------- */ /* 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. */ /* -------------------------------------------------------------------- */ int bNeedToFreeTransformer = FALSE; if( pfnTransformer == NULL ) { bNeedToFreeTransformer = TRUE; pTransformArg = GDALCreateGenImgProjTransformer( NULL, NULL, hDS, NULL, FALSE, 0.0, 0); pfnTransformer = GDALGenImgProjTransform; } /* -------------------------------------------------------------------- */ /* Establish a chunksize to operate on. The larger the chunk */ /* size the less times we need to make a pass through all the */ /* shapes. */ /* -------------------------------------------------------------------- */ if( poBand->GetRasterDataType() == GDT_Byte ) eType = GDT_Byte; else eType = GDT_Float32; nScanlineBytes = nBandCount * poDS->GetRasterXSize() * (GDALGetDataTypeSize(eType)/8); nYChunkSize = 10000000 / nScanlineBytes; if( nYChunkSize > poDS->GetRasterYSize() ) nYChunkSize = poDS->GetRasterYSize(); pabyChunkBuf = (unsigned char *) VSIMalloc(nYChunkSize * nScanlineBytes); if( pabyChunkBuf == NULL ) { CPLError( CE_Failure, CPLE_OutOfMemory, "Unable to allocate rasterization buffer." ); return CE_Failure; } /* ==================================================================== */ /* Loop over image in designated chunks. */ /* ==================================================================== */ CPLErr eErr = CE_None; pfnProgress( 0.0, NULL, pProgressArg ); for( iY = 0; iY < poDS->GetRasterYSize() && eErr == CE_None; iY += nYChunkSize ) { int nThisYChunkSize; int iShape; 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 ); if( eErr != CE_None ) break; for( iShape = 0; iShape < nGeomCount; iShape++ ) { gv_rasterize_new_one_shape( pabyChunkBuf, iY, nThisYChunkSize, nBandCount, eType, poDS, (OGRGeometry *) pahGeometries[iShape], padfGeomBurnValue + iShape*nBandCount, bBurnValueSource, pfnTransformer, pTransformArg ); } eErr = poDS->RasterIO( GF_Write, 0, iY, poDS->GetRasterXSize(), nThisYChunkSize, pabyChunkBuf, poDS->GetRasterXSize(), nThisYChunkSize, eType, nBandCount, panBandList, 0, 0, 0 ); if( !pfnProgress((iY+nThisYChunkSize)/((double)poDS->GetRasterYSize()), "", 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 layer 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 georferenced 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 pTransformerArg callback data for transformer. * @param padfLayerBurnValues the array of values to burn into the raster. * There should be nBandCount values for each layer. * @param papszOption special options controlling rasterization: *
*
"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.
*
"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.
*
* @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. */ #ifdef OGR_ENABLED CPLErr GDALRasterizeLayers( GDALDatasetH hDS, int nBandCount, int *panBandList, int nLayerCount, OGRLayerH *pahLayers, GDALTransformerFunc pfnTransformer, void *pTransformArg, double *padfLayerBurnValues, char **papszOptions, GDALProgressFunc pfnProgress, void *pProgressArg ) { GDALDataType eType; unsigned char *pabyChunkBuf; int iY; GDALDataset *poDS = (GDALDataset *) hDS; if( pfnProgress == NULL ) pfnProgress = GDALDummyProgress; /* -------------------------------------------------------------------- */ /* Do some rudimentary arg checking. */ /* -------------------------------------------------------------------- */ if( nBandCount == 0 || nLayerCount == 0 ) return CE_None; // prototype band. GDALRasterBand *poBand = poDS->GetRasterBand( panBandList[0] ); const char *pszOpt = CSLFetchNameValue( papszOptions, "BURN_VALUE_FROM" ); int bBurnValueSource = ( pszOpt && EQUAL(pszOpt,"Z")) ? 1 : 0; // Z = 1 other values undefined currently /* -------------------------------------------------------------------- */ /* Establish a chunksize to operate on. The larger the chunk */ /* size the less times we need to make a pass through all the */ /* shapes. */ /* -------------------------------------------------------------------- */ int nYChunkSize, nScanlineBytes; const char *pszYChunkSize = CSLFetchNameValue( papszOptions, "CHUNKYSIZE" ); if( poBand->GetRasterDataType() == GDT_Byte ) eType = GDT_Byte; else eType = GDT_Float32; nScanlineBytes = nBandCount * poDS->GetRasterXSize() * (GDALGetDataTypeSize(eType)/8); if ( pszYChunkSize && (nYChunkSize = atoi(pszYChunkSize)) ) ; else nYChunkSize = GDALGetCacheMax() / nScanlineBytes; if( nYChunkSize > poDS->GetRasterYSize() ) nYChunkSize = poDS->GetRasterYSize(); pabyChunkBuf = (unsigned char *) VSIMalloc(nYChunkSize * nScanlineBytes); if( pabyChunkBuf == NULL ) { CPLError( CE_Failure, CPLE_OutOfMemory, "Unable to allocate rasterization buffer." ); return CE_Failure; } /* ==================================================================== */ /* Read thespecified layers transfoming and rasterizing */ /* geometries. */ /* ==================================================================== */ CPLErr eErr = CE_None; int iLayer; const char *pszBurnAttribute = CSLFetchNameValue( papszOptions, "ATTRIBUTE" ); pfnProgress( 0.0, NULL, pProgressArg ); for( iLayer = 0; iLayer < nLayerCount; iLayer++ ) { int iBurnField = -1; double *padfBurnValues = NULL; OGRLayer *poLayer = (OGRLayer *) pahLayers[iLayer]; if ( !poLayer ) { CPLError( CE_Warning, CPLE_AppDefined, "Layer element number %d is NULL, skipping.\n", 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; if ( pszBurnAttribute ) { iBurnField = poLayer->GetLayerDefn()->GetFieldIndex( pszBurnAttribute ); if ( iBurnField == -1 ) { CPLError( CE_Warning, CPLE_AppDefined, "Failed to find field %s on layer %s, skipping.\n", 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. */ /* -------------------------------------------------------------------- */ int bNeedToFreeTransformer = FALSE; if( pfnTransformer == NULL ) { char *pszProjection; bNeedToFreeTransformer = TRUE; OGRSpatialReference *poSRS = poLayer->GetSpatialRef(); if ( !poSRS ) { CPLError( CE_Warning, CPLE_AppDefined, "Failed to fetch spatial reference on layer %s " "to build transformer, skipping.\n", poLayer->GetLayerDefn()->GetName() ); continue; } poSRS->exportToWkt( &pszProjection ); pTransformArg = GDALCreateGenImgProjTransformer( NULL, pszProjection, hDS, NULL, FALSE, 0.0, 0 ); pfnTransformer = GDALGenImgProjTransform; CPLFree( pszProjection ); } OGRFeature *poFeat; poLayer->ResetReading(); /* -------------------------------------------------------------------- */ /* Loop over image in designated chunks. */ /* -------------------------------------------------------------------- */ for( iY = 0; iY < poDS->GetRasterYSize() && eErr == CE_None; iY += nYChunkSize ) { int nThisYChunkSize; 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 ); if( eErr != CE_None ) break; while( (poFeat = poLayer->GetNextFeature()) != NULL ) { OGRGeometry *poGeom = poFeat->GetGeometryRef(); double dfBurnValue; if ( pszBurnAttribute ) { dfBurnValue = poFeat->GetFieldAsDouble( iBurnField ); padfBurnValues = &dfBurnValue; } gv_rasterize_new_one_shape( pabyChunkBuf, iY, nThisYChunkSize, nBandCount, eType, poDS, poGeom, padfBurnValues, bBurnValueSource, pfnTransformer, pTransformArg ); delete poFeat; } eErr = poDS->RasterIO( GF_Write, 0, iY, poDS->GetRasterXSize(), nThisYChunkSize, pabyChunkBuf, poDS->GetRasterXSize(), nThisYChunkSize, eType, nBandCount, panBandList, 0, 0, 0 ); poLayer->ResetReading(); if( !pfnProgress((iY+nThisYChunkSize)/((double)poDS->GetRasterYSize()), "", pProgressArg) ) { CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" ); eErr = CE_Failure; } } if ( bNeedToFreeTransformer ) { GDALDestroyTransformer( pTransformArg ); pTransformArg = NULL; pfnTransformer = NULL; } } /* -------------------------------------------------------------------- */ /* cleanup */ /* -------------------------------------------------------------------- */ VSIFree( pabyChunkBuf ); return eErr; } #endif /* def OGR_ENABLED */