如何从 RasterBrick 中提取数据? [英] How to extract data from a RasterBrick?
问题描述
我有一个由 7 年的月降雨数据组成的 RasterBrick，所以它有 7 层，每层有 12 个槽:
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.
推荐答案
让我们在一个只有七个月的愚蠢日历中使用三年内的 3x4 栅格网格:
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"
并提出一些测试点:
> 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))
现在我不认为 extract 中的层索引是矢量化的，所以你必须在循环中进行...
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....
另一种通过单个 extract 调用完成的方法是计算所有层的值，然后使用 2 列矩阵子集进行提取:
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.
这篇关于如何从 RasterBrick 中提取数据?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
