r - 计算R中具有不同原点和范围的2个栅格层之间的重叠区域

Question

我完全不熟悉关于 stackoverflow 的问题，并且或多或少是 R(以及一般编程)的新手，所以请耐心等待。

我有仅显示存在的物种范围的 ASCII 文件。在搜索了互联网的各个角落后，我设法上传、转换为栅格、沿所需边界(在我的情况下是澳大利亚的海岸线)进行 mask ，并绘制它们，以便我可以在未投影的 map 上可视化范围。

完成了定性方面的工作后，我需要进入定量方面；也就是说，我需要计算物种之间的共鸣程度。为了做到这一点，我需要首先计算重叠区域，这是我遇到问题的地方。到目前为止，这是我设法做的事情:

> 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