r - How to extract data from a RasterBrick?

Question

I have a RasterBrick consisting of monthly rainfall data over 7 years, so it has 7 layers with 12 slots each:

rainfall <- brick("Rainfall.tif")
    > rainfall
    class       : RasterBrick
    dimensions  : 575, 497, 285775, 7  (nrow, ncol, ncell, nlayers)
    resolution  : 463.3127, 463.3127  (x, y)
    extent      : 3763026, 3993292, -402618.8, -136213.9  (xmin, xmax, ymin, ymax)
    coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs 
    data source : in memory
    names       : layer.1.1, layer.2.1, layer.1.2, layer.2.2,   layer.1,   layer.2,     layer 
    min values  :  239.6526,  499.8343,  521.0316,  617.2896,  596.0397,  663.6633,  298.0572 
    max values  :  691.9075, 1158.2064, 1184.9858, 1198.7121, 1241.8077, 1114.7598,  832.6042 

From this I would like to extract a value for rainfall at points distributed both spatially and temporally. These points are in a data frame:

points <- read.csv("Points.csv")
    > head(points)
        ID      x          y      ncell  jday  FRP_max    FRI   year   month
       69211  3839949  -171684.6    17    59      NA  230.2500  2001     2
       69227  3808720  -238808.7    16    52      NA        NA  2001     2
       69237  3793373  -267563.1     1    52      NA        NA  2001     2
       69244  3986574  -292118.7     1    43      NA        NA  2001     2
       32937  3864736  -164296.8   106    77    94.8  249.1524  2001     3
       32938  3871463  -163123.4    31    82      NA  253.5081  2001     3

I can handle the spatial aspect by converting the data frame to a spatial data frame and using the extract function:

points.sp <- points
coordinates(points.sp) <- ~ x + y
rainfall.points <- extract(rainfall, points.sp)

However, I can't work out how to make sure the rainfall values are being extracted from the correct raster layer from within the raster brick. I've tried various ways of indexing using the "year" and "month" columns from my data frame but nothing has worked. Any tips would be much appreciated!

This is my first post so apologies if there's too much/not enough info. Let me know if seeing more of my code would be useful.

Answer 1

lets take a grid of 3x4 rasters over three years in a silly calendar that only has seven months in it:

d = array(1:(3*4*7*3),c(3,4,7*3))
b = brick(d)

Now lets give the brick layers names by year and month:

names(b) = paste("rain",outer(1:7,2001:2003,paste,sep="-"),sep="-")
> names(b)
 [1] "rain.1.2001" "rain.2.2001" "rain.3.2001" "rain.4.2001" "rain.5.2001"
 [6] "rain.6.2001" "rain.7.2001" "rain.1.2002" "rain.2.2002" "rain.3.2002"
[11] "rain.4.2002" "rain.5.2002" "rain.6.2002" "rain.7.2002" "rain.1.2003"
[16] "rain.2.2003" "rain.3.2003" "rain.4.2003" "rain.5.2003" "rain.6.2003"
[21] "rain.7.2003"

and make some test points:

> pts = data.frame(x=runif(3),y=runif(3), month=c(5,1,3),year = c(2001,2001,2003))
> pts
          x         y month year
1 0.2513102 0.8552493     5 2001
2 0.4268405 0.3261680     1 2001
3 0.7228359 0.7607707     3 2003

Now construct the layer name for the points, and match to the names:

pts$layername = paste("rain",pts$month,pts$year,sep=".")
pts$layerindex = match(pts$layername, names(b))

Now I don't think the layer index in extract is vectorised, so you have to do it in a loop...

> lapply(1:nrow(pts), function(i){extract(b, cbind(pts$x[i],pts$y[i]), layer=pts$layerindex[i], nl=1)})
[[1]]
     rain.5.2001
[1,]          57

[[2]]
     rain.1.2001
[1,]           5

[[3]]
     rain.3.2003
[1,]         201

Or in a simple vector:

> sapply(1:nrow(pts), function(i){extract(b, cbind(pts$x[i],pts$y[i]), layer=pts$layerindex[i], nl=1)})
[1]  57   5 201

I'd do some checks to make sure those values are what you expect from those inputs before doing it on anything major though. Its easy to get indexes the wrong way round....

Another way to do it with a single extract call is to compute the values for all layers and then extract with a 2-column matrix subset:

> extract(b, cbind(pts$x, pts$y))[
      cbind(1:nrow(pts),match(pts$layername, names(b)))
     ]
[1]  57   5 201

Same numbers, comfortingly.