如何用R将ROMS曲线网格变量保存为可识别的NetCDF文件？

处理ROMS海洋模型的曲线网格（curvilinear grid）数据时，要保存变量子集为NetCDF并让stars、gdalwarp等工具识别网格，核心是遵循CF规范：保留二维坐标变量，并给数据变量添加coordinates属性关联坐标。以下是完整实现步骤：

1. 加载依赖包 #

library(ncdf4)
library(stars)

2. 读取原始ROMS数据的关键内容 #

从原始NetCDF中提取需要的温度数据、二维坐标矩阵和网格维度：

filename <- "http://dap.saeon.ac.za/thredds/dodsC/SAEON.EGAGASINI/2019.Penven/DAILY_MEANS/1_12_degree/roms_avg_Y2014M11.nc.1"
outfile <- "woesII_singleDay_temp_vertlevl60.nc"

# 打开原始文件并读取核心内容
dta <- nc_open(filename)
lons <- ncvar_get(dta, 'lon_rho')          # 二维经度矩阵（RHO点）
lats <- ncvar_get(dta, 'lat_rho')          # 二维纬度矩阵（RHO点）
xi_rho <- ncvar_get(dta, "xi_rho")         # 曲线网格x维度索引
eta_rho <- ncvar_get(dta, "eta_rho")       # 曲线网格y维度索引
# 提取第60层、第1个时间步的温度数据
temperature <- ncvar_get(dta, "temp", start = c(1,1,60,1), count = c(-1,-1,1,1))
# 读取原始温度变量的元数据（复用保证一致性）
temp_attrs <- ncatt_get(dta, "temp")
nc_close(dta)

3. 定义NetCDF的维度与变量 #

严格按照ROMS的结构定义维度和变量，确保二维坐标的完整性：

# 定义曲线网格的索引维度
xi_dim <- ncdim_def(name = "xi_rho", units = "", vals = xi_rho, unlim = FALSE)
eta_dim <- ncdim_def(name = "eta_rho", units = "", vals = eta_rho, unlim = FALSE)

# 定义二维坐标变量（必须与原始文件结构一致）
lon_var <- ncvar_def(
  name = "lon_rho", 
  units = "degrees_east", 
  dim = list(xi_dim, eta_dim), 
  missval = NA,
  longname = "Longitude of RHO-points"
)
lat_var <- ncvar_def(
  name = "lat_rho", 
  units = "degrees_north", 
  dim = list(xi_dim, eta_dim), 
  missval = NA,
  longname = "Latitude of RHO-points"
)

# 定义温度数据变量
data_var <- ncvar_def(
  name = "temp", 
  units = temp_attrs$units,  # 复用原始单位
  dim = list(xi_dim, eta_dim), 
  missval = NA,
  longname = temp_attrs$longname  # 复用原始长名称
)

4. 创建NetCDF文件并写入数据（核心步骤） #

重点是给温度变量添加coordinates属性，这是工具识别曲线网格的关键：

# 创建NetCDF文件
ncfile <- nc_create(outfile, vars = list(data_var, lon_var, lat_var))

# 写入数据到文件
ncvar_put(ncfile, data_var, temperature)
ncvar_put(ncfile, lon_var, lons)
ncvar_put(ncfile, lat_var, lats)

# 核心：给温度变量添加coordinates属性，关联二维坐标变量
ncatt_put(ncfile, data_var, "coordinates", "lon_rho lat_rho")

# 可选：复用原始变量的其他属性（如填充值）
if (!is.null(temp_attrs$`_FillValue`)) {
  ncatt_put(ncfile, data_var, "_FillValue", temp_attrs$`_FillValue`)
}
# 添加网格类型属性，增强兼容性
ncatt_put(ncfile, data_var, "grid_type", "curvilinear")

# 添加全局元数据
ncatt_put(ncfile, 0, "title", "ROMS Surface Temperature (Level 60) on Curvilinear Grid")
ncatt_put(ncfile, 0, "source", "Derived from SAEON EGAGASINI ROMS daily mean output")
ncatt_put(ncfile, 0, "CRS", "EPSG:4326")

# 关闭文件
nc_close(ncfile)

5. 验证文件正确性 #

用stars读取验证 #

ncvar2 <- read_ncdf(
  outfile, 
  var = "temp", 
  curvilinear = c("lon_rho", "lat_rho"), 
  make_time = FALSE
) %>% st_set_crs(4326)

# 查看结果
print(ncvar2)
# 绘图验证
plot(ncvar2)

用GDAL验证（需系统安装GDAL） #

system(paste0("gdalinfo ", outfile))

关键注意事项 #

• 必须保留二维的lon_rho和lat_rho变量，不能简化为一维坐标，这是曲线网格的核心特征
• coordinates属性是CF规范要求，是stars、gdal等工具识别曲线网格的关键标识
• 复用原始ROMS文件的维度名称（xi_rho/eta_rho）和元数据，避免兼容性问题

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