Error combining multiple spatRaster files along time dimension

Question

I am working with temperature data from PRISM. Thus far, I have been able to download data for individual days, then crop to shape files, and extract information accordingly. However, I am looking to combine multiple spatRaster files along the time dimension, so I end up with one data frame.

Thus far, I have tried two approaches. Firstly, I have tried combining my two days' worth of data into a single .rds file (stars object), then converting to spatRaster:

file1<- readRDS("file1.rds")
file2<- readRDS("file2.rds")
combinedfiles<- c(file1, file2, along = 3)

the combined file has the following description/attributes:

stars object with 3 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
                                  Min. 1st Qu. Median    Mean 3rd Qu.  Max.  NA's
PRISM_tmin_stable_4kmD2_202...  -1.447    8.65 10.344 11.7053  16.168 20.97 60401
dimension(s):
        from   to offset    delta refsys x/y
x          1 1405   -125  0.04167  NAD83 [x]
y          1  621  49.94 -0.04167  NAD83 [y]
new_dim    1    2     NA       NA     NA    

however, I cannot seem to convert this to a spatRaster successfully:

temprast<-stack(combinedfiles)
Error in stack.default(combinedfiles) : 
  at least one vector element is required

I also tried converting each individual day to a spatRaster, then combining- this also did not work:

temprast<-rast(file1)
temprast2 <-rast(file2)

raster_brick <- brick(temprast, temprast2)
Error in .local(x, ...) : 
  unused argument (new("SpatRaster", pnt = new("Rcpp_SpatRaster", .xData = <environment>)))

The ultimate goal is to be able to crop this raster file to a shape file. It would be much easier if I combined the .rds files first in chronological order, before then dealing with shape file cropping etc. but if necessary, I guess my last option is to individually crop each .rds file to a shape file, extract data, and then combine (potentially using bind_rows)? But this seems very inefficient.

Answer 1

You are using stack where you intend to use raster::stack You are using raster::brick where you should use a function from "terra", namely c()

A sensible approach could be

ff <- list.files(pattern=".bil$")
library(terra)
x <- terra::rast(ff)

Answer 2

PRISM files are all structured identically: raw binary 32-bit numeric data (the BIL format) of 621 rows and 1405 columns. This is a format typically used by GIS software (from way back when) and both the terra and stars packages can deal with that format easily. There is, however, a lot of overhead going on in both packages and no good support for the date information contained in the file name. Since you are running into errors, you might be better off with a bare-bones base R approach using readBin() and then set some additional information yourself. After merging all the files you can import the combined data in terra or stars and do the cropping and further analysis.

As a function, merging all the data files in a directory would look like this:

mergePRISM <- function(dir) {
  # Need package abind
  requireNamespace("abind")
    
  # The files to merge
  lf <- list.files(dir, "\\.bil$", full.names = TRUE)
  
  # The dates of each of the files
  dates <- sapply(strsplit(basename(lf), "_"), function(p) p[5])
  
  # Read the data, then rearrange to go from row-major (BIL) to column-major (R)
  data <- lapply(lf, function(fn) {
    d <- readBin(fn, what = "numeric", n = 621 * 1405, size = 4)
    dim(d) <- c(1405, 621)
    aperm(d, c(2, 1))
  })
  
  # Merge the data
  data <- abind(data, along = 3)
  
  # Set the dimension names
  lon <- seq(from = -125, to = -125 + 1404 * 0.0416666666667, by = 0.0416666666667)
  lat <- seq(from = 49.9166666666664, to = 49.9166666666664 - 620 * 0.0416666666667, by = -0.0416666666667)
  dimnames(data) <- list(lat, lon, dates)
  data
}

Note that you might have to set the endian argument of readBin() if your data is garbled (on my Macbook I don't have to).

You can do this for a month's worth of data or even a year but I wouldn't make the time-series longer than that because of RAM issues.

You can turn the result into a stars object with st_as_stars() and then further process, such as cropping to a shapefile.