Cálculo de la desviación estándar de un raster de forma eficiente en memoria en R

Question

Estoy tratando de calcular standard deviation entre dos NDVI rasters apilados utilizando el raster::calc función en R . Estos dos rásteres apilados son 355 MB en tamaño, y tengo 4.2 GB de RAM disponible fuera de 16 GB . Ahora el calc me da un error:

  Error: memory exhausted (limit reached?)

¿Cómo se puede resolver el problema de la memoria?

library(raster)

NDVI1 =  raster("path/NDVI1.tif")
NDVI2 =  raster("path/NDVI2.tif")

NDVI1= resample(NDVI1, NDVI1)

NDVI_Stack = stack(NDVI1, NDVI2)

NDVI_Std = calc(NDVI_Stack, fun = sd)

Answer 1

Eso es impar; sugiere que algo se ha perdido en la evaluación de las necesidades de memoria. Puedes comprobar por ti mismo lo que ocurre bajo el capó con

canProcessInMemory(NDVI_Stack, verbose=TRUE)

Y puedes permitir que raster para usar menos memoria, por ejemplo, así

rasterOptions(memfrac=.3)

(y/o utilizando tamaños de trozos más pequeños), véase ?rasterOptions . Y ver la diferencia

canProcessInMemory(NDVI_Stack, verbose=T)

E intenta calc de nuevo.

Reducir la RAM que puede utilizar el raster lo ralentizará, por lo que generalmente no lo harías. De hecho, otros la aumentarán para acelerar las cosas.

Normalmente, nada de esto debería ser necesario. Aunque hay algunos casos de esquina, me sorprende un poco que esto ocurra, y sería útil que editaras tu pregunta para mostrar lo que canProcessInMemory(NDVI_Stack, verbose=TRUE) y sessionInfo() volver.

Si los valores de NDVI1 están en la memoria, puede ser útil escribirlo en el disco para que haya más RAM disponible

 NDVI1 <- writeRaster(NDVI1, "temp.tif")

O de un solo golpe

 NDVI1 = resample(NDVI1, NDVI1, filename="temp.tif")

Por cierto, dices que tus rásteres tienen un tamaño de 355 MB, pero supongo que se trata del tamaño en disco, que puede ser mucho menor que el tamaño en memoria, debido a la compresión y al uso de, por ejemplo, valores de tamaño byte en lugar de los valores de punto flotante 8 veces mayores (dobles) en que se convierten en R.

Por último, si pruebas el versión de desarrollo de terra que puedes instalar así:

 `install.packages('terra', repos='https://rspatial.r-universe.dev')

Puedes hacer

library(terra)
NDVI1 =  rast("path/NDVI1.tif")
NDVI2 =  rast("path/NDVI2.tif")
NDVI1 = resample(NDVI1, NDVI1)
NDVI_Stack = c(NDVI1, NDVI2)

Y calcular el muestra sd (el denominador es n-1 ) así:

NDVI_Std_sample = app(NDVI_Stack, fun = sd)

O el población sd (el denominador es n ) así:

NDVI_Std_pop = stdev(NDVI_Stack)