R中NetCDF提取数据无法转为Raster对象的技术问询

我想用Copernicus数据库的NetCDF文件绘制北纬10°、东经40°点的Hovmöller图，数据覆盖1993年至2021年6月的月度表层数据（仅1个深度层）。通过RNetCDF包提取目标变量SelData后，执行DataRast <- raster(SelData)时触发以下错误：

> DataRast <- raster(SelData) Error in (function (classes, fdef, mtable)  :  unable to find an inherited method for function ‘raster’ for signature ‘"numeric"’

尝试将SelData转为矩阵后创建Raster对象（DataRast<- raster(as.matrix(SelData))），但生成的矩阵全为NA值，问题仍未解决。

1. 类型不匹配：raster()函数无法直接处理一维数值向量，而SelData提取的是单个经纬度点的全时间序列数据，属于一维结构，直接传入会触发类型错误。
2. 维度顺序错误：NetCDF变量的维度顺序为time(0), depth(1), lat(2), lon(3)，但原代码提取时用了lon, lat, depth, time的顺序，导致数据结构混乱，转矩阵后出现全NA。
3. 冗余操作：绘制Hovmöller图无需先转成栅格对象，直接整理成数据框更高效，栅格转换属于多余步骤。

步骤1：修正数据提取维度顺序 #

按照NetCDF的维度定义调整var.get.nc的start和count参数，提取[时间×深度]的二维数组：

• start设为c(1, 1, ReqLatI, ReqLonI)（对应time、depth、lat、lon的起始索引）
• count设为c(NA, NA, 1, 1)（提取全时间、全深度、单个经纬度点）

步骤2：直接构建数据框 #

用expand.grid生成时间-深度的网格组合，再绑定数据值，跳过冗余的栅格转换步骤。

步骤3：适配单深度层绘图逻辑 #

由于只有1个深度层，可选择绘制更直观的时间序列折线图，或保留原轮廓图格式。

library(RNetCDF)
library(dplyr)
library(ggplot2)

FileNC <- "cmems_mod_glo_phy_my_0.083deg_P1M-m_1744109649217.nc"

# 目标深度层（1=表层）
ReqLev <- 1
# 目标经纬度
ReqLat <- 10  # 北纬
ReqLon <- 40  # 东经

# 读取NetCDF文件
DataFile <- open.nc(FileNC)
print.nc(DataFile)

# 目标变量（海表温度）
VarName <- 'thetao'

# 提取基础维度数据
Lat <- var.get.nc(DataFile, 'latitude')
Lon <- var.get.nc(DataFile, 'longitude')
Depth <- var.get.nc(DataFile, 'depth')
Time <- var.get.nc(DataFile, 'time')

# 转换时间格式
TimeUnit <- att.get.nc(DataFile, "time", "units")
Date <- utcal.nc(TimeUnit, Time, "c")
DateCh <- format(Date, format="%Y-%m")

# 获取变量属性
VarLongName <- att.get.nc(DataFile, VarName, "long_name")
VarUnits    <- att.get.nc(DataFile, VarName, "units")
Credit <- att.get.nc(DataFile, 'NC_GLOBAL', 'credit')

# 确认维度顺序
var.inq.nc(DataFile, 'thetao')
dim.inq.nc(DataFile, 0)  # 维度0：time
dim.inq.nc(DataFile, 1)  # 维度1：depth
dim.inq.nc(DataFile, 2)  # 维度2：lat
dim.inq.nc(DataFile, 3)  # 维度3：lon

# 找到最接近目标经纬度的索引
ReqLatI <- which.min(abs(Lat-ReqLat))
ReqLonI <- which.min(abs(Lon-ReqLon))

# 按正确维度顺序提取数据：[time, depth]二维数组
SelData <- var.get.nc(DataFile, VarName,
                      start=c(1, 1, ReqLatI, ReqLonI),
                      count=c(NA, NA, 1, 1))

# 关闭NetCDF文件，释放资源
close.nc(DataFile)

# 构建时间-深度网格数据框
GridDF <- expand.grid(Time = Time, Depth = Depth)
GridDF$Data <- as.vector(SelData)

# 转换时间为日期格式
GridDF <- GridDF %>%
  mutate(DateCT = utcal.nc(TimeUnit, Time, "c"),
         Date = as.Date(DateCT))

# 筛选表层数据（可选）
GridDF_surface <- GridDF %>% filter(Depth == Depth[ReqLev])

# 选项1：单深度层时间序列折线图
ggplot(data=GridDF_surface, aes(x=Date, y=Data)) +
  geom_line(color='steelblue') +
  geom_point(shape=3, color='darkred', size=1) +
  coord_cartesian(expand=FALSE) +
  labs(y = paste(VarLongName, " (", VarUnits, ")", sep="")) +
  theme_light() +
  scale_x_date(date_breaks = "3 year", date_labels = "%Y") +
  labs(title = paste(VarLongName, " - Time Series"),
       subtitle = paste('Lat=', Lat[ReqLatI], 'º, Lon=', Lon[ReqLonI],'º', sep=''),
       x = "Date",
       caption = paste('Source:', Credit))

# 选项2：保留Hovmöller轮廓图格式
ggplot(data=GridDF, aes(x=Date, y=Depth, z=Data)) +
  geom_contour_filled() +
  geom_point(shape=3, color='azure3', size=0.5) +
  coord_cartesian(expand=FALSE) +
  scale_y_reverse() +
  labs(fill = VarUnits) +
  theme_light() +
  scale_x_date(date_breaks = "3 year", date_labels = "%Y") +
  labs(title = VarLongName,
       subtitle = paste('Hovemöller Diagram (Lat=', Lat[ReqLatI],
                        'º, Lon=', Lon[ReqLonI],'º)', sep=''),
       x = "Date",
       y = "Depth",
       caption = paste('Source:', Credit))

• 调整了NetCDF数据提取的维度顺序，确保得到二维数组而非一维向量。
• 移除冗余的栅格转换步骤，直接用expand.grid构建数据框，避免类型错误。
• 添加了关闭NetCDF文件的操作，释放系统资源。
• 针对单深度层的情况，提供两种绘图选项适配不同需求。

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