如何在R中利用shapefile为netcdf文件添加区域标识变量

1. 依赖包安装与加载 #

首先安装并导入操作所需的核心包，覆盖NetCDF读写、矢量处理、栅格计算三类功能：

# 首次使用安装依赖包
install.packages(c("ncdf4", "sf", "terra"))
# 加载包
library(ncdf4)
library(sf)
library(terra)

2. 生成NetCDF匹配的栅格模板 #

读取叶绿素NetCDF文件，提取经纬度范围、分辨率等信息，生成空栅格模板，保证后续生成的区域变量和原文件坐标完全对齐：

# 读取NetCDF文件
nc <- nc_open("你的叶绿素NetCDF文件路径.nc")
# 提取经纬度坐标
lon <- nc$dim$lon$val
lat <- nc$dim$lat$val
# 生成对齐模板
r_template <- rast(
  nrows = length(lat),
  ncols = length(lon),
  xmin = min(lon), xmax = max(lon),
  ymin = min(lat), ymax = max(lat),
  crs = "EPSG:4326" # 若你的文件使用其他坐标系，修改为对应EPSG编码即可
)
# 关闭临时读取连接
nc_close(nc)

3. 读取并处理Longhurst海洋区域shapefile #

读取shapefile后统一坐标系，给每个区域分配你需要的数值编码：

# 读取shapefile
longhurst <- st_read("你的Longhurst shapefile文件路径.shp")
# 统一坐标系，和栅格模板保持一致
longhurst <- st_transform(longhurst, crs(r_template))

# 分配区域编码
# 方案1：自动按属性表顺序从0开始编码
longhurst$region_code <- 0:(nrow(longhurst)-1)

# 方案2：自定义特定区域编码，示例如下
# longhurst$region_code <- dplyr::case_when(
#   longhurst$ProvName == "Mediterranean Sea" ~ 0,
#   longhurst$ProvName == "North Sea" ~ 1,
#   TRUE ~ 2:(nrow(longhurst)-1) # 其余区域按顺序编码
# )

4. 矢量区域栅格化 #

将shapefile的矢量区域按照分配的编码，转成和NetCDF对齐的栅格数据：

region_rast <- rasterize(
  x = vect(longhurst),
  y = r_template,
  field = "region_code",
  background = NA # 未落在任何区域的点位赋值为NA，可自行修改为其他数值
)

5. 新增region变量写入NetCDF #

将生成的区域编码矩阵写入原NetCDF文件，完成变量新增：

# 以可写模式打开NetCDF文件
nc <- nc_open("你的叶绿素NetCDF文件路径.nc", write = TRUE)
# 提取区域编码矩阵，调整维度匹配NetCDF的经纬度顺序
region_vals <- as.matrix(region_rast, wide = TRUE)
# 定义region变量
region_var <- ncvar_def(
  name = "region",
  units = "无单位",
  dim = list(nc$dim$lon, nc$dim$lat),
  missval = NA,
  longname = "所属Longhurst海洋区域编码"
)
# 新增变量到文件
ncvar_add(nc, region_var)
# 写入数值
ncvar_put(nc, "region", region_vals)
# 关闭文件保存修改
nc_close(nc)

注意事项 #

• 若你的NetCDF包含时间维度，可将region_vals矩阵按时间维度复制拓展后再写入，即可实现每个时间层都对应相同的区域编码
• 可自行导出longhurst[,c("ProvName","region_code")]的对应关系表，方便后续数据分析时对照区域名称和编码
• 若出现区域错位问题，优先检查两个文件的坐标系是否完全一致，以及经纬度的排列顺序是否匹配

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