基于hgt.mon.mean.nc数据绘制500hpa高度场

读取文件格式

import xarray as xr
file_path = 'D:/data/hgt.mon.mean.nc'  
  
# 使用xarray打开NetCDF文件  
ds = xr.open_dataset(file_path)  
print(ds)

<xarray.Dataset>
Dimensions:  (level: 17, lat: 73, lon: 144, time: 883)
Coordinates:
  * level    (level) float32 1e+03 925.0 850.0 700.0 ... 50.0 30.0 20.0 10.0
  * lat      (lat) float32 90.0 87.5 85.0 82.5 80.0 ... -82.5 -85.0 -87.5 -90.0
  * lon      (lon) float32 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
  * time     (time) datetime64[ns] 1948-01-01 1948-02-01 ... 2021-07-01
Data variables:
    hgt      (time, level, lat, lon) float32 ...
Attributes:
    description:     Data from NCEP initialized reanalysis (4x/day).  These a...
    platform:       Model
    Conventions:    COARDS
    NCO:            20121012
    history:        Created by NOAA-CIRES Climate Diagnostics Center (SAC) fr...
    title:          monthly mean hgt from the NCEP Reanalysis
    dataset_title:  NCEP-NCAR Reanalysis 1
    References:     http://www.psl.noaa.gov/data/gridded/data.ncep.reanalysis... 

从文件信息可以看出，hgt.mon.mean.nc的文件分别有四个变量

level 高度 总共17层，范围为10-1000hpa

lat 纬度 范围-90-90

lon 经度 范围0-360

time 时间 1948年1月1日-2021年7月1日

import xarray as xr  
import matplotlib.pyplot as plt
import numpy as np  

file_path = 'D:/data/hgt.mon.mean.nc'  
  
# 使用xarray打开NetCDF文件  
ds = xr.open_dataset(file_path)  
  
# 查找包含高度场数据的变量，通常是'hgt'  
hgt_var = ds['hgt']  
  
# 检查是否有压力层（level）维度，并找到500hPa对应的索引  
if 'level' in hgt_var.dims:  
    pressure_levels = ds['level'].values  # 假设压力层存储在'lev'变量中  
    level_index = np.argmin(np.abs(pressure_levels - 500))  
    hgt_500hPa = hgt_var.isel(level=level_index)  
else:  
    raise ValueError("数据集中没有压力层维度，无法直接选择500hPa层")  
  
# 打印时间坐标以确认索引  
print(hgt_500hPa.time)  
  
# 选择一个具体的时间点来绘制，例如2020年1月  
# 注意：这里的时间格式需要与你的数据集中的时间格式相匹配  
time_index = '2020-01'  
try:  
    # 尝试选择时间，并立即使用squeeze去除可能的单一维度  
    hgt_500hPa_jan = hgt_500hPa.sel(time=time_index, method='nearest').squeeze()  
      
    # 检查数据维度  
    if hgt_500hPa_jan.ndim != 2:  
        raise ValueError("提取后的数据维度不是二维的，尽管已经尝试使用squeeze")  
  
    # 绘制等值线图  
    plt.figure(figsize=(10, 5))  
    contour_levels = np.arange(hgt_500hPa_jan.min().values, hgt_500hPa_jan.max().values + 10, 10)  
    cf = plt.contourf(hgt_500hPa_jan.lon, hgt_500hPa_jan.lat, hgt_500hPa_jan, levels=contour_levels, cmap='viridis', extend='both')  
    cs = plt.contour(hgt_500hPa_jan.lon, hgt_500hPa_jan.lat, hgt_500hPa_jan, colors='black', linewidths=0.5)  
    plt.colorbar(cf, label='Height (m)')  
    plt.clabel(cs, inline=True, fontsize=10)  
    plt.title('500hPa Mean Height Field Contours in January 2020')  
    plt.xlabel('Longitude')  
    plt.ylabel('Latitude')  
    plt.show()  
except KeyError:  
    print(f"时间索引 '{time_index}' 不存在于数据集中")  
except ValueError as e:  
    print(e)