[gdal-dev] GDAL : Access window out of range in RasterIO and blank buffer in warp for same block extents
asmita
asmitawankhede at gmail.com
Mon Feb 13 16:08:12 PST 2017
I was using the WarpingReader but since it somehow broke the warp part in the
upgrade from GDAL 1.1 to 2.1.2 , I tried to investigate the read with
simple read, however I got this error :
Access window out of range in RasterIO(). Requested
(0,24104) of size 430x450 on raster of 11783x14544.
However the Original problem is the warp operation itself is not outputting
any data in the buffer &rawReadBuf[0]
( so I tried the SimpleWarpReader, but that is not how it suppose to be).
When I call the Warpregiontobuffer , I see the same values, but it does not
error out however it does not even populated the output buffer:
*** srcExtents.width 430
*** srcExtents.height 450
*** srcextents.beginX() 0
*** srcextents.beginY() 24104
In warpOperation::WarpRegionToBuffer however, I dont see that error ( mainly
I believe because its defaulting to ReadBlock ?) Anyways, the point is that
how should I proceed with to get successful read of the intermediate file
created from two or more tiffs to create a mosaic?
What Warp to use is the decision made in ImageReader:
class ImageReader
{
const ImageInfo &imageInfo;
xyzDeleteGuard<xyzGDALReader> reader;
public:
explicit ImageReader(const ImageInfo &imageInfo_) : imageInfo(imageInfo_)
{
// Refer GDALCreateApproxTransformer documentation in gdal_warp.c
// The maximum cartesian error in the "output" space that is to be
accepted
// in the linear approximation.
// The value 0.125 is as default in xyzGDALWarpingReader.
const double max_error = 0.125;
if (imageInfo.srcDS.needReproject()) {
notify(NFY_NOTICE, "Reprojecting source using %s",
GDALResampleAlgorithmName(imageInfo.resampleAlg));
const double* noData = NULL;
if (imageInfo.rasterType == xyzRasterProduct::Heightmap)
noData = &imageInfo.noDataValue;
reader = TransferOwnership(
new xyzGDALWarpingReader(
imageInfo.srcDS, imageInfo.resampleAlg,
imageInfo.numbands, max_error,
imageInfo.projectionType, noData));
} else if (imageInfo.srcDS.needSnapUp()) {
// we snap with the warper, it knows to simplify the
// reproject logic when it only needs a snapup
notify(NFY_NOTICE, "Snapping source using %s",
GDALResampleAlgorithmName(imageInfo.resampleAlg));
reader = TransferOwnership(
new xyzGDALWarpingReader(imageInfo.srcDS,
imageInfo.resampleAlg,
imageInfo.numbands,
max_error,
imageInfo.projectionType,
imageInfo.rasterType == xyzRasterProduct::Heightmap ?
&imageInfo.noDataValue
: NULL));
} else {
// image that exactly matches our projection and pixel size
// requirements
reader = TransferOwnership(
new xyzGDALSimpleReader(imageInfo.srcDS, imageInfo.numbands));
}
}
// Read tile from source.
// Do all coordinate calculations, filling, warping, palette lookup,
// datatype conversion, band replication, etc.
// Will throw exception if unable to complete read
template <class TileType>
void ReadTile(uint32 row, uint32 col, TileType &tile);
};
WarpingReader and SimpleReader Implementation:
#include "xyzGDALReader.h"
#include <xyzException.h>
//
****************************************************************************
// *** xyzGDALReader
//
****************************************************************************
xyzGDALReader::xyzGDALReader(const xyzGDALDataset &srcDS_, uint numbands)
: srcDS(srcDS_),
numBands(numbands),
gdalDatatype(srcDS->GetRasterBand(1)->GetRasterDataType()),
topToBottom(srcDS.normalizedGeoExtents().topToBottom()),
paletteSize(0),
palette(0)
{
uint gdalRasterCount = srcDS->GetRasterCount();
if (numBands == 0) {
throw xyzException
(xyz::tr("Internal Error: xyzGDALReader has no bands"));
}
if (numBands > gdalRasterCount) {
throw xyzException
(xyz::tr("Internal Error: xyzGDALReader w/ too many bands"));
}
if (!StorageEnumFromGDT(gdalDatatype, storage)) {
throw xyzException(xyz::tr("Unsupported GDAL datatype '%1'")
.arg(GDALGetDataTypeName(gdalDatatype)));
}
if (!srcDS.geoExtents().leftToRight()) {
throw xyzException(xyz::tr("Right to left images not supported"));
}
// populate list of gdal bands
gdalBandIndexes.resize(numBands);
const uint num_orig_bands = srcDS->GetRasterCount();
if (numBands == 1 || numBands == num_orig_bands) {
// In both these cases none of the bands need to be dropped.
for (uint i = 0; i < numBands; ++i) {
gdalBandIndexes[i] = i+1;
}
} else if (num_orig_bands == 4 && numBands == 3) { // 4 band images
// Ignore if only one alpha band or 4th band if none is alpha band
for (uint i = 0, j = 0, num_alpha = 0; i < num_orig_bands;) {
++i;
if (srcDS->GetRasterBand(i)->GetColorInterpretation() ==
GCI_AlphaBand) {
if (++num_alpha > 1) {
throw xyzException(xyz::tr("The image has more than one alpha
band."));
}
} else if (j < numBands) {
gdalBandIndexes[j] = i;
++j;
}
}
} else {
throw xyzException(xyz::tr("Band list has wrong number of bands"));
}
// get palette information
GDALRasterBand *band = srcDS->GetRasterBand(1);
if (band->GetColorInterpretation() == GCI_PaletteIndex) {
if (numBands > 1) {
throw xyzException
(xyz::tr("Unsupported: palette image with multiple bands"));
}
GDALColorTable *ctbl = band->GetColorTable();
if (!ctbl) {
throw xyzException(xyz::tr("Unable to get color table"));
}
paletteSize = (uint)ctbl->GetColorEntryCount();
palette = ctbl->GetColorEntry(0);
if (!paletteSize || !palette) {
throw xyzException(xyz::tr("Unable to get color table"));
}
}
}
//
****************************************************************************
// *** xyzGDALSimpleReader
//
****************************************************************************
void
xyzGDALSimpleReader::FetchPixels(const xyzExtents<uint32> &srcExtents)
{
if (srcDS->RasterIO(GF_Read,
srcExtents.beginX(), srcExtents.beginY(),
srcExtents.width(), srcExtents.height(),
&rawReadBuf[0],
srcExtents.width(), srcExtents.height(),
gdalDatatype,
numBands, &gdalBandIndexes[0],
0, 0, 0) != CE_None) {
throw xyzException(xyz::tr("Failed to read source data"));
}
}
//
****************************************************************************
// *** xyzGDALWarpingReader
//
****************************************************************************
TransformArgGuard&
TransformArgGuard::operator=(void *a)
{
if (arg && (arg != a)) {
destroyFunc(arg);
}
arg = a;
return *this;
}
TransformArgGuard::~TransformArgGuard(void)
{
if (arg) {
destroyFunc(arg);
}
}
int
SnapTransformFunc( void *pTransformArg, int bDstToSrc,
int nPointCount,
double *padfX, double *padfY, double *padfZ,
int *panSuccess )
{
if (!bDstToSrc) {
notify(NFY_FATAL, "Internal error: Wrong Snap direction");
}
SnapArg *snapArg = (SnapArg*)pTransformArg;
for (int i = 0; i < nPointCount; ++i) {
padfX[i] *= snapArg->xScale;
padfY[i] *= snapArg->yScale;
panSuccess[i] = 1;
}
return TRUE;
}
xyzGDALWarpingReader::xyzGDALWarpingReader(const xyzGDALDataset &srcDS,
GDALResampleAlg resampleAlg,
uint numbands,
double transformErrorThreshold,
xyzTilespace::ProjectionType
projType,
const double* noData)
: xyzGDALReader(srcDS, numbands),
snapArg(srcDS.geoExtents().geoTransform(),
srcDS.alignedGeoExtents().geoTransform()),
genImgProjArg(GDALDestroyGenImgProjTransformer),
approxArg(GDALDestroyApproxTransformer)
{
std::string srcSRS(srcDS->GetProjectionRef());
std::string dstSRS;
if (projType == xyzTilespace::FLAT_PROJECTION) {
dstSRS = GetWGS84ProjString();
} else if (projType == xyzTilespace::MERCATOR_PROJECTION) {
dstSRS = GetMercatorProjString();
} else {
throw xyzException(xyz::tr("Unknown projection type"));
}
// ***** Make a dummy output dataset *****
// We need a dataset in order to get the transformer. The
// transformer will only use the geoTransform out of the dataset,
// but you can't pass the geoTransform directly.
GDALDriverH memDriver = GDALGetDriverByName("MEM");
if (!memDriver) {
throw xyzException(xyz::tr("Unable to get GDAL MEM driver"));
}
xyzDeleteGuard<GDALDataset>
dstDS(TransferOwnership(
(GDALDataset*)
GDALCreate(memDriver, "unused filename",
1 /* numCols */,
1 /* numRows */,
numBands,
gdalDatatype,
0 /* unused create options */)));
if (!dstDS) {
throw xyzException(xyz::tr("Unable to create MEM dataset"));
}
assert(CE_None == dstDS->SetProjection(dstSRS.c_str()));
dstDS->SetGeoTransform(const_cast<double*>
(srcDS.alignedGeoExtents().geoTransform()));
GDALTransformerFunc transformFunc = 0;
void* transformArg = 0;
if (srcDS.needSnapUp()) {
notify(NFY_NOTICE, "Using snap transform ...");
transformArg = &snapArg;
transformFunc = SnapTransformFunc;
} else if (srcDS.needReproject()) {
// ***** make a transform w/ the dummy dstDS *****
genImgProjArg = GDALCreateGenImgProjTransformer
(srcDS, srcSRS.c_str(),
(GDALDataset*)dstDS, dstSRS.c_str(),
FALSE, /* use src GCP if necessary */
0.0, /* GCP error threshold */
0 /* maximum order for GCP polynomials */);
if (!genImgProjArg) {
throw xyzException(xyz::tr
("Unable to determine reproject transform"));
}
transformFunc = GDALGenImgProjTransform;
transformArg = genImgProjArg;
// ***** see if an approximate transformation is good enough *****
if (transformErrorThreshold != 0.0) {
notify(NFY_NOTICE, "Approximating transform ...");
approxArg = GDALCreateApproxTransformer(GDALGenImgProjTransform,
genImgProjArg,
transformErrorThreshold);
transformFunc = GDALApproxTransform;
transformArg = approxArg;
} else {
notify(NFY_NOTICE, "Using exact transform ...");
}
} else {
throw xyzException
(xyz::tr
("Internal Error: xyzGDALWarpingReader used when not needed"));
}
// ***** initialize GDALWarpOptions *****
// see gdalwarper.h for explaination of options
GDALWarpOptions *options = GDALCreateWarpOptions();
xyzCallGuard<GDALWarpOptions*> optionGuard(GDALDestroyWarpOptions,
options);
options->eResampleAlg = resampleAlg;
options->eWorkingDataType = gdalDatatype;
options->hSrcDS = srcDS;
// pass our dummy, GDALWarpOperation validates that it exists, even
though,
// we'll never use it by calling, WarpRegionToBuffer */
options->hDstDS = (GDALDataset*)dstDS;
options->nBandCount = numBands;
options->panSrcBands = (int *)CPLMalloc(sizeof(int) * numBands);
options->panDstBands = (int *)CPLMalloc(sizeof(int) * numBands);
for (uint i = 0; i < numBands; ++i) {
options->panSrcBands[i] = gdalBandIndexes[i];
options->panDstBands[i] = gdalBandIndexes[i];
}
options->pfnTransformer = transformFunc;
options->pTransformerArg = transformArg;
if (noData != NULL) {
if (options->padfSrcNoDataReal == NULL) {
options->padfSrcNoDataReal = (double *)
CPLMalloc(sizeof(double) * options->nBandCount);
options->padfSrcNoDataImag = (double *)
CPLMalloc(sizeof(double) * options->nBandCount);
}
if (options->padfDstNoDataReal == NULL) {
options->padfDstNoDataReal = (double *)
CPLMalloc(sizeof(double) * options->nBandCount);
options->padfDstNoDataImag = (double *)
CPLMalloc(sizeof(double) * options->nBandCount);
}
for (uint i = 0; i < numBands; ++i) {
options->padfSrcNoDataImag[i] = 0.0;
options->padfDstNoDataImag[i] = 0.0;
options->padfSrcNoDataReal[i] = *noData;
options->padfDstNoDataReal[i] = *noData;
}
}
// ***** finally initialize the GDALWarpOperation object *****
if (warpOperation.Initialize(options) != CE_None) {
throw xyzException(xyz::tr("Unable to initialize GDAL warper"));
}
}
void
xyzGDALWarpingReader::FetchPixels(const xyzExtents<uint32> &srcExtents)
{
/*** I TRIED CHUNK AS WELL ****/
if ( warpOperation.ChunkAndWarpImage( srcExtents.beginX(),
srcExtents.beginY(),
GDALGetRasterXSize(dstDS),
GDALGetRasterYSize(dstDS)) !=
CE_None ) {
throw khException(kh::tr("Failed chunk and warp image"));
}
if (warpOperation.WarpRegionToBuffer(srcExtents.beginX(),
srcExtents.beginY(),
srcExtents.width(),
srcExtents.height(),
&rawReadBuf[0], /*****THIS ONE IS
EMPTY*****/
gdalDatatype) != CE_None) {
throw xyzException(xyz::tr("Failed to read source data"));
}
}
RasterBand ReadIO and ReadBlock :
//
****************************************************************************
// *** xyzVRRasterBand
//
****************************************************************************
xyzVRRasterBand::xyzVRRasterBand(GDALDataset *dataset,
int gdalBandNum,
const xyzSize<uint32> &rastersize,
const xyzSize<uint32> &blocksize,
GDALDataType gdalDatatype,
GDALColorInterp colorInterp_,
const xyzTransferGuard<GDALColorTable>
&ctbl_,
bool haveNoData_,
double noDataValue_) :
colorInterp(colorInterp_),
ctbl(ctbl_),
haveNoData(haveNoData_),
noDataValue(noDataValue_)
{
// initialize parts of my base class
// these really should be args to the base class constructor
poDS = dataset;
nBand = gdalBandNum;
nRasterXSize = rastersize.width;
nRasterYSize = rastersize.height;
eDataType = gdalDatatype;
eAccess = GA_ReadOnly;
nBlockXSize = blocksize.width;
nBlockYSize = blocksize.height;
}
CPLErr
xyzVRRasterBand::IRasterIO( GDALRWFlag eRWFlag,
int nXOff, int nYOff, int nXSize, int nYSize,
void *pData, int nBufXSize, int nBufYSize,
GDALDataType eBufType,
int nPixelSpace, int nLineSpace
// GDALRasterIOExtraArg*
psExtraArg
) {
if( poDS->RasterIO(eRWFlag,
nXOff, nYOff, nXSize, nYSize,
pData, nBufXSize, nBufYSize,
eBufType,
1, &nBand,
nPixelSpace, nLineSpace,
nLineSpace * nBufYSize) != CE_None) {
CPLError(CE_Failure, CPLE_AppDefined,
"Internal Error: RasterIO failed.");
return CE_Failure;
}
return CE_None;
}
CPLErr xyzVRRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff, void
*pImage) {
int nXOff = nBlockXOff * nBlockXSize,
nYOff = nBlockYOff * nBlockYSize,
nXSize = nBlockXSize,
nYSize = nBlockYSize;
if( nXOff + nXSize > nRasterXSize )
nXSize = nRasterXSize - nXOff;
if( nYOff + nYSize > nRasterYSize )
nYSize = nRasterYSize - nYOff;
int nPixelSpace = GDALGetDataTypeSize(eDataType) / 8;
int nLineSpace = nPixelSpace * nBlockXSize;
if (IRasterIO(GF_Read,
nXOff, nYOff, nXSize, nYSize,
pImage, nXSize, nYSize,
eDataType,
nPixelSpace, nLineSpace) == CE_None){
notify(NFY_NOTICE, "Raster IO successful");
return CE_None;
}
CPLError(CE_Failure, CPLE_AppDefined,
"Internal Error: RasterIO failed in block reading");
return CE_Failure;
}
--
View this message in context: http://osgeo-org.1560.x6.nabble.com/GDAL-Access-window-out-of-range-in-RasterIO-and-blank-buffer-in-warp-for-same-block-extents-tp5307794.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
More information about the gdal-dev
mailing list