读取多NetCDF文件绘制海冰浓度时间序列报错求解决方案

错误原因 #

你的代码核心问题是错误地将单个变量的DataArray逐一收集后尝试合并，而非保留每个文件的Dataset结构进行合并。每个NetCDF文件中的time本身是坐标变量，当你把time、xgrid、ygrid、cdr_seaice_conc_monthly作为独立DataArray加入列表后，xr.concat会尝试将这些不同变量的数组沿time维度合并，导致time变量/坐标重复冲突，触发ValueError: time already exists as coordinate or variable name。

另外需要注意：xgrid和ygrid是极射投影坐标（单位为米），不是地理经纬度，直接用longitude=-90、latitude=-70调用sel会定位错误，需要先转换经纬度到投影坐标。

解决方案一：使用xr.open_mfdataset（推荐，高效简洁） #

xarray提供了open_mfdataset专门用于批量读取并合并多NetCDF文件，无需手动循环，自动处理时间维度合并：

import xarray as xr
import matplotlib.pyplot as plt
from pyproj import Transformer

# 批量读取文件，只加载需要的变量
directory = r'E:\Southern Hemisphere'
ds = xr.open_mfdataset(
    f'{directory}/*.nc',
    combine='by_coords',  # 按坐标自动合并
    parallel=True,  # 并行读取（可选，加快大文件读取）
    variables=['cdr_seaice_conc_monthly', 'time', 'xgrid', 'ygrid']
)

# 转换目标经纬度到投影坐标（xgrid/ygrid的投影是NSIDC Sea Ice Polar Stereographic South）
proj_transformer = Transformer.from_crs(
    "EPSG:4326",  # WGS84经纬度
    "EPSG:3031",  # 南极极射投影
    always_xy=True
)
# 目标经纬度：lon=-90, lat=-70
target_x, target_y = proj_transformer.transform(-90, -70)

# 提取指定投影坐标的海冰浓度时间序列
seaice_conc = ds['cdr_seaice_conc_monthly'].sel(
    xgrid=target_x,
    ygrid=target_y,
    method='nearest'
)

# 过滤无效值（陆地区域、缺失数据等）
seaice_conc = seaice_conc.where(seaice_conc <= 100)

# 绘制时间序列
plt.figure(figsize=(12,6))
seaice_conc.plot(x='time', label='Sea Ice Concentration')
plt.xlabel('时间')
plt.ylabel('海冰浓度')
plt.title('南半球指定位置海冰浓度月度时间序列（1978-2022）')
plt.legend()
plt.show()

解决方案二：手动循环合并Dataset（适合自定义处理场景） #

如果需要手动控制每个文件的读取逻辑，可以收集每个文件的Dataset（仅保留所需变量）后合并：

import os
import xarray as xr
import matplotlib.pyplot as plt
from pyproj import Transformer

directory = r'E:\Southern Hemisphere'
file_names = [os.path.join(directory, f) for f in os.listdir(directory) if f.endswith('.nc')]
# 按文件名排序（确保时间顺序正确，避免文件乱序）
file_names.sort()

datasets = []
for fname in file_names:
    # 打开文件并仅加载需要的变量
    with xr.open_dataset(fname) as ds_single:
        ds_subset = ds_single[['cdr_seaice_conc_monthly', 'time', 'xgrid', 'ygrid']]
        datasets.append(ds_subset)

# 沿时间维度合并所有Dataset
ds_combined = xr.concat(datasets, dim='tdim').rename({'tdim': 'time'})

# 转换经纬度到投影坐标
proj_transformer = Transformer.from_crs("EPSG:4326", "EPSG:3031", always_xy=True)
target_x, target_y = proj_transformer.transform(-90, -70)

# 提取时间序列并过滤无效值
seaice_conc = ds_combined['cdr_seaice_conc_monthly'].sel(
    xgrid=target_x,
    ygrid=target_y,
    method='nearest'
).where(ds_combined['cdr_seaice_conc_monthly'] <= 100)

# 绘图
plt.figure(figsize=(12,6))
plt.plot(ds_combined['time'], seaice_conc)
plt.xlabel('时间')
plt.ylabel('海冰浓度')
plt.title('南半球指定位置海冰浓度月度时间序列（1978-2022）')
plt.show()

关键说明 #

1. 文件排序：务必对文件列表排序，确保时间序列的顺序正确，避免因文件名乱序导致时间轴错位。
2. 投影坐标转换：NSIDC的南半球海冰数据使用EPSG:3031南极极射投影，必须将地理经纬度转换为该投影下的坐标才能正确定位。
3. 数据质量控制：数据中的flag_values（251-255）代表无效值（陆地区域、缺失数据等），通过where过滤后可获得有效浓度值（0-100）。

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