To achieve this task, we can use the calc
function from the raster package. We also need to know how to subset the RasterBrick
object.
Data preparation
library(raster)
set.seed(123)
r <- raster(ncol=10, nrow=10)
r_brick <- brick(sapply(1:45, function(i) setValues(r, rnorm(ncell(r), i, 3))))
Example of calc function
calc
can apply a function for all layers in a RasterBrick
object. The end result is a raster layer.
# Calculate mean
r_mean <- calc(r_brick, mean)
# Calculate median
r_median <- calc(r_brick, median)
# Calculate sd
r_sd <- calc(r_brick, sd)
Notice that r_mean
, r_median
, and r_sd
are all RasterLayer
.
Example of subset a RasterBrick
We can use the index to subset the layer. For example,
r_sub <- r_brick[[1:3]]
r_sub
is the first three layers of r_brick
A function to Conduct the Analysis
Knowing the technique of calc
and subset, we can design a function to conduct the analysis.
The first thing to do is create a vector serving as a reference to year and index.
# Create the index
ind <- 1:45
names(ind) <- 1971:2015
Calling the year number to ind
will return the index. For example,
# Get the index of 2015
ind[as.character(2015)]
#2015
# 45
Now design the function, which has five arguments
end_year
: The end year of analysis
n_year
: The last n
year in terms of the end year
FUN
: A function, such as mean
, median
, and sd
index
: The year index (ind
)
ras_brick
: RasterBrick
to work with
# Define the function
raster_stat <- function(end_year, n_year, FUN, index, ras_brick){
# Subset the raster
index_temp <- index[as.character((end_year - n_year + 1):end_year)]
ras_brick_temp <- ras_brick[[index_temp]]
# Calculate the statistics
ras_result <- calc(ras_brick_temp, FUN)
# Set the name
names(ras_result) <- paste("Y", end_year, n_year, substitute(FUN), sep = "_")
return(ras_result)
}
Now we can test the function.
raster_stat(2015, 5, FUN = sd, index = ind, ras_brick = r_brick)
#class : RasterLayer
#dimensions : 10, 10, 100 (nrow, ncol, ncell)
#resolution : 36, 18 (x, y)
#extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
#data source : in memory
#names : Y_2015_5_sd
#values : 12.05333, 14.61298 (min, max)
Notice that the result of the raster_stat
function has the name Y_2015_5_sd
. This is helpful to identify which end_year
, n_year
, and FUN
was applied.
Apply the function
We can use for loop to apply raster_stat
through all end_year
and n_year
. Here is an example calculating the mean
.
# Set the range of end_year and n_year
end_year_vec <- 2000:2015
n_year_vec <- c(5, 10 , 15, 20, 30)
# Create an empty list to store result
r_mean_list <- list()
for (i in end_year_vec){
for(j in n_year_vec){
result_temp <- raster_stat(end_year = i, n_year = j, FUN = mean,
index = ind, ras_brick = r_brick)
# Add the raster layer to the result_list
r_mean_list[[names(result_temp)]] <- result_temp
}
}
All the results are stored in r_mean_list
with a unique name. We can use the same approach for median
and sd
.