为什么terra::cellSize()和raster::area()会给出不同的栅格单元面积估计？

Question

为什么terra::cellSize()和raster::area()会给出不同的栅格单元面积估计？

7

我刚刚注意到 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包（v0.3.0）于2021年7月8日创建}

过去几年中，我一直在使用raster库（并且非常喜欢它），最近被新的terra包提供的速度和其他改进所震撼。我假设terra::cellSize()（或者说terra::expanse()）提供的面积估计比raster::area()更精确，但在更新以前的面积估计之前，我想更多地了解发生了什么变化。

感谢你所做的一切，https://github.com/rspatial！

- Christopher L. Crawford

1

感谢您在这里发布。GitHub 是解决 bug 的最佳平台，Stack Overflow 是解决编程问题的最佳平台。虽然这不属于以上两者，但它可能是在这里最有用的。 - Robert Hijmans

太好了 - 这正是我所希望的。再次感谢你的回答！ - Christopher L. Crawford

1个回答

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

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

- Robert Hijmans · Accepted Answer

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 还计算了投影（即非经纬度）栅格的真实单元格大小；而不是假设名义分辨率是恒定的。

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))