按网格区域归一化 R stars 对象？

Question

第一个post:)

我一直在将我的 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（矩阵，lat/lon 中的区域）转换为数字向量时，我可能会弄乱 lat/lon 网格与其区域之间的匹配。我做了一些粗略的测试，看看这些区域是否符合我的预期，但无法避免这可能导致难以捕捉的错误的担忧。
我以正确的单位从“x_area”开始似乎不太干净，只是在计算过程中删除然后重新设置单位。

有人可以为我正在尝试做的事情建议一个“更干净”的实现，即在整个过程中保持单元的同时将网格乘以或除以其面积吗？或者让我相信我的代码没问题？

谢谢！

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

按网格区域归一化 R stars 对象？

Normalizing an R stars object by grid area?

r

spatial

r-stars