为什么terra::cellSize()和raster::area()会生成不同的栅格单元面积估计值? [英] Why do terra::cellSize() and raster::area() produce different estimates of raster cell area?
本文介绍了为什么terra::cellSize()和raster::area()会生成不同的栅格单元面积估计值?的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!
问题描述
我刚刚注意到terra::cellSize()
生成的单元格面积估计值与raster::area()
生成的单元格面积估计值不匹配。
library(raster)
#> Loading required package: sp
library(terra)
#> terra version 1.3.4
# make test raster with raster::raster()
a <- raster::raster(ncols = 100, nrows = 100,
xmn = -84, xmx = -83,
ymn = 42, ymx = 43)
# make test raster with terra::rast()
b <- terra::rast(ncols = 100, nrows = 100,
xmin = -84, xmax = -83,
ymin = 42, ymax = 43)
# calculate cell areas (km2)
a_area <- raster::area(a) # km by default
b_area <- terra::cellSize(b, unit = "km")
# sum across cells
a_sum <- raster::cellStats(a_area, "sum")
b_sum <- terra::global(b_area, fun = "sum")
a_sum
#> [1] 9088.98
b_sum
#> sum
#> area 9130.795
# note that this terra workflow yields the same answer as terra::expanse()
terra::expanse(b, unit = "km")
#> [1] 9130.795
sessionInfo()
#> R version 4.0.2 (2020-06-22)
#> Platform: x86_64-apple-darwin17.0 (64-bit)
#> Running under: macOS 10.16
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
#>
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] terra_1.3-4 raster_3.4-5 sp_1.4-5
#>
#> loaded via a namespace (and not attached):
#> [1] Rcpp_1.0.6 codetools_0.2-18 lattice_0.20-41 digest_0.6.27
#> [5] grid_4.0.2 magrittr_2.0.1 evaluate_0.14 highr_0.8
#> [9] rlang_0.4.10 stringi_1.5.3 rmarkdown_2.6 rgdal_1.5-19
#> [13] tools_4.0.2 stringr_1.4.0 xfun_0.20 yaml_2.2.1
#> [17] compiler_4.0.2 htmltools_0.5.1.1 knitr_1.31
packageVersion("raster")
#> [1] '3.4.5'
packageVersion("terra")
#> [1] '1.3.4'
由reprex package(v0.3.0)于2021-07-08创建
在过去的几年里,我一直在使用raster
(并且非常喜欢它),最近我被较新的terra
包提供的速度和其他改进所震撼。我假设terra::cellSize()
(或terra::expanse()
)提供的面积估计比raster::area()
更准确,但在更新以前的面积估计之前,我想了解更多有关变化的信息。
感谢您所做的一切,https://github.com/rspatial!
推荐答案
raster
使用单元格的宽度(经度)和高度(纬度)的乘积。terra
更精确,它计算单元格的球形面积(由其四角定义)。这一点在高纬度地区最为重要,因为那里的单元格宽度变化最大,并且单元格的底部和顶部可能会有所不同。因此,高纬度地区和垂直分辨率较低的单元格之间的差异将最大。
插图如下:
library(raster)
library(terra)
a <- raster::raster()
b <- terra::rast()
values(a) <- 1:ncell(a)
values(b) <- 1:ncell(b)
a_rast <- rast(area(a))
a_terra <- cellSize(b, unit = "km")
dif <- a_terra - a_rast
reldif <- 100 * dif / a_terra
plot(reldif)
在两极附近,这些细胞的差异几乎是1%,而在赤道是~0。
terra
中更大的更改是,它还计算投影(即不是纵向)栅格的实际像元大小。来自?cell Size的示例
r <- rast(ncols=90, nrows=45, ymin=-80, ymax=80)
m <- project(r, "+proj=merc")
bad <- init(m, prod(res(m)) / 1000000, names="naive")
good <- cellSize(m, unit="km", names="corrected")
plot(c(good, bad), nc=1, mar=c(2,2,1,6))
不幸的是,这在当时似乎有点问题,特别是对于大型栅格;希望在下一个版本中可以解决这一问题。
这篇关于为什么terra::cellSize()和raster::area()会生成不同的栅格单元面积估计值?的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!
查看全文