在R中转换.GRD二进制文件为NetCDF并进行0.5°分辨率重采样

1. 读取无配套信息的IMD GRD文件并转为NetCDF #

IMD的这类GRD是纯二进制文件，需手动指定数据维度、地理范围等参数，步骤如下：

• 加载依赖包：raster、ncdf4
• 读取二进制数据：IMD温度数据以16位整数存储，需除以10得到实际温度值
• 构建栅格栈：按日期维度排列，指定投影与地理范围
• 导出为NetCDF格式

代码示例：

# 加载包
library(raster)
library(ncdf4)

# 文件路径与基础参数
grd_file <- "Maxtemp_MaxT_1955.GRD"
n_lat <- 31          # 纬度网格数
n_lon <- 31          # 经度网格数
n_days <- 365        # 平年用365，闰年改为366
scale_factor <- 0.1  # 整数转实际温度的系数

# 读取二进制数据并转换为温度值
raw_data <- readBin(grd_file, what = "integer", n = n_lat * n_lon * n_days, endian = "little")
temp_data <- raw_data * scale_factor

# 重塑为[经度, 纬度, 日期]的数组
temp_array <- array(temp_data, dim = c(n_lon, n_lat, n_days))

# 定义地理范围与投影
ext <- extent(67.5, 97.5, 7.5, 37.5) # xmin, xmax, ymin, ymax
crs <- CRS("+proj=longlat +datum=WGS84")

# 构建栅格栈
raster_stack <- stack()
for (d in 1:n_days) {
  # 提取单日数据并转置，适配栅格的纬度顺序
  day_raster <- raster(t(temp_array[,,d]), xmn=ext@xmin, xmx=ext@xmax, ymn=ext@ymin, ymx=ext@ymax, crs=crs)
  # 反转纬度（IMD数据纬度从北到南存储，栅格默认从南到北）
  day_raster <- flip(day_raster, direction='y')
  raster_stack <- addLayer(raster_stack, day_raster)
}

# 为栅格层命名为对应日期
dates <- seq(as.Date("1955-01-01"), as.Date("1955-12-31"), by="day")
names(raster_stack) <- dates

# 导出为NetCDF
writeRaster(raster_stack, filename="Maxtemp_MaxT_1955.nc", format="CDF", overwrite=TRUE)

2. 将栅格重采样至0.5°×0.5°分辨率 #

通过创建0.5°分辨率的模板栅格，使用双线性插值（适合温度这类连续数据）完成重采样：

代码示例：

# 读取NetCDF格式的栅格栈（或直接使用上一步的raster_stack）
raster_stack <- stack("Maxtemp_MaxT_1955.nc")

# 创建0.5°分辨率的模板栅格
template <- raster(ext=extent(67.5, 97.5, 7.5, 37.5), res=0.5, crs=CRS("+proj=longlat +datum=WGS84"))

# 执行重采样
resampled_stack <- resample(raster_stack, template, method="bilinear")

# 验证分辨率（输出应为0.5 0.5）
res(resampled_stack)

后续：年度求和与分区统计 #

完成重采样后，可继续执行年度求和与多边形分区统计：

# 计算年度温度总和
annual_sum <- sum(resampled_stack, na.rm=TRUE)

# 读取要素多边形（示例为shapefile格式）
poly <- shapefile("your_polygons.shp")

# 分区统计取均值
zonal_stats <- zonal(annual_sum, poly, fun="mean", na.rm=TRUE)
print(zonal_stats)

注意事项 #

• 若读取后数据异常，可将readBin的endian参数改为"big"尝试
• 闰年需调整n_days为366，日期序列同步更新
• 离散类数据重采样可改用method="ngb"（最近邻插值）

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