r - 使用 {stars} 包对 Sentinel-2 数据进行远程 (vsicurl) 计算？

Question

我正在使用 rstac 访问所需边界框和日期范围内的 Sentinel-2 数据并计算 NDVI。当我使用 {terra} 包时，这对我来说相对[*]干净且直接，但我想使用 {stars} 语法(更多关于为什么更进一步)。

首先，通过快速 {rstac} 查询来获取数据的 URL:

library(rstac)
library(sf)
library(stars)
library(terra) 

bbox <- st_bbox(c(xmin=-86.94663, ymin=33.43930, 
                  xmax=-86.67684, ymax=33.62239),
                crs=4326) # Birmingham, AL
matches <- 
  stac("https://planetarycomputer.microsoft.com/api/stac/v1/") |>
  stac_search(collections = "sentinel-2-l2a",
              bbox = bbox,
              datetime = "2019-06-01/2019-08-01") |>
  get_request() |>
  items_sign(sign_fn = sign_planetary_computer())

这会返回大量匹配项和适当的元数据，为简单起见，我只从属性元数据中选择一个 eo:cloudcover 较低的 (#59):

best_img <- matches$features[[59]] 

现在我将使用 vsicurl无需下载整个文件即可访问红色和近红外波段的机制。这些图像比我的搜索框大得多，因此我还想裁剪掉那些我不会使用的像素，以避免不必要的计算。

我的第一步很丑陋。要使用 {terra} 裁剪我的图像，我需要一个 SpatVec cookie-cutter 来传递到 crop() 。我已经有bbox上面作为 sf 类型边界框，我执行以下操作以将其放入与 Sentinel2 资源匹配的投影中，但这感觉非常 hacky。我想要一个简洁的、纯粹的版本，但这可行:

red <- read_stars( paste0("/vsicurl/", best_img$assets$B04$href) )
bbox_proj <- bbox |> st_as_sfc() |> st_transform(st_crs(red)) |> vect()

无论如何，在手中裁剪矢量，terra 中的 NDVI 计算非常优雅且快速(在使用最少 RAM 的良好网络连接上):

red <- rast( paste0("/vsicurl/", best_img$assets$B04$href) ) |> crop(bbox_proj)
nir <- rast( paste0("/vsicurl/", best_img$assets$B08$href) ) |> crop(bbox_proj)
ndvi_fun <- function(x, y) (x - y) / (x + y)
ndvi <- lapp(c(nir, red), fun = ndvi_fun)
ndvi  |> plot()

所以我的主要问题是使用 {stars} 进行相同计算的等效语法是什么？到目前为止，我只提出了下面的代码，特别是仅在使用 local 时才有效。我必须首先创建副本(因此毫不奇怪，速度也慢得多!)

# ugh why can't we combine these in a single read_stars?
red <- read_stars( paste0("/vsicurl/", best_img$assets$B04$href) )
nir <- read_stars( paste0("/vsicurl/", best_img$assets$B08$href) )

bbox_proj <- bbox |> st_as_sfc() |> st_transform(st_crs(red))

# combine 'by hand' and then crop... 
remote <- c(r1,r2, along=3) |> st_crop(bbox_proj)

# ugh! ugh! why do we have to use local copy for this to work??
stars::write_stars(remote, "test.tif")
local <- read_stars("test.tif")

# Um, I think this is correct NDVI method, hope I didn't reverse the bands...
# also surprised this is considerably slower and uses much more RAM
calc_ndvi <- function(x) (x[2] - x[1])/(x[2] + x[1])
ndvi <- st_apply(local, 1:2,  FUN = calc_ndvi)
plot(ndvi, col =  rgb(0, (0:100)/100, 0))

我肯定在我的 star 语法中遗漏了一些东西，这导致它变得更慢，表达起来更冗长，并且只有在 st_apply() 时才起作用。运行于 local复制而不是 remote对象。

风格和偏好

如果它在 {terra} 中工作，那么问为什么在 {stars} 中这样做也许是合理的——部分原因是我学习 star，但我也是一名讲师，总是觉得教我的学生两个 sf 很麻烦和 terra 语法。如果您在不重新投影边界框 CRS 的情况下尝试上述裁剪，terra 也不会警告错误匹配的 CRS，这是学生常见的错误。在这两种情况下，我发现裁剪边界框的重新投影也比我喜欢的更麻烦。特别是访问文件“两次”似乎很尴尬，一次读取 crs，然后再次裁剪，我希望可以有更优雅的语法，但还没有弄清楚。

Answer 1

这并不能回答你的问题，但这里有一个更简洁的 terra 方法，使用 project<SpatExtent>我在terra 1.6.31(目前是开发版本)中引入的方法。

library(rstac)
library(terra) 
#terra 1.6.31

bbox <- c(xmin=-86.94663, ymin=33.43930, xmax=-86.67684, ymax=33.62239)

matches <- stac("https://planetarycomputer.microsoft.com/api/stac/v1/") |>
        stac_search(collections = "sentinel-2-l2a",
            bbox = bbox, datetime = "2019-06-01/2019-08-01") |>
        get_request() |>  items_sign(sign_fn = sign_planetary_computer())

best_img <- matches$features[[59]] 

获取第一个数据源，并将经/纬度搜索范围投影到数据的坐标引用系。请注意(新)参数 xy=TRUE 的使用指示 bbox 数字的顺序为 (xmin,ymin,xmax,ymax )，而“terra”默认情况下需要 (xmin,xmax,ymin,ymax)。

red <- rast( paste0("/vsicurl/", best_img$assets$B04$href) ) 
e <- ext(bbox, xy=TRUE) |> 
     project("+proj=longlat", crs(red))

裁剪第一个数据源并下载并裁剪其他数据源

red <- crop(red, e)
nir <- rast( paste0("/vsicurl/", best_img$assets$B08$href) ) |> crop(e)

并使用数据

ndvi_fun <- function(x, y) (x - y) / (x + y)
ndvi <- lapp(c(nir, red), fun = ndvi_fun)

以上使用lapp很棒，但是，对于这样的简单函数，下面的速度更快

ndvi <- (red-nir) / (red+nir)

如果您要使用lapp ，你可以考虑这样做

rednir <- paste0("/vsicurl/", c(best_img$assets$B04$href, best_img$assets$B08$href)) |> 
          rast() |> crop(e, names=c("red", "nir"))    
ndvi <- lapp(rednir, ndvi_fun)