从.nc文件绘制浮游生物生物量时间序列图及空图问题排查

我有一个存储浮游生物生物量的.nc文件，包含lon（经度，180个值）、lat（纬度，125个值）、time（时间，480个值）三个维度。时间为数值型，取值1到480，对应1980-2000年。我想绘制生物量随年份变化的图，有两个疑问：

1. 是否需要对同一时间点的所有经纬度数据取平均值？
2. 运行以下R代码后得到空图，问题出在哪里？

用户提供的R代码：

library(ncdf4)
library(dplyr)
library(ggplot2)

our_nc_data <- nc_open("corr_eco3m_bs44_conversion_factor_v7_withDetritus_MEECE.nc")

lat <- ncvar_get(our_nc_data, "latitude")
lon <- ncvar_get(our_nc_data, "longitude")
time <- ncvar_get(our_nc_data, "time")

nt <- dim(time)

start_date <- as.Date("1980-01-01")
time_seq <- seq(start_date, by = "15 days", length.out = nt)

biomass_array <- ncvar_get(our_nc_data, "Zo")
fillvalue <- ncatt_get(our_nc_data, "Zo", "_FillValue")

biomass_array[biomass_array == fillvalue$value] <- NA

lonlattime <- as.matrix(expand.grid(lon, lat, time_seq))

biomass_vec_long <- as.vector(biomass_array)

biomass_obs <- data.frame(cbind(lonlattime, biomass_vec_long))

colnames(biomass_obs) <- c("Long", "Lat", "Date", "Biomass_Zo")

biomass_final <- na.omit(biomass_obs)

biomass_final <- biomass_final[-c(1:2)]

biomass_final$Date <- as.Date(biomass_final$Date)
biomass_final$Biomass_Zo <- as.double(biomass_final$Biomass_Zo)

# Take the mean 
biomass_final <- biomass_final %>%
  group_by(Date) %>%
  summarize(Mean_biomass = mean(Biomass_Zo, na.rm = TRUE))

write.csv(as.data.frame(lswt_final), "Zo_biomass.csv", row.names=T)

# Plot the data using ggplot2
ggplot(biomass_final, aes(x = Date, y = Mean_biomass)) +
  geom_line(size = 0.5) +
  labs(x = "Date", y = "Biomass_Zo") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

示例数据：

our_nc_data <- list(
  filename = "corr_eco3m_bs44_conversion_factor_v7_withDetritus_MEECE.nc",
  writable = FALSE,
  id = 262144L,
  error = FALSE,
  safemode = FALSE,
  format = "NC_FORMAT_CLASSIC"
  )

biomass_final <- structure(list(Date = structure(c(3652, 3667, 3682, 3697, 3712, 
3727), class = "Date"), Mean_biomass = c(7.40254209928812e+35, 
159.139770071623, 262.03385048119, 275.716041716774, 48.992811774022, 
7.42591774539943)), row.names = c(NA, -6L), class = c("tbl_df", 
"tbl", "data.frame"))

1. 是否需要对同一时间点的经纬度取平均？ #

是的。如果你的目标是展示整个研究区域内浮游生物生物量随时间的整体变化趋势，必须对同一时间点的所有经纬度网格单元的生物量取平均值。原始数据是每个经纬度网格在对应时间点的生物量，直接绘图会因数据量过大（总计900万条记录）无法呈现清晰趋势，取平均后每个时间点对应一个区域平均生物量，才能得到直观的时间变化曲线。

2. 空图问题的原因及修正 #

你的代码存在三个关键问题：

• 变量名错误：write.csv中使用了未定义的lswt_final，正确变量应为biomass_final；
• 未彻底处理填充值：示例数据中出现的7.40254209928812e+35是未被正确替换的填充值，极大值会拉宽y轴范围，导致正常数值的线条被压缩至不可见；
• 维度顺序错位：expand.grid(lon, lat, time_seq)的顺序可能与biomass_array的维度顺序不匹配，导致数据错位，产生大量无效值。

修正后的代码 #

library(ncdf4)
library(dplyr)
library(ggplot2)

# 读取nc文件
our_nc_data <- nc_open("corr_eco3m_bs44_conversion_factor_v7_withDetritus_MEECE.nc")

# 获取维度数据
lat <- ncvar_get(our_nc_data, "latitude")
lon <- ncvar_get(our_nc_data, "longitude")
time <- ncvar_get(our_nc_data, "time")

# 生成时间序列：1980-2000共21年，480个时间点（间隔约15天）
start_date <- as.Date("1980-01-01")
time_seq <- seq(start_date, by = "15 days", length.out = length(time))

# 读取生物量变量，先通过print(our_nc_data)确认维度顺序（通常为lon, lat, time或lat, lon, time）
biomass_array <- ncvar_get(our_nc_data, "Zo")
fillvalue <- ncatt_get(our_nc_data, "Zo", "_FillValue")$value

# 处理填充值：用浮点误差范围匹配，避免精准匹配的遗漏
biomass_array[abs(biomass_array - fillvalue) < 1e-6] <- NA

# 直接对三维数组按时间维度（第3维）计算区域平均，高效且避免维度错位
# 若维度顺序为(lat, lon, time)，则将MARGIN改为c(1,2)
mean_biomass <- apply(biomass_array, MARGIN = 3, FUN = mean, na.rm = TRUE)

# 构建最终数据框
biomass_final <- data.frame(
  Date = time_seq,
  Mean_biomass = mean_biomass
)

# 过滤极端异常值：根据示例数据，合理值在几十到几百，过滤大于1000的值
biomass_final <- biomass_final %>% filter(Mean_biomass < 1000)

# 保存结果
write.csv(biomass_final, "Zo_biomass.csv", row.names = FALSE)

# 绘制时间变化图
ggplot(biomass_final, aes(x = Date, y = Mean_biomass)) +
  geom_line(size = 0.5) +
  labs(x = "日期", y = "浮游动物生物量") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

# 关闭nc文件，释放资源
nc_close(our_nc_data)

关键修正说明 #

• 高效计算均值：用apply直接操作三维数组，避免转成数据框的冗余步骤，同时规避维度顺序错位问题；
• 精准处理填充值：通过浮点误差范围匹配，确保所有填充值被替换为NA；
• 过滤异常值：移除明显不合理的极大值，保证y轴范围正常，线条可见；
• 规范资源管理：添加nc_close关闭nc文件，避免资源泄漏。

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