stack raster using GDAL and numpy dstack

Question

I am testing Savitzky-Golay filter with my MODIS dataset. Currently I am testing for 10 rasters (but later on will use my whole dataset) only just to test how should I create my script. So my initial workflow is to:

1. Access my raster files
2. each file I need to convert to 2D array using GDAL or GDAL Numeric
3. Append all those created 2D arrays to a list
4. Implement numpy.dstack() to stack my 2D arrays and hopefully I can begin some analysis.

However, I getting error in the 4th step (it says memory error) -- converting my list of 2D arrays to a stack. Could you guide me some advice how to proceed. I am just building my script so for now its not so much.

import os, sys
import gdal, numpy as np
import gdalnumeric as gd
from scipy.signal import savgol_filter

rasters = list()

ws = 'G:/Test/Raster'
for folder, subs, files in os.walk(ws):
    for filename in files:
        aSrc = gd.LoadFile(os.path.join(folder,filename))
        rasters.append(aSrc)

stackRast = np.dstack(rasters)

Answer 1

It is indeed a MemoryError, and there is not much you can do about it unfortunately. This is a common problem in Python when working on spatial data, as those libraries don't support out-of-core operations (read-write from disk instead of memory).

I would suggest to see if something like this improves your situation:

import os, sys
import gdal, numpy as np
import gdalnumeric as gd
from scipy.signal import savgol_filter

ws = 'G:/Test/Raster'
try:
    del stackRast
except NameError:
    pass

for folder, subs, files in os.walk(ws):
    for filename in files:
        aSrc = gd.LoadFile(os.path.join(folder,filename))
        if 'stackRast' not in locals():
            stackRast = aSrc
        else:
            stackRast = np.dstack([stackRast, aSrc])

I cannot test this myself, but hopefully the junk collection will free up some memory every loop. I reckon there may be a more efficient way to do stackRast = np.dstack([stackRast, aSrc]), but cannot think of it at the moment

Answer 2

build a gdal VRT file containing a pointer to their address