根据R中2个栅格之间的土地覆盖类型计算R2值？

Question

我在 r 中有三个栅格

> lpjre
class      : RasterLayer 
dimensions : 2803, 5303, 14864309  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : 60.85, 105.0417, 15.95833, 39.31667  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : xxx 
values     : 0, 21  (min, max)

> gcre
class      : RasterLayer 
dimensions : 2803, 5303, 14864309  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : 60.85, 105.0417, 15.95833, 39.31667  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : 0, 39.72  (min, max)

和土地覆盖栅格

> tif4
class      : RasterLayer 
dimensions : 2803, 5303, 14864309  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : 60.85, 105.0417, 15.95833, 39.31667  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : C:/Users/XXXX/landusemaskedme.tif 
names      : landusemaskedme 
values     : 1, 12  (min, max)
attributes :
       ID                          zn
 from:  1 evergreen needleleaf forest
  to : 12                   croplands

我根据不同的土地覆盖等级在lpjre和gcre之间绘制散点图

plot(lpjre[tif4==2],gcre[tif4==2]) 

我想知道如何根据土地覆盖类型计算lpjre 和gcre 之间的r2 值？

我尝试这段代码并给出错误：

> cor(values(lpjre)[tif4==1], values(gcre)[tif4==1], use="complete.obs", method = 'pearson')
Error in values(gcre)[tif4 == 1] : invalid subscript type 'S4'

可重现的栅格：

library(raster)
ras1 <- raster(matrix(c(1,1,1,2,2,2)))
ras2 <- raster(matrix(c(1,1,1,2,2,2)))

#Generating landcover example data
raster2 <- raster(matrix(c(1,1,1,2,2,2,3,3,3),ncol =3))
raster2 <- as.factor(raster2)

rat <- levels(raster2 )[[1]]
rat[["landcover"]] <- c("land","ocean/lake", "rivers")
levels(raster2 ) <- rat

Answer 1

您可以为此使用矢量转换。首先是一个可重现的例子：

library(raster)
set.seed(42)

ras1 <- raster(nrow = 10, ncol = 10)
ras2 <- raster(nrow = 10, ncol = 10)
lcc <- raster(nrow = 10, ncol = 10)
ras1[] <- runif(100)
ras2[] <- runif(100)
lcc[] <- sample(c(1,2), 100, replace = TRUE)

现在您可以使用[] 将栅格值转换为矢量。从那里开始，您甚至可以在栅格之间使用逻辑选择，假设它们具有相应的几何特征。

cor(ras1[lcc[] == 2], ras2[lcc[] == 2], use = "complete.obs", method = "pearson")
# [1] -0.1644459