如何从RasterBrick中提取数据？

我有一个RasterBrick,包含7年以上的月降雨量数据,所以它有7层,每层有12个插槽:

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 

由此我想提取在空间和时间上分布的点的降雨值.这些点位于数据框中:

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

我可以通过将数据帧转换为空间数据框并使用提取函数来处理空间方面:

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

但是,我无法确定如何确保从栅格砖中正确的栅格图层中提取降雨量值.我已经尝试了使用数据框中的"年"和"月"列进行索引的各种方法,但没有任何工作.任何提示将不胜感激!

这是我的第一篇文章,如果有太多/不够的信息,请道歉.如果看到更多我的代码会有用,请告诉我.

让我们在一个只有七个月的傻日历中，在三年内获取3x4栅格网格：

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

现在，让我们按年和月为砖层命名：

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"

并提出一些测试点：

> 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

现在为这些点构造图层名称，并与名称匹配：

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

现在我不认为层索引extract是矢量化的，因此您必须循环执行...

> 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

或在一个简单的向量中：

> 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

我会做一些检查，以确保这些值是您从这些输入中获得的期望值，然后再对任何主要问题进行处理。容易以错误的方式获取索引。

一次extract调用的另一种方法是计算所有图层的值，然后使用2列矩阵子集进行提取：

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

同样的数字，令人欣慰。