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

让我们在一个只有七个月的傻日历中，在三年内获取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

同样的数字，令人欣慰。