Jaccard 索引：R 不让我将光栅转换为数字？

Question

我正在尝试为我的数据计算 Jaccard 索引，该索引用于两个存在/不存在栅格（两个不同物种）。我还有第三个栅格，其中包含物种范围的重叠。

我想计算Jaccard指数，试过下面的代码

g1_binary <- as.integer(pg1_rast > tr1)
pg3_binary <- as.integer(pg3_rast > tr3)

# Find the overlap between the two rasters
overlap_binary <- pg1_binary * pg3_binary

# Convert the overlap raster to a RasterLayer object
overlap_rast <- raster(overlap_binary)
extent(overlap_rast) <- extent(pg1_rast)

# Calculate area of overlap, area of pg1, and area of pg2
a <- sum(as.matrix(overlap_rast)) * res(overlap_rast)[1] * res(overlap_rast)[2]
b <- sum(as.matrix(pg1_binary)) * res(pg1_rast)[1] * res(pg1_rast)[2]
c <- sum(as.matrix(pg3_binary)) * res(pg3_rast)[1] * res(pg3_rast)[2]

# Calculate degree of overlap
degree_overlap <- as.numeric(a / (b + c))

但是，我不断收到错误消息：

Error in as.numeric(a/(b + c)) : 
  cannot coerce type 'S4' to vector of type 'double'

我已经尝试了一百万种不同的方法，但似乎无法将我的数据转换为 Jaccard 索引的数字。

我也尝试将其计算为：

overlap_rast/(pg1_rast+pg3_rast)

虽然我什至不知道这是否可行，但它仍然给我同样的错误，即我的栅格无法转换为数值。请帮忙！

我试过上面的方法，但它给了我显示的错误。

Answer 1

示例数据（通过在代码中生成或使用 R 附带的文件，在询问 R 问题时始终包含一些数据。有关示例，请参见 R 帮助文件）

library(terra)
set.seed(1)
r1 <- rast(nrow=10, ncol=10, vals=runif(100))
r2 <- rast(nrow=10, ncol=10, vals=runif(100))

解决方案

b1 <- r1 > 0.5
b2 <- r2 > 0.5
 
b12 <- b1 * b2

x <- c(b1, b2, b12)
g <- global(x, sum, na.rm=TRUE)[[1]]
J <- g[3] / (sum(g[1:2]) - g[3])
J
#[1] 0.3246753

这是假设所有单元格都具有相同的大小，就像您在脚本中所做的那样。在这种情况下，无需将

a

、

b

和

c

与常数单元格大小相乘。这对于较小的区域来说很好，但在其他情况下可能是错误的。

要考虑可变的单元格大小，您可以这样做

a <- cellSize(b1, mask=FALSE, unit="km")
xa <- c(b1, b2, b12) * a
ga <- global(xa, sum, na.rm=TRUE)[[1]]
Ja <- ga[3] / (sum(ga[1:2]) - ga[3])
Ja
#[1] 0.3352633

另请注意，您有

J = a / (b + c)

但 Jaccard 指数可以计算为

J = a / (b + c - a)

Answer 2

谢谢大家！实际上，我从此处的另一个非常有用的页面上发现，您可以通过首先转换为矩阵，然后转换为数字来更改它，例如使用 as.matrix() 然后使用 as.numeric() 来总结所有内容。以防万一对任何人有帮助！

Answer 3

Jaccard指数计算需要数值，不能使用'as.numeric()'函数将'S4'对象强制转换为数值。

要修复此错误，您可以使用“raster”包中的“getValues()”函数从“S4”对象中提取数值。试试这个代码：

library(raster)

g1_binary <- as.integer(pg1_rast > tr1)
pg3_binary <- as.integer(pg3_rast > tr3)

# Find the overlap between the two rasters
overlap_binary <- pg1_binary * pg3_binary

# Convert the overlap raster to a RasterLayer object
overlap_rast <- raster(overlap_binary)
extent(overlap_rast) <- extent(pg1_rast)

# Get the values from the rasters
overlap_values <- getValues(overlap_rast)
pg1_values <- getValues(pg1_rast)
pg3_values <- getValues(pg3_rast)

# Calculate area of overlap, area of pg1, and area of pg2
a <- sum(overlap_values) * res(overlap_rast)[1] * res(overlap_rast)[2]
b <- sum(pg1_values) * res(pg1_rast)[1] * res(pg1_rast)[2]
c <- sum(pg3_values) * res(pg3_rast)[1] * res(pg3_rast)[2]

# Calculate degree of overlap
degree_overlap <- as.numeric(a / (b + c)) # Mind the formula!

注：杰卡德指数是两个集合的交集大小与集合并集大小的比值。在这种情况下，并集的大小是两个集合的大小之和减去交集的大小，因此公式 J = a / (b + c - a)。