按网格区域规范化 R 星对象？答案

【问题标题】：Normalizing an R stars object by grid area?按网格区域规范化 R 星对象？
【发布时间】：2021-05-26 02:28:33
【问题描述】：

第一篇文章 :)

我一直在将我的 R 代码从 sp() 转换为 sf()/stars()，而我仍在努力掌握的一件事是计算网格中的区域。

这里有一个示例代码来解释我的意思。

library(stars)
library(tidyverse)

# Reading in an example tif file, from stars() vignette
tif = system.file("tif/L7_ETMs.tif", package = "stars")
x = read_stars(tif)
x

# Get areas for each grid of the x object. Returns stars object with "area" in units of [m^2]
x_area <- st_area(x)
x_area

我尝试松散地采用此小插图 (https://github.com/r-spatial/stars/blob/master/vignettes/stars5.Rmd) 中的代码来将 x 中的每个值除以其网格区域，但它没有按预期工作（可能是因为我的对象是星星而不是 sf？）

x$test1 = x$L7_ETMs.tif / x_area  # Some computationally intensive calculation seems to happen, but doesn't produce the results I expect?

x$test1 = x$L7_ETMs.tif / x_area$area # Throws error, "non-conformable arrays"

以下似乎有效。

x %>%
  mutate(test1 = L7_ETMs.tif / units::set_units(as.numeric(x_area$area), m^2))

以下是我对这段代码的担忧。

我担心当我将 x_area$area（一个矩阵，纬度/经度的区域）转换为数字向量时，我可能会弄乱网格与其区域之间的纬度/经度匹配。我做了一些粗略的测试，看看这些区域是否符合我的预期，但我无法避免担心这可能会导致难以捕捉的错误。
我从“x_area”开始以正确的单位开始似乎并不干净，只是在计算过程中删除然后再次设置单位。

有人可以为我正在尝试做的事情提出一个“更清洁”的实施方案，即在始终保持单位的同时将网格乘以或除以其面积？或者说服我我的代码没问题？

谢谢！

【问题讨论】：

标签： r spatial r-stars

【解决方案1】：

我不知道如何改进stars代码，但是您可以将得到的结果与此进行比较

tif <- system.file("tif/L7_ETMs.tif", package = "stars")
library(terra)
r <- rast(tif)

a <- cellSize(r, sum=FALSE)
x <- r / a

对于平面数据，您可以在可以安全地假设没有失真的情况下执行此操作（通常情况并非如此，但它可能是这种情况）

y <- r / prod(res(r))

【讨论】：

感谢您的 cmets！我帖子中代码的“平面单元”方面是我（懒惰）尝试创建可重现示例的意外结果，我的真实数据是经度/纬度，网格区域有变化。另外，感谢您指出“terra”包（您的？），我正在阅读文档并发现它非常有趣（包括可能是关于栅格与空间向量有何不同的最好、最直接的解释，我读过！），一定会探索我能用它做什么！
我发现terra::cellSize() 是我要交叉检查我的st_area() 计算的函数。另请注意，st_area() 要求 sf/stars 对象具有 CRS 以返回 m^2 单位的区域。我原始帖子中的“不合格数组”错误实际上是由于矩阵维度不匹配问题，因为这些区域只有 2 个维度（x，y），而我上面示例中的 x 有第三个维度（@987654327 @)。