将砖块中的每个栅格层与数据帧 R 中的值相关联
Correlate each raster layer in a brick with values in a dataframe R
我在 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 和栅格堆栈 中每一列之间的相关性(比方说 Spearman)的五个栅格]rs?
谢谢大家的帮助!
我不确定你到底想做什么。 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 中 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'))
我在 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 和栅格堆栈 中每一列之间的相关性(比方说 Spearman)的五个栅格]rs?
谢谢大家的帮助!
我不确定你到底想做什么。 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 中 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'))