0

第一篇文章:)

我一直在将我的 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 中的每个值除以其网格区域，它没有按预期工作（也许是因为我的对象是星星而不是科幻？）

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

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

  1. 我担心当我将 x_area$area（一个矩阵，纬度/经度的区域）转换为数字向量时，我可能会弄乱网格与其区域之间的纬度/经度匹配。我做了一些粗略的测试，看看这些区域是否符合我的预期，但我无法避免担心这可能会导致难以捕捉的错误。

  2. 我以正确的单位从“x_area”开始，只是在计算过程中删除然后再次设置单位，这似乎并不干净。

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

谢谢！

4

1 回答 1

0

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

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))
于 2021-02-24T02:56:04.137 回答