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

Posted 2023-03-23

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

【英文标题】：Calculating area of overlap between 2 raster layers with different origins and extents in R

【发布时间】：2013-07-29 14:45:21

【问题描述】：

我对关于 *** 的问题完全陌生，而且或多或少是 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) 的值的单位是什么？我想单位并不重要，因为我最终会取一个比率（通过将重叠除以更受限制的物种的范围），但我很想知道以供将来参考。

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

【问题讨论】：

我可能误解了确切的问题，所以请对我的评论持保留态度。您能否提供包含您的在线数据的 ASCII 格式示例？我怀疑我们可以更直接地从位置数据计算区域之间的重叠。此外，如果可能，请发布您正在创建的情节样本。一张图讲一千个字……

好吧，显然我错过了 area() 的描述，它告诉你单位是 km^2。

@Dinre - 显然我没有在我的问题中添加图像的声誉。这是一个链接 imgur.com/tHL99j6"><imgsrc="i.imgur.com/tHL99j6.png" title="Hosted by imgur.com" />。在我有足够的资金将澳大利亚添加到image. 至于 ASCII 格式，我不知道怎么描述；它是一系列 2s 和 -3.4e+38 排列的模式，其中 -3.4e+38 代表 NODATA_VALUE 和 2 我假设代表物种的存在.

【参考方案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

【讨论】：

这比我做的更直观，但它创建了一个一般的重叠区域；而我正在研究的其中一个物种在其范围内有中断。我需要它看起来像这样： [IMG]i.imgur.com/tHL99j6.png[/IMG] 另外，如果我要重新投影栅格，我会像你对 CRS() 所做的那样，只指定一个不同的投影吗？我可以用它来计算面积还是仍然是一个近似值？

要重新投影，您需要像 projectRaster 这样的转换函数，第一个 CRS 应用程序只是告诉 R 什么是“从”坐标系，它可以是任何东西，并且需要指定如何到达“到”坐标系“ 坐标系。我确信在这种情况下 area() 的近似值会很好，栅格的重新投影会增加更多问题。您可以使用多边形数据集屏蔽像素，并根据测试对像素区域求和以获得仅陆地部分。

我知道你现在在寻找什么，你想要组合范围的区域。我已经更改了上面的代码来回答您的问题。不过，总的来说，栅格可能不是回答这个问题的最佳工具。我会使用矢量格式（sp 包中的SpatialPolygons 是一个很好的实现）。至于面积：area 函数仅适用于未投影的数据。如果您将其投影到面积投影（是的，使用带有 proj4 规范的CRS），面积只是细胞的分辨率乘以阳性细胞的数量。

@Noah - 这正是我所需要的！谢谢！我也会尝试使用 SpatialPolygons。为什么矢量比光栅更好，只是出于好奇（为了精确）？

@gen 根据定义，矢量形状是封闭区域的规范，这就是您要描述的内容。此外，使用向量可以让您访问向量数学，其中包括非常简单的联合（重叠）方程。