从 rasterBrick 或 rasterStack 中提取值时出错答案

【问题标题】：Error when extracting values from a rasterBrick or rasterStack从 rasterBrick 或 rasterStack 中提取值时出错
【发布时间】：2015-06-11 01:01:42
【问题描述】：

我无法从 rasterStack 或 rasterBrick 类的多波段栅格中提取值或点。 'extract' 适用于单个栅格，但在应用于 rasterStack 或 brick 时会出现错误。

> all.var
class       : RasterBrick 
dimensions  : 89, 180, 16020, 34  (nrow, ncol, ncell, nlayers)
resolution  : 2, 2  (x, y)
extent      : -179, 181, -89, 89  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       :    period_max, pct_signif_periodmax, pct_coi_periodmax,    pct_ispos_signif, events_pos_periodmax, events_neg_periodmax, events_pos_all, events_neg_all,  maxpower_pos,  maxpower_neg, maxpower_events_pos, maxpower_events_neg, maxpower_pos_norm, maxpower_neg_norm, maxpower_events_pos_norm, ... 

> point
Lon      Lat
1  166.2790 -10.2690
2   26.9000 -33.6000
3  153.6209 -28.7001
4  113.8333 -28.6833
5  153.6335 -28.6591
6  153.5836 -28.4643
7   32.6833 -27.5333
8   32.6880 -27.5260
9   32.6880 -27.5260
10  32.6880 -27.5260


> point.extract<-extract(all.var, point, buffer=50000,na.rm=TRUE,fun=mean)
Error in apply(x, 2, fun2) : dim(X) must have a positive length

这适用于单个栅格，但会在堆栈/砖块上失败，并且仅在我使用缓冲区参数时才会引发错误。

【问题讨论】：

标签： r r-raster

【解决方案1】：

这是一个说明错误的工作 R 示例：

library(raster)
b <- brick(nrow=89, ncol=180, nl=34, xmn=-179, xmx=181, ymn=-89, ymx=89, crs="+proj=longlat +datum=WGS84")
b[] <- 1

p <- matrix(c(166.2790,-10.2690,26.9000,-33.6000,153.6209,-28.7001,113.8333,-28.6833,153.6335,-28.6591,153.5836,-28.4643,32.6833,-27.5333,32.6880,-27.5260,32.6880,-27.5260,32.6880,-27.5260), ncol=2, byrow=TRUE)

v <- extract(b, p, buffer=50000, na.rm=TRUE, fun=mean)

这确实给出了您报告的错误，可能是由于光栅包中的错误。这是一个解决方法：

v <- extract(b, p, buffer=15000000)
# get the mean for each point (buffer) by layer
vv <- lapply(v, function(x) ifelse(is.matrix(x), colMeans(x, na.rm=TRUE), x))
# combine
do.call(rbind, vv)

【讨论】：

罗杰。这种解决方法效果很好。但是，输出（即 v）是一个向量，而不是代表每一层的 34 列的数据帧/矩阵。我错过了什么吗？
这为您提供了每个点的平均值（因此向量的长度为 10）。要获得每一层的平均值，您可以执行 rowMeans(do.call(cbind, v))
不完全。这得到了全局行的意思。理想情况下，输出应该是 10 行（xy 坐标）和 34 列（层）维度的矩阵
这不会得到缓冲区中每个点和层的提取值的平均值（即原始代码中的 fun=mean）。
....例如，如果您将上例中的缓冲区大小增加到 150000m，则行尺寸不是 10