如何校正Raster Stack方向并按行政区提取年降雨数据至Excel？

一、NetCDF转Raster Stack的方向校正 #

问题原因 #

手动构建栅格时，若NetCDF中的纬度是从北到南降序排列，直接按原矩阵创建栅格会导致图像上下颠倒；若经度顺序异常则会左右颠倒。手动循环构建栅格容易忽略维度顺序的匹配问题。

解决方案1：直接用stack()读取NetCDF（推荐） #

raster包的stack()函数可自动识别NetCDF的维度和坐标信息，避免手动构建的错误：

library(raster)

# 直接读取NetCDF为RasterStack，指定目标变量名
raster_stack <- stack("RF25_ind1955_rfp25.nc", varname = "RAINFALL")

# 查看栅格的坐标与维度信息，确认方向问题
print(raster_stack)

# 根据图像方向问题校正：
# 若上下颠倒（纬度方向错误），翻转y轴
raster_stack_corrected <- flip(raster_stack, direction = "y")
# 若左右颠倒（经度方向错误），翻转x轴
# raster_stack_corrected <- flip(raster_stack, direction = "x")

# 验证校正结果
plot(raster_stack_corrected)

解决方案2：手动构建时修正维度顺序 #

如果必须手动提取NetCDF数据构建栅格，先检查纬度顺序再调整矩阵：

library(ncdf4)
library(raster)

ncfile <- nc_open("RF25_ind1955_rfp25.nc")
variable_name <- "RAINFALL"
rainfall_data <- ncvar_get(ncfile, variable_name)
lon <- ncvar_get(ncfile, "LONGITUDE")
lat <- ncvar_get(ncfile, "LATITUDE")
nc_close(ncfile) # 操作完成后关闭NetCDF文件

# 检查纬度顺序：若为从北到南（降序），反转矩阵行
if (lat[1] > lat[length(lat)]) {
  rainfall_data <- rainfall_data[nrow(rainfall_data):1, , ]
}

# 构建RasterStack
raster_stack <- stack()
for (i in 1:dim(rainfall_data)[3]) {
  raster_layer <- raster(rainfall_data[,,i],
                         xmn = min(lon), xmx = max(lon),
                         ymn = min(lat), ymx = max(lat),
                         crs = "+proj=longlat +datum=WGS84")
  raster_stack <- addLayer(raster_stack, raster_layer)
}

plot(raster_stack)

二、提取年降雨量并导出为Excel #

操作步骤及代码 #

假设你有1955-2023年的日降雨量tif文件，按以下流程处理：

library(raster)
library(sf)
library(writexl)
library(lubridate)

# 1. 读取所有日降雨量文件（假设文件命名为rain_YYYYMMDD.tif格式）
rain_files <- list.files(path = "你的日降雨文件文件夹路径", 
                         pattern = "\\.tif$", full.names = TRUE)
rain_stack <- stack(rain_files)

# 2. 匹配日期并按年份分组
# 从文件名提取日期（需根据你的实际命名规则调整正则表达式）
dates <- ymd(sub("rain_(\\d{8})\\.tif", "\\1", basename(rain_files)))
years <- year(dates) # 提取年份用于分组计算

# 3. 计算每年总降雨量
annual_rain <- stackApply(rain_stack, indices = years, fun = sum, na.rm = TRUE)
names(annual_rain) <- unique(years) # 给年图层命名为对应年份

# 4. 读取行政区Shapefile并统一坐标系统
admin_shp <- st_read("你的行政区Shapefile路径.shp")
admin_shp <- st_transform(admin_shp, crs = crs(annual_rain)) # 确保与栅格CRS一致

# 5. 提取每个行政区的年降雨量（zonal统计，高效计算区域总和/均值）
# 若需要区域总降雨量用sum，平均降雨量用mean，按需调整fun参数
annual_rain_stats <- zonal(annual_rain, admin_shp, fun = sum, na.rm = TRUE)

# 6. 合并结果并导出为Excel
# 移除几何列，转为纯数据框后导出
rain_result <- cbind(st_drop_geometry(admin_shp), annual_rain_stats[, -1])
write_xlsx(rain_result, "行政区年降雨量统计.xlsx")

注意事项 #

• 文件名的日期格式需统一，若命名规则不同，需修改sub()中的正则表达式以正确提取日期。
• 若数据量极大，推荐使用terra包替代raster，处理速度更快。
• zonal()函数的fun参数根据需求选择：区域总降雨量用sum，区域平均降雨量用mean。

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