如何在R中读取无后缀时序气候数据生成日最高温栅格时间序列

依赖包说明 #

用到两个核心R包，提前安装加载即可：

• data.table：高速读取大体积文本格点数据，性能远优于基础R的读表函数，适合年长时序格点场景
• terra：新一代栅格数据处理包，原生支持带时间维度的栅格时间序列对象，读写、运算效率高

提前运行以下命令安装包（已安装可跳过）：

install.packages(c("data.table", "terra"))

完整实现步骤 #

1. 读取格点经纬度元信息 #

数据前2行从第4列开始依次对应每个格点的经度、纬度，先单独读取这部分元数据：

library(data.table)
library(terra)

# 把此处文件名替换为你本地的无扩展名数据文件路径
raw_file_path <- "your_tmax_raw_file"

# 读取前2行元数据，自动识别分隔符
meta_data <- fread(
  file = raw_file_path,
  nrows = 2,
  header = FALSE,
  sep = "auto"
)

# 提取第4列到末尾的经度、纬度向量
lon <- as.numeric(meta_data[1, 4:ncol(meta_data)])
lat <- as.numeric(meta_data[2, 4:ncol(meta_data)])
grid_point <- data.frame(lon = lon, lat = lat)

# 校验经纬度长度匹配
if (length(lon) != length(lat)) stop("经纬度格点数量不匹配，请检查源数据前2行格式")

2. 读取逐日气温观测值 #

从第3行开始读取数据，前3列为年、月、日，后续列对应每个格点的日最高温值：

# 跳过前2行元数据，读取所有逐日数据
tmax_raw <- fread(
  file = raw_file_path,
  skip = 2,
  header = FALSE,
  sep = "auto"
)

# 拆分日期列与气温值矩阵
date_col <- tmax_raw[, 1:3]
colnames(date_col) <- c("year", "month", "day")
tmax_mat <- as.matrix(tmax_raw[, 4:ncol(tmax_raw)])

# 生成标准日期向量
date_seq <- as.Date(paste(date_col$year, date_col$month, date_col$day, sep = "-"))

# 校验：气温值列数与格点数量一致
if (ncol(tmax_mat) != nrow(grid_point)) stop("气温值列数与格点数量不匹配，请检查源数据")

3. 构建栅格时间序列对象 #

将逐日的格点值转为空间栅格，拼接为带时间维度的栅格栈：

# 初始化栅格存储列表
rast_list <- vector("list", length = nrow(tmax_mat))

# 循环处理每日数据
for (i in seq_along(rast_list)) {
  # 拼合当日坐标与气温值
  day_data <- grid_point
  day_data$tmax <- tmax_mat[i, ]
  # 转为WGS84坐标系的空间点对象
  day_vect <- vect(day_data, geom = c("lon", "lat"), crs = "EPSG:4326")
  # 基于格点自动推断分辨率，栅格化当日气温
  day_rast <- rasterize(day_vect, rast(day_vect), field = "tmax")
  rast_list[[i]] <- day_rast
  
  # 每处理100天打印一次进度
  if (i %% 100 == 0) message(sprintf("已完成%d天数据处理，进度%.1f%%", i, i/length(rast_list)*100))
}

# 拼接为多图层栅格，绑定时间维度
tmax_rast_ts <- rast(rast_list)
time(tmax_rast_ts) <- date_seq
names(tmax_rast_ts) <- paste0("tmax_", date_seq)

结果校验与导出 #

• 校验：随机抽取某一天的图层绘图，检查空间分布是否符合预期，例如抽取第100天的数据：

  plot(tmax_rast_ts[[100]], main = paste0(time(tmax_rast_ts)[100], " 日最高温(℃)"))
  

• 导出：支持两种常用导出方式
  • 导出为带时间维度的NetCDF文件，方便后续气候数据处理：

    writeCDF(tmax_rast_ts, filename = "tmax_daily_timeseries.nc", varname = "tmax", unit = "deg_c", overwrite = TRUE)
    

  • 按日导出为GeoTIFF文件，适配GIS软件打开：

    writeRaster(tmax_rast_ts, filename = paste0("daily_tmax/tmax_", date_seq, ".tif"), overwrite = TRUE)
    

注意事项：如果源数据用特殊值标记缺测（比如-9999、-99999），可以在栅格化前将值矩阵中对应的缺测值替换为NA，避免后续计算出错。

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