使用terra包读取EPSG3035栅格时出现投影数据库错误排查

问题重现 #

使用terra批量读取EPSG:3035投影的栅格文件时，出现以下PROJ相关警告：

proj_create_from_database: datum not found
proj_create_from_database: ellipsoid not found
proj_create_from_database: prime meridian not found

执行代码如下：

Land_variables_files <- list.files("Land_Variables1", pattern = '.tif$', full.names = T)
Land_variables_temp <- rast(Land_variables_files)

同时，show(Land_variables_temp)显示部分图层的最小值/最大值为?，但单独读取或绘图时数值正常；crs(Land_variables_temp)返回正确的EPSG:3035，当前GDAL环境为gdal 3.5.2、proj 8.2.1、geos 3.9.3。

原因分析 #

1. 投影元数据缺失：虽然QGIS中设置了EPSG:3035投影，但导出部分栅格时，可能未写入完整的基准面、椭球、本初子午线等投影元数据。PROJ在批量读取多栅格时会逐个解析文件元数据，缺失信息时触发警告；单独读取时terra会自动补全默认EPSG参数，因此无警告。
2. Min/Max显示问题：这是terra默认行为——批量读取多栅格时，不会自动计算所有图层的统计值，因此显示?；单独读取时会实时计算统计值，所以数值显示正常。

解决方案 #

1. 检查并修复单个栅格的投影元数据 #

先排查具体哪个栅格的元数据存在问题：

# 用terra查看单个文件的详细元数据
describe(Land_variables_files[1])
# 或用GDAL命令行查看更完整的投影信息
system(paste0("gdalinfo ", Land_variables_files[1]))

若发现某文件元数据缺失，重新在QGIS中导出该栅格：

• 确认选择EPSG:3035投影
• 在保存选项中勾选“写入完整的投影信息”（确保导出完整WKT格式投影，而非仅EPSG代码）

2. 强制指定投影读取 #

读取时直接指定CRS，跳过文件元数据的解析步骤：

Land_variables_temp <- rast(Land_variables_files, crs = "EPSG:3035")

3. 批量修复栅格投影元数据 #

用terra重新写入所有栅格，强制写入完整的EPSG:3035投影信息：

# 遍历修复每个文件
for (f in Land_variables_files) {
  r <- rast(f)
  # 生成修复后的文件名
  fixed_f <- gsub(".tif", "_fixed.tif", f)
  writeRaster(r, fixed_f, crs = "EPSG:3035", overwrite = TRUE)
}
# 读取修复后的文件
Land_variables_files_fixed <- list.files(pattern = '_fixed.tif$', full.names = T)
Land_variables_temp <- rast(Land_variables_files_fixed)

4. 解决Min/Max显示问题 #

手动计算所有图层的统计值，或让terra加载所有数据到内存以自动计算：

# 方法1：计算全局统计值
stats <- global(Land_variables_temp, c("min", "max"))
print(stats)

# 方法2：加载所有数据到内存，之后show()会显示统计值
Land_variables_temp <- readAll(Land_variables_temp)
show(Land_variables_temp)

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