Search code examples
rimageimage-processingrasterr-raster

R- Raster math while preserving integer format

I have a some large rasters (~110 MB each) I want to perform some raster calculations on. For the purposes of this example, I want to average the files SNDPPT_M_sl1_1km_ll.tif and SNDPPT_M_sl2_1km_ll.tif, available at this website. In reality, the math is a bit more complex (some multiplication and division of several rasters).

Both input rasters are integer (INT1U) data, and I would like the output to also be INT1U. However, whenever I try to do a raster calculation, it creates intermediate temporary files in floating point format which are very large in size. I am working on a laptop with about 7 GB of free hard drive space, which gets filled before the calculation is complete.

# load packages
require(raster)

## script control
# which property?
prop <- "SNDPPT"

# load layers
r.1 <- raster(paste0("1raw/", prop, "_M_sl1_1km_ll.tif"))
r.2 <- raster(paste0("1raw/", prop, "_M_sl2_1km_ll.tif"))

# allocate space for output raster - this is about 100 MB (same size as input files)
r.out <- writeRaster(r.1, 
                     filename=paste0("2derived/", prop, "_M_meanTop200cm_1km_ll.tif"), 
                     datatype="INT1U")

# perform raster math calculation
r.out <- integer(round((r.out+r.2)/2))  

# at this point, my hard drive fills due to temporary files > 7 GB in size

Is anyone aware of a workaround to perform raster math in R with integer input and output files while minimizing or avoiding the very large intermediate files?

Solution

  • The trick here could be to use raster::overlay to make the computation and save the results as a compressed tiff at the same time. Something like this should work: 

    library(raster)
#> Loading required package: sp
# load layers
r.1 <- raster("C:/Users/LB_laptop/Downloads/SNDPPT_M_sl1_1km_ll.tif")
r.2 <- raster("C:/Users/LB_laptop/Downloads/SNDPPT_M_sl1_1km_ll.tif")

out <- raster::overlay(r.1, r.2,
                       fun = function(x, y) (round((x + y) / 2)),
                       filename = "C:/Users/LB_laptop/Downloads/SNDPPT_out.tif", 
                       datatype = "INT1U", 
                       options  = "COMPRESS=DEFLATE")
> out
class       : RasterLayer 
dimensions  : 16800, 43200, 725760000  (nrow, ncol, ncell)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : -180, 180, -56.00083, 83.99917  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : C:\Users\LB_laptop\Downloads\SNDPPT_out.tif 
names       : SNDPPT_out 
values      : 0, 242  (min, max)

    HTH.