使用R语言Terra包为Sentinel-3 OLCI NetCDF分配CRS失败求助

问题背景 #

这是ACOLITE输出的Sentinel-3 OLCI产品，原生NetCDF文件无CRS信息。在R中使用terra包处理时，即使准备了带WGS84（EPSG:4326）CRS的GeoTIFF模板，调用terra::project()仍报错：

Error: [project] input raster CRS not set.

解决方案 #

1. 为原始NetCDF栅格设置CRS #

Sentinel-3 OLCI数据的原生地理坐标系为WGS84（EPSG:4326），直接为栅格对象赋值CRS：

s3 = rast("..../S3A_OLCI_2022_07_14_15_37_07_L2R.nc")
# 为输入栅格设置CRS
crs(s3) <- "EPSG:4326"
# 验证CRS是否设置成功
crs(s3)

2. 确认并设置栅格的地理范围（可选但关键） #

如果NetCDF文件未自动识别地理范围，需从文件中的经纬度变量提取极值来设置：

# 从NetCDF读取经纬度变量（变量名可能为lat/lon或latitude/longitude，根据ACOLITE输出调整）
lat <- rast("..../S3A_OLCI_2022_07_14_15_37_07_L2R.nc", var="lat")
lon <- rast("..../S3A_OLCI_2022_07_14_15_37_07_L2R.nc", var="lon")

# 设置s3的地理范围
ext(s3) <- c(min(lon[], na.rm=TRUE), max(lon[], na.rm=TRUE), min(lat[], na.rm=TRUE), max(lat[], na.rm=TRUE))

3. 执行投影操作 #

完成CRS和范围设置后，即可使用模板栅格进行投影：

s3_tmp = rast(".../S3A_OLCI_2022_07_14_15_37.tif")
s3_proj <- terra::project(s3, s3_tmp)

关键说明 #

terra::project()要求输入栅格必须具备明确的CRS信息，仅模板栅格有CRS无法完成坐标转换。因此先为无CRS的输入栅格补全坐标系和地理范围是解决报错的核心。

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