在R中计算两个栅格图层之间重叠区域的面积，这两个图层具有不同的起始点和范围。

Question

在R中计算两个栅格图层之间重叠区域的面积，这两个图层具有不同的起始点和范围。

4

我完全是在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的重叠部分的栅格（当我试图在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)的值的单位是什么？我想这些单位并不太重要，因为最终我会通过将重叠部分除以更严格物种的范围来取比率，但我很想知道以备将来参考。

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

- g.m

我可能误解了确切的问题，所以请谨慎考虑我的评论。你能否提供一份包含你出席数据的ASCII格式样本？我怀疑我们可以更直接地通过位置数据计算区域之间的重叠情况。此外，如果可能，请发布你正在创建的图表样本。图片胜过千言万语... - Dinre

好的，显然我在area()的描述中错过了一点，它告诉你单位是km^2。 - g.m

@Dinre - 显然我没有足够的声望来添加图片到我的问题。这里是一个链接 <a href="http://imgur.com/tHL99j6"><img src="http://i.imgur.com/tHL99j6.png" title="Hosted by imgur.com" /></a>。在我还没来得及添加澳大利亚到图像之前，R已经崩溃了。至于ASCII格式，我不知道如何描述它；它是一系列的2和-3.4e+38按照一定的模式排列，其中-3.4e+38代表NODATA_VALUE，而2则表示物种的存在。 - g.m

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Noah · Accepted Answer

使用 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")

enter image description here

area函数可获得每个单元格的近似面积（为获取真实面积，应将其重新投影到面积坐标系统）。要获取重叠区域的总近似面积，请将所有单元格值相加，并通过组合栅格的值对面积求和进行索引：

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