3

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)
3
  • 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, 2015 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? Sep 21, 2015 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, 2015 at 10:05

2 Answers 2

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
  • 1
    and how should this be done? your answers lack explanation (IMHO)
    – LaughU
    Jun 29, 2018 at 11:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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