R raster包crop裁剪GeoTIFF栅格时值被四舍五入与原值不符

问题根因 #

该问题和crop()函数本身无关，核心是栅格读写时的数据类型默认匹配错误：
raster包对GeoTIFF的输出数据类型默认采用自动识别逻辑，经常会将存储浮点数值的栅格误判定为整型（INT1U/INT2S/INT4S等），写入文件时会自动对所有小数位做四舍五入取整；QGIS裁剪时默认完整保留原始栅格的元数据与存储位深，因此不会出现该问题。
另外代码中手动强制覆写栅格CRS的操作，也可能导致raster包读取原始栅格的缩放系数、偏移量等元数据时出现异常，进一步放大数值错误。
同位置栅格值对比如下：

排查与修复步骤 #

• 第一步：验证异常出现的环节
  加载栅格栈后先执行head(getValues(stack_chelsa[[1]]))抽取若干像元值，若此时数值已经是整数，说明是读入环节元数据丢失；若crop()之后、writeRaster()之前抽取的像元值仍为正常小数，可100%确定是写入环节数据类型配置错误。
• 第二步：删除手动覆写CRS的代码
  CHELSA 2.1版本的原始GeoTIFF已经内置正确的WGS84坐标参考信息，不需要手动执行crs(stack_chelsa) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"强制覆写，该操作会覆盖栅格自带的元数据，容易触发参数识别异常。
• 第三步：写入文件时显式指定输出数据类型
  调用writeRaster()时增加datatype = "FLT4S"参数，指定输出为32位浮点型，可完整保留所有小数位；也可以直接在rasterOptions()中全局设置默认输出数据类型，避免后续操作重复配置。

修正后的可运行代码 #

require(raster)
require(rgdal)
# 全局配置临时目录，同时指定默认输出数据类型为32位浮点
rasterOptions(tmpdir="C:/GIS/Temp", progress="text", timer=TRUE, datatype="FLT4S")
setwd("C:/GIS/Saved")
maska <- readOGR(dsn="C:/GIS/Mask", layer="Europe")
pliki_chelsa <- list.files("C:/GIS/Chelsa 2.1/1981-2010/bio/1",pattern='tif',full.names=TRUE)
stack_chelsa <- raster::stack(pliki_chelsa)
nazwy_chelsa <- names(stack_chelsa)
maska_chelsa <- crop(stack_chelsa, maska)
# 写入时显式指定数据类型，可增加LZW压缩减小文件体积
writeRaster(maska_chelsa, 
            filename=nazwy_chelsa, 
            format="GTiff", 
            prj=TRUE, 
            bylayer=TRUE, 
            suffix=nazwy_chelsa, 
            overwrite=TRUE,
            datatype = "FLT4S",
            options = c("COMPRESS=LZW"))

替代方案（兼容性更好） #

如果上述调整后仍存在数值异常，可换用raster包的继任者terra包完成操作，该包对GeoTIFF元数据的识别兼容性更强，默认会完整保留原始栅格的数据类型、缩放系数、坐标参考等信息，不需要手动配置参数即可避免数值取整问题，示例代码如下：

library(terra)
# 读取矢量与栅格
maska <- vect("C:/GIS/Mask/Europe.shp")
pliki_chelsa <- list.files("C:/GIS/Chelsa 2.1/1981-2010/bio/1",pattern='tif',full.names=TRUE)
stack_chelsa <- rast(pliki_chelsa)
# 裁剪与写出
maska_chelsa <- crop(stack_chelsa, maska)
writeRaster(maska_chelsa, 
            filename = paste0(names(stack_chelsa), ".tif"), 
            overwrite=TRUE, 
            gdal=c("COMPRESS=LZW"))

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