主页 / coding / 问题 / 79353000
Accepted
M. Beausoleil
Asked: 2025-01-14 01:28:18 +0800 CST

使用 terra 来计算每个像素的面积，正确的方法是什么？

  • 772

我有一个栅格，我想将其值更改为每公顷像素值，其中面积来自像素的大小。因此，如果一个像素是 10 公顷，其当前值为 100 个任意单位，那么现在它将是 10 个任意单位/公顷。

栅格示例：

r <- rast(nrows=18, ncols=36)
v <- 1:ncell(r)
v[200:400] <- NA
values(r) <- v

此函数应采用栅格并乘以 1/ha。

px2ha <- function(raster) {
  raster * 10000/prod(res(raster))
}

那么像这样使用 cellSize 怎么样？

r2/cellSize(r2, unit = 'ha')

我没有得到相同的答案。

r2 = project(x = r, y = 'epsg:2950') # Put raster in meters 
plot((r2))
plot(px2ha(r2))
plot(r2/cellSize(r2, unit = 'ha'))

可以使用 cellSize 获取每面积（公顷）的像素值吗？

1 个回答

  1. Best Answer

    正确的做法是使用terra::cellSize。如果您的栅格数据是经度/纬度，则应使用它。

    library(terra)
r <- rast(nrows=18, ncols=36)
a <- cellSize(r, unit="ha")
minmax(a)
#         area
#min  10803350
#max 122787719

    注意值的范围。像元大小会发生变化。投影后也是如此。数字不一样是因为栅格会随着投影而发生根本性变化（这就是为什么您通常要避免该步骤）。

    r2 <- project(r, "epsg:2950")
b <- cellSize(r2, unit="ha")
minmax(b)
#         area
#min   3535976
#max 161759413

    除非你使用transform=FALSE

    d <- cellSize(r2, unit="ha", transform=FALSE)
minmax(d)
#         area
#min 162812887
#max 162812887

    这与（请注意，你的公式是错误的）相同

    prod(res(r2)) / 10000
#[1] 162812887

    但使用此数字并不正确，因为它假设没有失真。此数字在坐标 (0,0) 附近应该是正确的，但距离该位置越远，失真就越大。cellSize通过在计算面积之前将单元格转换回经度/纬度来纠正此问题。

    在这个例子中，失真很大，但这有点不现实。这是因为对加拿大部分地区的投影被应用于整个世界。假设您的实际数据仅在魁北克和/或安大略省东部，那里的失真会较小。

    较旧的教科书和课程可能会建议您将经度/纬度数据投影到平面坐标参考系统来计算面积，但这反映了所用软件的局限性。这不是一个好方法，事实上应该采取相反的做法。

    • 3

相关问题