如何在R中高效读取与矢量数据相交的部分大型栅格数据？

问题背景 #

现有投影为A的大型栅格数据集（推荐采用云优化GeoTIFF格式），以及投影为B、仅覆盖栅格一小部分的矢量数据集。需要在R中找到高效获取两者相交部分栅格数据的方法，避免常规流程中重投影整个栅格的高耗时，同时不将完整栅格加载到内存。

常规方法的弊端：
常规流程会先读取完整栅格、重投影整个数据集，再进行裁剪，这种方式对大型栅格来说耗时极长，且会占用大量内存：

ras <- terra::rast(rasterpath)
vec <- sf::st_read(vectorpath)

ras_reproject <- terra::project(ras, terra::crs(vec))
ras_crop <- terra::crop(ras_reproject , vec)

解决方案：先裁剪再重投影 #

核心逻辑是先缩小处理范围，再进行重投影，具体步骤如下：

1. 使用terra包实现（推荐） #

library(terra)
library(sf)

# 读取矢量数据，以及栅格的元数据（不加载完整栅格到内存）
vec <- sf::st_read(vectorpath)
ras <- terra::rast(rasterpath)

# 将矢量数据重投影到栅格的坐标系（A）
vec_proj <- sf::st_transform(vec, terra::crs(ras))

# 仅裁剪栅格中与矢量相交的区域（此时仅加载裁剪区域的栅格数据）
ras_crop <- terra::crop(ras, vec_proj)

# 对裁剪后的小栅格重投影到矢量原坐标系（B）
ras_final <- terra::project(ras_crop, terra::crs(vec))

2. 云优化GeoTIFF的额外优势 #

如果你的栅格是云优化GeoTIFF（COG），terra会自动利用其分块存储特性，仅读取裁剪区域对应的瓦片数据，无需加载整个文件，能进一步提升处理效率。

3. 底层控制方案（gdalUtils） #

如果需要更精细的控制，可以使用gdalUtils调用GDAL工具，直接读取指定区域并完成重投影：

library(gdalUtils)
library(sf)
library(terra)

vec <- sf::st_read(vectorpath)
ras <- terra::rast(rasterpath)

# 将矢量重投影到栅格坐标系，获取其边界范围
vec_proj <- sf::st_transform(vec, terra::crs(ras))
ext <- sf::st_bbox(vec_proj)

# 调用GDAL直接读取指定区域并投影
gdal_translate(
  src_dataset = rasterpath,
  dst_dataset = "cropped_proj.tif",
  projwin = c(ext$xmin, ext$ymax, ext$xmax, ext$ymin),
  projwin_srs = terra::crs(ras),
  t_srs = terra::crs(vec)
)

# 读取处理后的结果
ras_final <- terra::rast("cropped_proj.tif")

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