使用Terra通过WCS读取大栅格文件时遇到的问题

在R语言中使用terra包通过WCS服务访问栅格数据时，出现了和覆盖区域大小相关的异常：小范围的WCS请求可通过rast()正常读取，但大范围请求会报错；而stars包的read_stars()对两种范围的请求均能正常处理。此前同类代码运行正常，暂未明确具体原因。

复现代码 #

library(terra)
library(stars)

# 小范围WCS请求链接
smaller_dim <- 'https://www.mrlc.gov/geoserver/rcmap_anhb/wcs?request=GetCoverage&service=WCS&version=2.0.1&coverageid=rcmap_anhb__rcmap_annual_herbaceous_2020&subset=X(-1026935.24156938,-1015591.34251064)&subset=Y(2005364.48088881,2014307.50831563)'
# 大范围WCS请求链接
bigger_dim <- 'https://www.mrlc.gov/geoserver/rcmap_anhb/wcs?request=GetCoverage&service=WCS&version=2.0.1&coverageid=rcmap_anhb__rcmap_annual_herbaceous_2020&subset=X(-1054529.13659319,-1022600.05347719)&subset=Y(2020307.11783295,2056926.21900819)'

# 小范围请求：rast与read_stars均正常运行
x <- rast(smaller_dim)
y <- read_stars(smaller_dim)

# 大范围请求：read_stars正常，rast报错
xx <- rast(bigger_dim)
# 错误信息: [rast] 文件不存在: /vsicurl/https://www.mrlc.gov/geoserver/rcmap_anhb/wcs?request=GetCoverage&service=WCS&version=2.0.1&coverageid=rcmap_anhb__rcmap_annual_herbaceous_2020&subset=X(-1054529.13659319,-1022600.05347719)&subset=Y(2020307.11783295,2056926.21900819)
# 附加警告信息:
# `/vsicurl/https://www.mrlc.gov/geoserver/rcmap_anhb/wcs?request=GetCoverage&service=WCS&version=2.0.1&coverageid=rcmap_anhb__rcmap_annual_herbaceous_2020&subset=X(-1054529.13659319,-1022600.05347719)&subset=Y(2020307.11783295,2056926.21900819)' 未被识别为支持的文件格式。(GDAL错误4) 

yy <- read_stars(bigger_dim)

解决建议 #

1. 调整GDAL环境参数
   terra依赖GDAL处理远程数据，大范围请求可能触发vsicurl的内存或超时限制。可通过设置环境变量调整：

# 设置GDAL缓存大小为500MB
Sys.setenv(GDAL_CACHEMAX = 500)
# 设置超时时间为300秒
Sys.setenv(VSICURL_TIMEOUT = 300)

设置后重新尝试读取大范围数据。

2. 转换stars对象为terra栅格
   利用read_stars()可正常读取的特性，先读取为stars对象，再转换为terra的SpatRaster：

yy <- read_stars(bigger_dim)
xx <- rast(yy)

3. 分块读取后合并
   将大范围请求拆分为多个小范围子请求，读取完成后合并为完整栅格：

# 示例：拆分X轴范围为两个子区域
x_part1 <- rast('https://www.mrlc.gov/geoserver/rcmap_anhb/wcs?request=GetCoverage&service=WCS&version=2.0.1&coverageid=rcmap_anhb__rcmap_annual_herbaceous_2020&subset=X(-1054529.13659319,-1038564.59503519)&subset=Y(2020307.11783295,2056926.21900819)')
x_part2 <- rast('https://www.mrlc.gov/geoserver/rcmap_anhb/wcs?request=GetCoverage&service=WCS&version=2.0.1&coverageid=rcmap_anhb__rcmap_annual_herbaceous_2020&subset=X(-1038564.59503519,-1022600.05347719)&subset=Y(2020307.11783295,2056926.21900819)')
xx <- merge(x_part1, x_part2)

4. 更新依赖包与GDAL
   旧版本的terra或GDAL可能存在WCS请求的大小限制问题，尝试更新到最新版本：

# 更新terra包
install.packages("terra")

若使用系统级GDAL，需自行更新对应系统的GDAL库。

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