[gdal-dev] ReadAsArray misbehaving...

[Even Rouault]

Riann,

hum, I think you've hit a classical problem with the use of the GDAL python 
bindings. The root of your problem is likely the following line:

 reader = RasterReaderWriter(gdal.Open(iFilename).GetRasterBand(1))

Currently the GDAL python bindings do a poor job with life cycle management of 
native underlying objects. When chaining the calls, some objects can go out-
of-scope and are then destroyed. Specifically, the band object doesn't keep the 
dataset object alive, so the RasterReaderWriter() constructor will get a band 
object with a ghost dataset. This can typically result in application crashes, 
or as you've observed various corruptions.

All in all, I'd advise you rewriting the above as :

ds = gdal.Open(iFilename)
band = ds.GetRasterBand(1)
reader =  RasterReaderWriter(band)
....
.....
ds = None # we can drop the reference on the dataset when we don't need to use 
the reader.

An alternative if you don't want to keep the reference on the dataset in the 
calling code, is to pass both the dataset and the band to RasterReaderWrite() 
so that it can stores them in member variables of the class.

Best regards,

Even

Le vendredi 03 septembre 2010 12:39:09, Riaan van den Dool a écrit :
> I have recently started a new open-source project that will use gdal
> extensively, named scikits.eartho.
> 
> The vision is to implement some advanced algorithms and ideas that we
> develop and work with at the South African Satellite Applications Centre
> (SAC) for use in Python (SciPy).
> 
> I am experiencing some trouble with RasterBand.ReadAsArray.
> 
> My python code:
> 
> import osgeo.gdal as gdal
> 
> class RasterReaderWriter:
>     def __init__(self,  raster):
>         self.raster = raster
>         self.tileWidth =1024
>         self.tileHeight = 1024
>         self.tileNum = 0
> 
>     def setTileSize(self,  width,  height):
>         self.tileWidth =width
>         self.tileHeight = height
> 
>     def resetIterator(self):
>         self.tileNum = 0
> 
>     def readNextTile(self):
>         tileRaster,  offsetX,  offsetY,  width,  height  =
> self.readTile(self.tileNum)
>         self.tileNum = self.tileNum + 1
>         return tileRaster,  offsetX,  offsetY,  width,  height
> 
>     def readTile(self,  tileNum):
>         tilesInWidth = self.raster.XSize / self.tileWidth + 1
>         offsetX = tileNum % tilesInWidth * self.tileWidth
>         offsetY = tileNum / tilesInWidth * self.tileHeight
>         width = min(self.tileWidth,  self.raster.XSize - offsetX)
>         height = min(self.tileHeight,  self.raster.YSize - offsetY)
>         print offsetX,  offsetY,  width,  height
>         tileRaster = self.raster.ReadAsArray(offsetX,  offsetY, width,
> height)
>         return tileRaster,  offsetX,  offsetY,  width,  height
> 
> 
> 
> 
> reader = RasterReaderWriter(gdal.Open(iFilename).GetRasterBand(1))
> print reader.raster.XSize
> print reader.raster.YSize
> tile, offsetX, offsetY, width, height = reader.readNextTile()
> 
> Before the call to readNextTile() the raster XSize and YSize is correct.
> After the call the raster seems to have been corrupted as the XSize and
> YSize have invalid values.
> 
> Any help would be appreciated!
> 
> Riaan