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.

Cruzeiro

Where's the "12 slots" in your brick? I see 7 layers. – Splat

Yes, exactly what @Splat asked. Each of those 7 layers (they're the ones named layer.1.1, layer.2.1, etc.) will only have one value per pixel/cell, so it sounds like you're not quite clear on what data you've actually got in the brick. – Angiosperm

I did wonder about this but when I expand rainfall in the global environment it says: Large RasterBrick (2000425 elements, 15.3 mb) .. @ file :Formal class ' .RasterFile' [package "raster"] with 12 slots – Cruzeiro

"Formal class 'RasterBrick' [package "raster"] with 12 slots"? The "slots" there are R's internal structure parts - they store the resolution, the coordinates, the layer names etc. Its a coincidence it matches 12 months in a year! Your brick only has 7 layers, and each layer is a single number. – Splat

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.

Splat

Thanks! I think I'm being hindered from getting what I want by only having years, rather than months, in my raster stack so I'm going to try and solve that problem. – Cruzeiro

Answer 2

Each layer in your RasterBrick will have a unique name, so you can use match('name', names(b)) to find the numeric index of the layer you're interested in. Then use the layer= argument to extract() to point to the layer from which you'd like to extract (setting nl=1 to indicate that you only want that one layer).

Here's a reproducible example, in which I use cell numbers to do the extracting. (This will work exactly the same way when using SpatialPoints to indicate which values are to be grabbed.)

 ## An example SpatialBrick
b <- brick(system.file("external/rlogo.grd", package="raster"))
nlayers(b)
# [1] 3
names(b)
# [1] "red"   "green" "blue"

## Extract data from given cells in the "green" layer, 
ii <- match("green", names(b))
extract(b, 1000:1003, layer=ii, nl=1)
#      green
# [1,]   254
# [2,]   255
# [3,]   255
# [4,]   255

Galasyn