Raster operations in a raster list in R

Question

I have a raster list (21 in total) S1 that have the same projection, resolution but not exactly the same extent.

Due to this difference, I cannot stack them without solving the extent/dimensions issue. However, this would require the use of raster::resample and I want to avoid it due to the large processing time that it requires.

Those rasters cover my study area with a different extent and orientation (see image | gray area == study area)

My objective is to create a raster that covers the study area and where the pixel value is the maximum value of that pixel using as reference all the rasters of the raster list

Any experience with this? Is it mandatory to solve the extent issue and stack the rasters?

Answer 1

Depends on how fine the resolution of your output raster needs to be. If the number of pixels in your study area isn't too large, you could create a regular grid of vector points where you want your cell centers to be, and then use them with raster::extract(). Toy example:

library(sp)
library(raster)

# mock up a grid
grid <- GridTopology(c(0.5,0.5), c(1,1), c(5,5))
grid <- SpatialGrid(grid)
pt_grid <- as(grid, 'SpatialPoints')
pt_grid <- shift(pt_grid, 4, 6)

# mock up a couple of non-aligned rasters
r1 <- raster(matrix(sample(1:150, size = 150, replace = F), nrow = 15),
             xmn = 0, xmx = 10, ymn = 0, ymx = 15)
r2 <- raster(matrix(sample(1:100, size = 100, replace = F), nrow = 10),
             xmn = 3.3, xmx = 13.3, ymn = 5.2, ymx = 15.2)

plot(r1, alpha = 0.5)
plot(r2, alpha = 0.5, add = T)
plot(pt_grid, add = T, pch = 19)

# extract data from raster list
dat <- sapply(list(r1, r2), function(x) {
  raster::extract(x, pt_grid)
})

# calc max and append to points
pt_grid$maxval <- apply(dat, MARGIN = 1, FUN = max, na.rm = TRUE)

# cast points object to rasterstack
pt_px <- as(pt_grid, 'SpatialPixelsDataFrame')
pt_r <- raster::raster(pt_px)

If you want a fine-grained raster with a lot of cells, you are better off resampling, but you can do things like use system2() with GDAL to do the work more quickly. There are other Q&A's on this site that can walk you through how to do that.

Answer 2

Outline solution:

• Convert your output raster to a set of points
• Apply extract(r[[i]], output_points) for all list elements i to get a list of values of each input raster at your output points. Keep NA values outside the raster.
• Apply max (with na.rm=TRUE) over corresponding values of the list created.
• Build the final output raster from the output raster geometry and the max values

I'm slightly worried that you have rasters in different orientations, because that implies the raster is in a different coordinate reference system. You may have to reproject the output grid points to each raster's CRS and extract.