计算R中具有不同来源和范围的2个栅格图层之间的重叠区域

Question

我对关于 stackoverflow 的问题完全陌生，而且或多或少是 R 的新手（以及一般的编程），所以请多多包涵。

我有仅显示存在的物种范围的 ASCII 文件。在搜索了互联网的远端之后，我设法上传、转换为栅格、沿着所需的边界（在我的例子中是澳大利亚的海岸线）进行遮罩，并绘制它们，以便我可以在未投影的地图上可视化范围。

完成了这方面的定性方面，我需要进入定量方面；也就是说，我需要计算物种之间的共感程度。为此，我需要首先计算重叠区域，这是我遇到问题的地方。以下是我到目前为止所做的：

> d
class       : RasterLayer 
dimensions  : 85, 270, 22950  (nrow, ncol, ncell)
resolution  : 0.08, 0.08  (x, y)
extent      : 119.4993, 141.0993, -36.65831, -29.85831  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : 2, 2  (min, max)

> b
class       : RasterLayer 
dimensions  : 140, 222, 31080  (nrow, ncol, ncell)
resolution  : 0.08, 0.08  (x, y)
extent      : 134.2456, 152.0056, -40.44268, -29.24268  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : 2, 2  (min, max)

x<-resample(b,d,method="ngb")
y<-mask(x,d)

>y
class       : RasterLayer 
dimensions  : 85, 270, 22950  (nrow, ncol, ncell)
resolution  : 0.08, 0.08  (x, y)
extent      : 119.4993, 141.0993, -36.65831, -29.85831  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : 2, 2  (min, max)

y 是 d 和 b 之间在 d 上被屏蔽的重叠的栅格（当我尝试在 b 上屏蔽时，我收到错误说范围不同）。我如何计算它的面积？ raster 包中的 area() 函数将其吐出：

area(y)
class       : RasterLayer 
dimensions  : 85, 270, 22950  (nrow, ncol, ncell)
resolution  : 0.08, 0.08  (x, y)
extent      : 119.4993, 141.0993, -36.65831, -29.85831  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : 63.65553, 68.75387  (min, max)

我不完全确定该怎么做。这甚至是获得区域的好/准确/正确的方法吗？为什么 y 和 b 的范围不同，而 d 和 y 的范围相同？此外，area(y) 的值的单位是什么？我想单位并不重要，因为我最终会取一个比率（通过将重叠除以更受限制的物种的范围），但我很想知道以供将来参考。

如果这是一个愚蠢的问题，我很抱歉。我很感激有人可能提出的任何意见。

Answer 1

获得重叠的最佳方法是使用intersect。您可以创建重叠值的砖块并使用类似any 的命令来获取重叠范围，假设每个范围内的值为 1 或 TRUE，范围外的值为 0、FALSE 或 NA：

library(raster)

wgs84 <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
d <- raster()
extent(d) <- extent(119.4993, 141.0993, -36.65831, -29.85831)
res(d) <- c(0.08, 0.08)
projection(d) <- CRS(wgs84)
values(d) <- sample(c(NA, 1), ncell(d), replace=TRUE)

b <- raster()
extent(b) <- c(134.2456, 152.0056, -40.44268, -29.24268 )
res(b) <- c(0.08, 0.08)
projection(b) <- CRS(wgs84)
values(b) <- sample(c(NA, 1), ncell(b), replace=TRUE)

y <- intersect(b, d)

x <- brick(resample(b, y, method = "ngb"),resample(d, y, method = "ngb"))
x2 <- any(x, na.rm = TRUE)

library(maps)
map(regions = "australia")
image(d, add = TRUE, col = "blue")
image(b, add = TRUE, col = "green")
plot(extent(y), add = TRUE)
image(x2, add = TRUE, col = "red")

area 函数为您获取每个 单元格的近似 区域（为了获得真实区域，您应该将其重新投影到区域坐标系）。要获得重叠的总近似面积，请将所有像元值相加，按组合栅格的值索引面积总和：

sum(values(area(x2))[which(values(x2))])
# [1] 361407.1