GDAL : Access window out of range in RasterIO and blank buffer in warp for same block extents

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

GDAL : Access window out of range in RasterIO and blank buffer in warp for same block extents

asmita
This post was updated on .
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 xyzException(xyz::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;

}




Reply | Threaded
Open this post in threaded view
|

Re: GDAL : Access window out of range in RasterIO and blank buffer in warp for same block extents

Johan de Braak
This post was updated on .
asmita wrote
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.
(...)
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?
I'm not sure if I understand fully what is going on, but the point of the error message is that you are trying to access more pixels than are available in your raster. What helps, in general, is to think of what the extent of the project you are working on is, and tell gdal (using -te or -projwin) what extent that is. If your input rasters have different extents / projections, use gdalwarp to fix this, before continuing with mosaic creation.