将砖块中的每个栅格图层与数据帧 R 中的值相关联答案

【问题标题】：Correlate each raster layer in a brick with values in a dataframe R将砖块中的每个栅格图层与数据帧 R 中的值相关联
【发布时间】：2018-11-06 12:57:03
【问题描述】：

我在 R 中有这个数据框：

library(raster)    

# create a random dataframe with yearly values for each column
df <- data.frame(year = seq(1981,2012), a = runif(32,1,33), b = rnorm(32, 6, 18), c = rnorm(32, 3, 12), 
                 d = rnorm(32, 0, 18))

然后是这个多层光栅：

rs <- stack()

for (i in 1:1:32){
  xy <- matrix(rnorm(400),20,20)
  # Turn the matrix into a raster
  rast <- raster(xy)
  # Give it lat/lon coords for 20-30°E, 43-49°N
  extent(rast) <- c(20,30,43,49)
  rs <- addLayer(rs, rast)
}


# create a Z field for raster just created
years <- seq(as.Date("1981-01-01"), as.Date("2012-12-31"), by = "years")
aa <- setZ(rs, years)
names(rs) <- years

我的问题是：如何获得表示数据帧 df 中的每一列与栅格堆栈 rs 之间的相关性（比如说 Spearman）的五个栅格？

感谢大家的帮助！

【问题讨论】：

标签： r dataframe correlation raster

【解决方案1】：

我不确定你到底想做什么。 df每列有32个值，RasterStack中有32层400个值，

也许您正在寻找df 中的列与层的平均值的相关性？你可以这样做：

你的数据

set.seed(0)
df <- data.frame(year = seq(1981,2012), a=runif(32,1,33), b=rnorm(32, 6, 18), c=rnorm(32, 3, 12), d=rnorm(32, 0, 18))
r <- raster(nrow=20, ncol=20, ext=extent(20,30,43,49))
rs <- stack(lapply(1:32, function(i) setValues(r, rnorm(400,20,20))))
years <- seq(as.Date("1981-01-01"), as.Date("2012-12-31"), by = "years")
names(rs) <- years

解决方案

x <- cellStats(rs, mean)
sapply(2:5, function(i) cor(x, df[,i]))
#[1]  0.123391584 -0.007801092 -0.124336155  0.060774465

【讨论】：

感谢您的提示。实际上，我希望从 df 中获取每列的相关栅格，并使用 rs 中的值。基本上，光栅堆栈中的每个像素与每列的 df 值的相关性，基于年份列。结果应该是一个有 5 层的堆栈（正如您将使用 corLocal 函数获得的那样）。但是，我会尝试在您的示例上构建一些东西。
要计算相关性，您需要两个长度相同的向量。您想将 1 个值与 32 个值进行比较。您可以用它计算一些统计数据，但不能计算相关系数。

【解决方案2】：

好吧，我想出了一个解决方案；不知道是否是最好的，但我认为正在工作。这是df中a列的示例；我为 a 列中的每一行创建了一个虚拟栅格图层；之后，我使用 corLocal 进行关联：

### create a raster layer for each row (year) for column 'a' in df
rs.r <- stack()
library(data.table)

### extract x and y coordinates for raster rs to create a raster stack 
cord <- rasterToPoints(rs[[1]], spatial = F)
cord<- cord[,1:2]
head(cord)

### create a raster where each layer is the value in column a from df
year.s <- unique(df$year)

for (i in 1:length(df$year)){
  print(df$year[i])
  re <- df$a[df$year==year.s[i]]
  c <- data.table(x = cord[,1], y = cord[,2], tt = re)
  m <- rasterFromXYZ(c)
  crs(m) <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 "
  rs.r <- addLayer(rs.r, m)
  crs(rs.r) <-" +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
}

names(rs.r) <- df$year ### set the names for the layers

ext <- extent(rs)
rs.r <- setExtent(rs.r, ext)


rs.r<- projectRaster(rs.r, rs,method = 'ngb')

spplot(corLocal(rs.r, rs, 'spearman'))

【讨论】：