R中如何正确读取.gz压缩的ASCII栅格(asc.gz)文件

问题说明 #

运行环境为R 4.1.1，需处理德国气象局公开的gzip压缩ASCII栅格（.asc.gz）月降水数据，要求全流程自动处理、无需手动解压，批量提取逐点位逐年月气象值。
运行现有公开方案时触发以下报错：

Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer",  : 
  Cannot create a RasterLayer object from this file.

unable to find an inherited method for function ‘raster’ for signature ‘"gzfile"’

已尝试方案的失效原因 #

• 直接传入gzfile()连接读取：raster::raster()没有针对gzfile连接的读取方法，且ASCII栅格读取需要随机访问文件内容，流式gz连接仅支持顺序读取，无法满足寻址要求。
• 调用unzip()解压读取：混淆压缩格式，unzip()仅支持.zip格式压缩包，无法识别gzip（.gz）格式，且代码逻辑中未将解压后的.asc文件作为读取目标。
• 读取原始二进制向量后传入连接解析：读取完二进制内容后连接指针已移至末尾，且raster()无法直接从消耗过的连接中解析栅格结构。

可行实现代码 #

前置依赖 #

用到raster包处理栅格，可选安装R.utils包简化gzip解压操作，也可使用base R原生函数实现无依赖解压。

单文件读取（R.utils简化版） #

library(raster)
library(R.utils)

# 目标压缩文件路径
ras_path <- file.path(myPath, 'rawData/DeutschWetter/01_Jan', "200501asc.gz")

# 生成后缀为.asc的临时文件
temp_asc <- tempfile(fileext = ".asc")

# 自动解压到临时文件，*remove=FALSE保证原压缩包不被修改/删除*
gunzip(ras_path, destname = temp_asc, remove = FALSE)

# 读取栅格
r <- raster(temp_asc)
# 赋值DWD数据默认坐标系（DHDN3 Gauss-Krüger 3度带，EPSG:31467），文件头未自带投影需手动添加
crs(r) <- CRS("+init=epsg:31467")

# 读取完成后删除临时文件
unlink(temp_asc)

批量全量自动读取 #

遍历所有目录下的.asc.gz文件，自动完成解压、读取、堆叠，全程无人工操作：

library(raster)
library(R.utils)

# 递归获取所有.asc.gz格式文件的完整路径
all_gz <- list.files(
  path = file.path(myPath, "rawData/DeutschWetter"),
  pattern = "\\.asc\\.gz$",
  recursive = TRUE,
  full.names = TRUE
)

# 循环读取所有栅格
ras_list <- lapply(all_gz, function(gz_path){
  temp_f <- tempfile(fileext = ".asc")
  gunzip(gz_path, destname = temp_f, remove = FALSE)
  ras <- raster(temp_f)
  crs(ras) <- CRS("+init=epsg:31467")
  unlink(temp_f)
  return(ras)
})

# 可直接堆叠为多层栅格，方便后续点位提取
ras_stack <- stack(ras_list)

无第三方包依赖实现（base R原生） #

不想安装R.utils的话，可通过原生gz连接分块读写完成解压，大文件下内存占用更低：

library(raster)

ras_path <- file.path(myPath, 'rawData/DeutschWetter/01_Jan', "200501asc.gz")
temp_asc <- tempfile(fileext = ".asc")

# 建立读写连接，分块传输
con_in <- gzfile(ras_path, open = "rb")
con_out <- file(temp_asc, open = "wb")
while(length(chunk <- readBin(con_in, what = "raw", n = 10^6)) > 0){
  writeBin(chunk, con_out)
}
close(con_in)
close(con_out)

# 读取栅格
r <- raster(temp_asc)
crs(r) <- CRS("+init=epsg:31467")
unlink(temp_asc)

提示：如果后续需要提取固定点位的时间序列值，可直接用raster::extract()传入点位空间对象和堆叠后的ras_stack，一次性输出所有点位全时段结果，效率远高于逐文件循环提取。

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