Python用Basemap画位势高度距平图和5880线

import xarray as xr
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.basemap import Basemap
plt.rcParams['font.sans-serif'] = ['KaiTi','SimHei']
plt.rcParams['font.size']= 12
plt.rcParams['axes.unicode_minus']=False
hgt=r"D:\论文\论文数据\nc文件1998hgt.mon.mean.nc"  #读取1998年目标文件
def DrawSeaLPJuly():
    Hgt = xr.open_dataset(hgt)["hgt"]
    time_label = (Hgt["time"].dt.month == [6,7,8])    #选定6，7，8三个月
    Hgt_main=Hgt[time_label].loc[dict(level=500)].mean("time").loc[60:0,70:180]  # 选定气压为500hpa，选定经纬度为0~60N，70~180E
    lat = np.array(Hgt_main['lat'])  # 设置好经纬度
    lon=np.array(Hgt_main['lon'])  # 设置好经纬度
   #  算气候状态（1981~2010）30年：
    i=1981
    hgt1=r"D:\论文\论文数据\nc文件1998hgt.mon.mean.nc"
    Hgt1= xr.open_dataset(hgt1)["hgt"]
    time_label = (Hgt1["time"].dt.month == [6, 7, 8])
    Hgt_total = Hgt1[time_label].loc[dict(level=500)].mean("time").loc[60:0,70:180]
    i=i+1
    while i<=2010:
        hgt1=r"D:\论文\论文数据\nc文件"+str(i)+"hgt.mon.mean.nc"
        Hgt1=xr.open_dataset(hgt1)["hgt"]
        time_label = (Hgt1["time"].dt.month == [6,7,8])
        Hgt_total = Hgt1[time_label].loc[dict(level=500)].mean("time").loc[60:0,70:180]+Hgt_total
        i=i+1
    Hgt_total=Hgt_total/30  #求平均
    con=Hgt_main-Hgt_total
    X, Y = np.meshgrid(lon, lat)
    fig = plt.figure(figsize=[12, 18])
    ax = fig.add_subplot(111)
    m = Basemap(lat_0=40, lon_0=130, projection="cyl",
                llcrnrlon=70,  # 地图左边经度
                llcrnrlat=0.0,  # 地图下方纬度
                urcrnrlon=180.0,  # 地图上方经度
                urcrnrlat=60.0,  # 地图上方纬度
                resolution=None,  # 分辨率，这里设置为无
                )
    m.readshapefile(r"D:\论文\论文数据\Shp\世界国家SHP\世界国家", "World Map", color="k", linewidth=1.2)  # 读取shp文件，想放多少就写多少这个代码
    m.readshapefile(r"D:\论文\论文数据\Shp\全国水系矢量数据\hyd1_4l", "river Map", color="lightskyblue", linewidth=1)  #river Map是名字，可任意换，但必须有
    levels = np.arange(-30, 70, 12)
    c = m.contourf(X, Y, con,cmap="RdBu_r",alpha=0.7)  # 绘制等值线
    a1=m.contour(X, Y, con, colors="k", linewidths=0.5)
    # ax.clabel(a2, fmt="%2.0f", fontsize=10, colors="k")  #色标
    a2=m.contour(X,Y, Hgt_main, range(5880, 5890, 10),linewidths=0.5,colors='r')   #画出特定等值线，range（）函数画出其实和重点线，以多少为间隔
    ax.clabel(a2, fmt="%2.0f", fontsize=10, colors="r")
    ax.clabel(a1, fmt="%2.0f", fontsize=10, colors="g")  #标出值
    # fig.colorbar(c, ax=ax, orientation='horizontal')
    plt.yticks([0, 10, 20, 30, 40, 50, 60], ["EQ", "10°N", "20°N", "30°N", "40°N", "50°N", "60°N"], fontsize=16)  #设置x和y坐标轴
    plt.xticks([70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180],
               ["70°E", "80°E", "90°E", "100°E", "110°E", "120°E", "130°E", "140°E", "150°E", "160°E", "170°E", "180°E"],
               fontsize=16)
    plt.ylabel("纬度", fontsize=16)
    plt.xlabel("经度", fontsize=16)
    m.drawparallels(np.arange(0, 60, 10))
    m.drawmeridians(np.arange(70, 181, 10))  # 画出格线
    plt.title("500hpa位势高度距平场", fontsize=18)  # 写上title

    plt.show()


if __name__ == '__main__':
    DrawSeaLPJuly()

本文数据为NCEP/NCAR 再分析数据，感觉这个真是比ERA5好用，就是分辨率和内容类型比ERA5差。链接如下：https://psl.noaa.gov/data/gridded/data.ncep.reanalysis.html。本文比较初级，就是画了等值线还有一个画出目标区域的线，shp文件网上有很多，想要可以发邮箱在评论区，我把我收藏的都给可以奉上。文中我把一个大文件裁剪成30个小文件，毫无必要，这个可以再优化下。希望能给大家提供帮助。