Calculating raster standard deviation in a memory efficient way in R

Question

I am trying to calculate standard deviation between two NDVI rasters stacked together using the raster::calc function in R. Both of these stacked rasters are 355 MB in size, and I have 4.2 GB of RAM available out of 16 GB. Now the calc function is giving me an error:

  Error: memory exhausted (limit reached?)

How can the memory issue be solved?

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

That is odd; it suggests that something is missed in assessing the memory needs. You can check for yourself what is going on under the hood with

canProcessInMemory(NDVI_Stack, verbose=TRUE)

And you can allow raster to use less memory, for example like this

rasterOptions(memfrac=.3)

(and/or using smaller chunk-sizes), see ?rasterOptions. And see the difference

canProcessInMemory(NDVI_Stack, verbose=T)

And try calc again.

Reducing the RAM that raster can use will slow it down, so you would not generally do that. In fact, others will increase it to speed things up.

None of this should normally be necessary. While there are some corner cases I am a bit surprised about this happening, and it would be useful if you edited your question to show what canProcessInMemory(NDVI_Stack, verbose=TRUE) and sessionInfo() return.

If the values of NDVI1 are in memory it may help to write it to disk so that more RAM stays available

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

Or in one swoop

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

By the way, you say that your rasters have a size of 355 MB, but I assume that this is the size on disk, which can be much smaller than the size in memory, because of compression and the use of, e.g., byte sized values instead of the 8 times larger (double) floating point values they become in R.

Finally, if you try the development version of terra, which you can install like this:

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

You can do

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

And compute the sample sd (denominator is n-1) like this:

NDVI_Std_sample = app(NDVI_Stack, fun = sd)

Or the population sd (denominator is n) like this:

NDVI_Std_pop = stdev(NDVI_Stack)