如何从 RasterBrick 中提取数据？

Question

我有一个 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)

但是，我无法确定如何确保从栅格砖内的正确栅格图层中提取降雨值。我尝试了使用数据框中的“年”和“月”列进行索引的各种方法，但没有任何效果。任何提示将不胜感激！

这是我的第一篇文章，如果信息过多/不足，敬请见谅。让我知道查看更多我的代码是否有用。

Answer 1

RasterBrick 中的每一层都有一个唯一的名称，因此您可以使用match('name', names(b)) 查找您感兴趣的层的数字索引。然后使用layer= 参数指向extract()到您要从中提取的层（设置nl=1 表示您只需要该层）。

这是一个可重现的示例，其中我使用单元格编号进行提取。 （这与使用SpatialPoints 指示要抓取哪些值时完全相同。）

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

Answer 2

让我们在一个只有七个月的愚蠢日历中在三年内采用 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

同样的数字，令人欣慰。