Gdal Dataset.ReadAsArray() crashes Python

Question

I am using Python 2.6.5 (32bit) with Numpy 1.3 and Gdal 1.9.1 installed on Windows 7 64bit. I am trying to read an 800 MB Imagine (.img) raster dataset into a Numpy array to do some raster algebra, but as soon as I run the following code, Python.exe crashes.

from osgeo import gdal

g = gdal.Open(r'path\to\dataset', gdal.GA_Readonly)
b = g.GetRasterBand(1)
data = b.ReadAsArray()

Python.exe crashes on the b.ReadAsArray() call. I did some Google searching and found articles dated from Gdal 1.6 that mentioned this issue with Windows 7 64bit, but they also mentioned that it had been fixed in the latest development versions at that time.

Has anyone else had this issue? Any workarounds?

UPDATE:

I decided to debug the code in PyDev to try to determine where it is failing. From what I can tell (still not getting any error messages) it is failing on line 22 of gdal_array.py.

_mod = imp.load_module('_gdal_array', fp, pathname, description)

When I step into the line of code above, it brings me into the init.py module of numpy. When I reach the end of the numpy.__init__.py module, it steps out back to the above line of code. Then, when I hit the step into button, which should bring me to the next line within gdal_array.py, the script just terminates with no error messages or anything.

UPDATE #2:

I uninstalled GDAL 1.9.1 and installed GDAL 1.6.1 from the Python Cheeseshop and OSGeo's Windows Binaries. Still had the same issue.

Answer 1

As suspected by many commenters, this was an issue with my install. Apparently, i was not paying close enough attention when installing GDAL and the Python Bindings.

I installed GDAL Core and plugins (dll's) from gisinternals.com, but somehow I didn't think to install the Python Bindings from there as well. The Python Bindings I installed were from a different site (can't remember which one at this point).

When I reinstalled GDAL and Python Bindings all from gisinternals.com, I was able to successfully ReadAsArray.

Thank you to all who commented and answered and I apologize for my ignorance.

Answer 2

It's possible it's a memory issue. When you use ReadAsArray it brings the data into memory, and while 800mb isn't massive, it's not tiny either. Have you tried reading the array in chunks?

data = b.ReadAsArray(x_offset, y_offset, x_size, y_size)

You should be able to loop through the array and process it chunk at a time, though depending on the processing you're doing you should probably look to reading in areas with overlap to avoid edge effects.

Answer 3

Sorry I'm late to this party, but your basic problem is that 32-bit Python can't store very large rasters in memory. You can read your large raster into memory in bite-size chunks, but then you're fairly limited in terms of what you can process efficiently without extremely inefficient/frequent read/writes to disk.

What I do instead (which sacrifices some efficiency due to disk read/writes) is to call (via EXE-wrapping) the 64-bit gisinternals.com version of whatever gdal method you need. Be careful with using python's subprocess module in a loop (i.e. you may want/need to call subprocess sequentially) as you may inadvertently generate too many open threads for your windows box and get ominous system warnings. You sacrifice a little bit in the way to disk read/writes with this gdal approach, but your processing efficiency only drops (i.e. relative to a lightening fast in-memory computation, if your box/library could support it) by about a factor or 10.