如何从 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)

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

这是我的第一个 post 如果有太多 much/not 足够的信息，我们深表歉意。让我知道查看更多我的代码是否有用。

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

相同的数字，令人欣慰。

如何从 RasterBrick 中提取数据？

How to extract data from a RasterBrick?

r

raster

spatial

r-raster