2

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)
  • 2
    What is the size of your raster files? You are probably running out of memory due to the big amount of data you are loading, hence the memory error. – chkaiser Sep 21 '15 at 9:55
  • as chkaiser mentioned - perhaps a problem because the rasters are in memory twice? (once during append and another during dstack?). maybe create a np array with the final shape and add each raster to the array as it's read? – fluidmotion Sep 21 '15 at 12:37
  • 1
    Hi chkaiser and fluidmotion, I tried to revised the script just to test if it really because of memory problem. I change the dstack argument to: stackRast = np.dstack((rasters[0], rasters[1], rasters[2], rasters[3], rasters[4], rasters[5], rasters[6], rasters[7], rasters[8], rasters[9])) and it worked. The shape of my array is (4800, 4800, 10) and if I extract an timeseries value for a particular index it returns an array([1719, 7068, 6958, 6822, 6214, 6534, 6972, 7715, 7113, 7441]). However, my method will not be possible if I will be dealing with approx. 800 rasters. Any advice? Thanks-Leo – user32145 Sep 22 '15 at 10:05
0

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

-3

build a gdal VRT file containing a pointer to their address

  • 1
    and how should this be done? your answers lack explanation (IMHO) – LaughU Jun 29 '18 at 11:41

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.