pila raster usando GDAL y numpy dstack

Question

Estoy probando el filtro Savitzky-Golay con mi conjunto de datos MODIS. Actualmente estoy probando para 10 rasters (pero más adelante utilizaré todo mi conjunto de datos) sólo para probar cómo debo crear mi script. Así que mi flujo de trabajo inicial es:

1. Acceder a mis archivos rasterizados
2. cada archivo que necesito convertir en un array 2D usando GDAL o GDAL Numeric
3. Añade todas las matrices 2D creadas a una lista
4. Implementar numpy.dstack() para apilar mis matrices 2D y espero poder empezar algún análisis.

Sin embargo, estoy obteniendo un error en el cuarto paso (dice que hay un error de memoria) -- convirtiendo mi lista de arrays 2D en una pila. Podría guiarme con algún consejo sobre cómo proceder. Estoy construyendo mi script así que por ahora no es mucho.

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

Efectivamente, se trata de un MemoryError, y lamentablemente no hay mucho que se pueda hacer al respecto. Este es un problema común en Python cuando se trabaja con datos espaciales, ya que esas bibliotecas no soportan operaciones fuera del núcleo (lectura-escritura desde el disco en lugar de la memoria).

Yo sugeriría ver si algo así mejora su situación:

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])

No puedo probar esto por mí mismo, pero espero que la recolección de basura libere algo de memoria en cada bucle. Creo que puede haber una forma más eficiente de hacer stackRast = np.dstack([stackRast, aSrc]) pero no puedo pensar en ello en este momento

Answer 2

Construir un archivo gdal VRT que contenga un puntero a su dirección