I'm trying to run Mort Canty's http://mcanty.homepage.t-online.de/ Python iMAD implementation on bitemporal RapidEye Multispectral images. Which basically calculates the canonical correlation for the two images and then substracts them. The problem I'm having is that the images are of 5000 x 5000 x 5 (bands) pixels. If I try to run this on the whole image I get a memory error.

Would the use of something like pyTables help me with this?

What Mort Canty's code tries to do is that it loads the images using gdal and then stores them in an 10 x 25,000,000 array.

# initial weights
wt = ones(cols*rows)      
# data array (transposed so observations are columns)
dm = zeros((2*bands,cols*rows))
k = 0
for b in pos:
band1 = inDataset1.GetRasterBand(b+1)
band1 = band1.ReadAsArray(x0,y0,cols,rows).astype(float)
dm[k,:] = ravel(band1)
band2 = inDataset2.GetRasterBand(b+1)
band2 = band2.ReadAsArray(x0,y0,cols,rows).astype(float)        
dm[bands+k,:] = ravel(band2)
k += 1

Even just creating a 10 x 25,000,000 numpy array of floats throws a memory error. Anyone have a good idea of how to get around this?

  • which operating system do you use? – Aragon May 14 '12 at 18:28
  • 4
    Some attractive workarounds: (1) read only pairs of bands at a time. This reduces RAM requirements 60% but doubles the amount of input. (2) Subsample each image either randomly or systematically. E.g., taking every other row and every other column will reduce RAM requirements by 75% but will barely change the precision of the correlation coefficients. (3) Often the original bands are just bytes, so if you can, read them as bytes rather than floats. This reduces RAM requirements by 75%. (4) Install more RAM! – whuber May 14 '12 at 21:14

You can try to use PyTables along with numpy. With PyTables you can store some of your temporary arrays in HDF format on disk during processing.

You can also use some tips from interesting document The NumPy array: a structure for efficient numerical computation

Another way is to optimize your algorithm in some way - to force garbage collection. See here.

| improve this answer | |

The tool "gdalwarp" from GDAL can do that well:

Shell style:


# in megabytes
CACHE="--config GDAL_CACHEMAX 8000 -wm 8000"

LIST="`ls ??????.tif`"

# fix some 'no data' on the fly
# for reprojection, using -r and -tr (don't use otherwise)
gdalwarp $CACHE $BIGTIFF -srcnodata 65535 -dstnodata 65535 \
         -r bilinear -tr $RES $RES $LIST out_$RES.tif
gdaladdo out_$RES.tif 2 4 8 16 32

Python style:

See these scripts:

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.