如何用{stars}包基于vsicurl远程计算Sentinel-2数据的NDVI？

核心思路 #

利用stars结合GDAL的/vsicurl/特性直接读取云端COG数据，无需本地存储；同时借助sf的CRS工具链统一语法，自动触发CRS不匹配警告，规避教学中常见的坐标错误。

完整代码实现 #

# 加载依赖包
library(rstac)
library(stars)
library(sf)

# 1. 定义研究边界（用sf语法统一处理CRS）
bbox_sf <- st_bbox(c(xmin = -122.4, ymin = 37.7, xmax = -122.3, ymax = 37.8), 
                   crs = st_crs(4326)) %>% 
  st_as_sfc()

# 2. 通过rstac获取Sentinel-2影像元数据
stac_result <- stac("https://earth-search.aws.element84.com/v0") %>%
  stac_search(collections = "sentinel-s2-l2a-cogs",
              bbox = st_bbox(bbox_sf),
              datetime = "2023-06-01/2023-06-30",
              limit = 1) %>%
  get_request()

# 提取红波段(B04)和近红外波段(B08)的远程URL，添加vsicurl前缀
assets <- stac_result$features[[1]]$assets
b04_url <- paste0("/vsicurl/", assets$B04$href)
b08_url <- paste0("/vsicurl/", assets$B08$href)

# 3. 远程读取并裁剪数据（自动检查CRS，不匹配时抛出警告）
crs_img <- st_crs(b04_url)
bbox_cropped <- st_bbox(st_transform(bbox_sf, crs_img))

b04 <- read_stars(b04_url, RasterIO = list(bbox = bbox_cropped))
b08 <- read_stars(b08_url, RasterIO = list(bbox = bbox_cropped))

# 4. 内存中直接计算NDVI
ndvi <- (b08 - b04) / (b08 + b04)

# 可选：转换为sf兼容格式用于可视化/分析
ndvi_sf <- st_as_sf(ndvi, as_points = FALSE)

关键优化点 #

• 无本地存储：通过/vsicurl/前缀直接读取云端数据，跳过下载步骤，节省磁盘空间与时间
• CRS安全处理：用sf的st_transform统一转换边界CRS，stars会自动校验CRS匹配度并抛出警告，帮助学生及时发现坐标错误
• 高效裁剪：在read_stars阶段通过RasterIO$bbox直接裁剪，避免读取全量影像的内存浪费
• 简洁运算：利用stars的矢量运算特性，直接对波段对象做算术操作，语法与sf一致

教学优势对比{terra} #

• 强制CRS匹配校验：避免学生因忽略坐标系统差异导致的分析偏差
• 语法统一：全程复用sf的坐标处理逻辑，降低跨包学习成本
• 内存效率相当：依托GDAL底层优化，远程读取与计算的内存消耗与{terra}持平

内容的提问来源于stack exchange，提问作者cboettig