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.
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.