My intention is to stack my 2D array created from my MODIS rasters using numpy dstack. After that, I hope to implement savitzky-golay filtering algorithm from scipy signal processing so I think it is better to stack my 2D array to access the time series values of a particular pixel. Right now, I am just testing for 10 MODIS rasters (but later on will use my whole dataset -- approx 800 rasters (15 years of data)) only to test how should I create my script. So my initial workflow is to:
1) Access my raster files
2) For each file I need to convert it to a 2D array using RasterToNumPyArray
3) Append all those created 2D arrays to a list, and
4) Implement numpy.dstack() to stack my 2D arrays and hopefully I can begin some analysis.
I am getting an error in the 4th step (it says memory error) -- converting my list of 2D arrays to a stack (please see script 1 below) but if I revise my script and load in dstack using the index of each element in my 2D array it runs okay. Could you guide how to approach my problem. I am just building my script so for now its not so much. please see below. Thanks.
Script 1
import arcpy
from arcpy.sa import *
import os, sys
import numpy as np
from scipy.signal import savgol_filter
arcpy.CheckOutExtension("Spatial")
rasters = []
ws = 'G:/Test/Raster'
for folder, subs, files in os.walk(ws):
for filename in files:
aSrc = arcpy.RasterToNumPyArray(os.path.join(folder,filename))
rasters.append(aSrc)
stackRast = np.dstack(rasters)
Script 2
import arcpy
from arcpy.sa import *
import os, sys
import numpy as np
from scipy.signal import savgol_filter
arcpy.CheckOutExtension("Spatial")
rasters = []
ws = 'G:/Test/Raster'
for folder, subs, files in os.walk(ws):
for filename in files:
aSrc = arcpy.RasterToNumPyArray(os.path.join(folder,filename))
rasters.append(aSrc)
stack = np.dstack((rasters[0], rasters[1], rasters[2], rasters[3], rasters[4], rasters[5],rasters[6], rasters[7], rasters[8], rasters[9]))