天气学原理插图复现（七）——北半球平均经圈环流

使用数据：omega.mon.ltm.1991-2020.nc、vwnd.mon.ltm.1991-2020.nc

使用库：Matplotlib、NumPy、netCDF4、Cartopy

彩色图像 

import matplotlib.pyplot as plt
import netCDF4 as nc
import numpy as np
import cartopy.mpl.ticker as cmt

plt.rcParams['font.sans-serif']= ['Microsoft YaHei'] # 设置“微软雅黑”，图上显示出中文
plt.rcParams['axes.unicode_minus'] = False # 设置中文后，解决坐标轴上负号显示问题


## 等温大气模型
Rd = 287
Tv = 273
g = 9.8
p0 = 1013.2
Hp = Rd*Tv/g
def p2z(p):
    '''p:array'''
    return Hp*np.log(p0/p)

f = nc.Dataset(r"vwnd.mon.ltm.1991-2020.nc",'r')
f_omega = nc.Dataset(r"omega.mon.ltm.1991-2020.nc",'r')
level = f.variables['level'][:12] # omega只有12层
lat = f.variables['lat'][:]
v = f.variables['vwnd'][:]
omega = f_omega.variables['omega'][:]

def omegaa(omega):
    i=len(omega)
    j=len(omega[0])
    k=len(omega[0][0])
    l=len(omega[0][0][0])
    omega02=[]
    for a in range(i):
        omega021=[]
        for b in range(j):
            omega022=[]
            for c in range(k):
                omega023 = omega[a][b][c][30:45]
                omega022.append(omega023)
            omega021.append(omega022)
        omega02.append(omega021)
    omega02=np.array(omega02)
    return omega02
 
omega00=omegaa(omega)


time = nc.num2date(f.variables['time'][:],f.variables['time'].units).data
months = np.array([i.month for i in time])
def pa(para_month):
    v = f.variables['vwnd'][:]
    omega = f_omega.variables['omega'][:]
    omega=omegaa(omega)
    index_month = np.where(months==para_month)[0][0] # 改月份

    #行相加且求平均并切片
    v = v[index_month][:12,:,:].mean(axis=2)[[0,1,2,3,4,5,6,7,9,11],:]
    omega = omega[index_month][:,:,:].mean(axis=2)[[0,1,2,3,4,5,6,7,9,11],:]

    z = p2z(level)
    oz=np.arange(100,1100,100)
    Lat, Z = np.meshgrid(lat,oz)
    return Lat,Z,omega,v

Lat1,Z1,omega1,v1=pa(1)
Lat2,Z2,omega2,v2=pa(7)


fig = plt.figure(figsize=[11,14])
#001
ax = plt.subplot(211)
#绘制阴影区
v1 = v1/np.mean(v1)
omega11=omega1.copy()

omega1 = omega1/np.mean(omega00)

ec = ax.contourf(Lat1[:,::-1],Z1, omega1[::-1,::-1]*1.4,levels=np.arange(omega1.min(),0,6),cmap = 'Oranges_r',extend='both')

cb = fig.colorbar(ec,extend='both',label='hPa/s')
ax2 = cb.ax
ax2.set_yticks(np.arange(omega1.min(),0.001,6),[ -4, -3, -2, -1, 0])
levels = np.arange(-10,20,2)
#二维列表逆序
ax.streamplot(Lat1[:,::-1],Z1,v1[::-1,::-1],omega1[::-1,::-1]*50,color='k',linewidth=2,arrowsize=0.5,density=3)
ax.quiver(Lat1[:,::-1],Z1,v1[::-1,::-1]*25,-omega1[::-1,::-1]*100,scale=150000,width=0.0025, headwidth=8, headlength=3, headaxislength=1)
ax.semilogy([90,90], [1000,100], color='C3')
ax.set_ylim(1000, 100)
ax.set_yticks(np.linspace(100, 1000, 10),[1000,900,800,700,600,500,400,300,200,100][::-1])
ax.set_xlim(-30, 60)
ax.set_xlabel('Latitude[°]')
ax.set_ylabel('level[hPa]')
ax.set_xticks(range(-30,61,10),['30°','20°','10°S', 'EQ','10°','20°', '30°', '40°','50°','60°N'])
ax.set_xlim([-30,60])

plt.title("1月")

#002
ax = plt.subplot(212)
#绘制阴影区
v2 = v2/np.mean(v2)
omega2 = omega2/np.mean(omega00)
ec2 = ax.contourf(Lat2[:,::-1],Z2, omega2[::-1,::-1]*5,levels=np.arange(round(omega2.min(),0),0,19.6),cmap = 'Oranges_r',extend='both')
cb2 = fig.colorbar(ec2,extend='both',label='hPa/s')
ax2 = cb2.ax
ax2.set_yticks(np.arange(round(omega2.min(),0),0,19.6),[-4, -3, -2, -1, 0])
levels = np.arange(-10,20,2)
#二维列表逆序
ax.streamplot(Lat2[:,::-1],Z2,v2[::-1,::-1],omega2[::-1,::-1]*10,color='k',linewidth=2,arrowsize=0.5,density=3)
ax.quiver(Lat2[:,::-1],Z2,v2[::-1,::-1]*150,-omega2[::-1,::-1]*170,scale=180000,width=0.0025, headwidth=8, headlength=3, headaxislength=1)
ax.semilogy([90,90], [1000,100], color='C3')
ax.set_ylim(1000, 100)
ax.set_yticks(np.linspace(100, 1000, 10),[1000,900,800,700,600,500,400,300,200,100][::-1])
ax.set_xlim(-30, 60)
ax.set_xlabel('Latitude[°]')
ax.set_ylabel('level[hPa]')
ax.set_xticks(range(-30,61,10),['30°','20°','10°S', 'EQ','10°','20°', '30°', '40°','50°','60°N'])
ax.set_xlim([-30,60])
plt.title("7月")
plt.suptitle("图4-1-10 75°—110°E的平均经圈环流",fontsize='xx-large')

plt.savefig("图4-1-10.png",dpi=800)

plt.show()

 黑白图像

import matplotlib.pyplot as plt
import netCDF4 as nc
import numpy as np
import cartopy.mpl.ticker as cmt

plt.rcParams['font.sans-serif']= ['Microsoft YaHei'] # 设置“微软雅黑”，图上显示出中文
plt.rcParams['axes.unicode_minus'] = False # 设置中文后，解决坐标轴上负号显示问题


## 等温大气模型
Rd = 287
Tv = 273
g = 9.8
p0 = 1013.2
Hp = Rd*Tv/g
def p2z(p):
    '''p:array'''
    return Hp*np.log(p0/p)

f = nc.Dataset(r"vwnd.mon.ltm.1991-2020.nc",'r')
f_omega = nc.Dataset(r"omega.mon.ltm.1991-2020.nc",'r')
level = f.variables['level'][:12] # omega只有12层
lat = f.variables['lat'][:]
v = f.variables['vwnd'][:]
omega = f_omega.variables['omega'][:]

def omegaa(omega):
    i=len(omega)
    j=len(omega[0])
    k=len(omega[0][0])
    l=len(omega[0][0][0])
    omega02=[]
    for a in range(i):
        omega021=[]
        for b in range(j):
            omega022=[]
            for c in range(k):
                omega023 = omega[a][b][c][30:45]
                omega022.append(omega023)
            omega021.append(omega022)
        omega02.append(omega021)
    omega02=np.array(omega02)
    return omega02
 
omega00=omegaa(omega)


time = nc.num2date(f.variables['time'][:],f.variables['time'].units).data
months = np.array([i.month for i in time])
def pa(para_month):
    v = f.variables['vwnd'][:]
    omega = f_omega.variables['omega'][:]
    omega=omegaa(omega)
    index_month = np.where(months==para_month)[0][0] # 改月份

    #行相加且求平均并切片
    v = v[index_month][:12,:,:].mean(axis=2)[[0,1,2,3,4,5,6,7,9,11],:]
    omega = omega[index_month][:,:,:].mean(axis=2)[[0,1,2,3,4,5,6,7,9,11],:]

    z = p2z(level)
    oz=np.arange(100,1100,100)
    Lat, Z = np.meshgrid(lat,oz)
    return Lat,Z,omega,v

Lat1,Z1,omega1,v1=pa(1)
Lat2,Z2,omega2,v2=pa(7)


fig = plt.figure(figsize=[11,14])
#001
ax = plt.subplot(211)
#绘制阴影区
v1 = v1/np.mean(v1)
omega11=omega1.copy()

omega1 = omega1/np.mean(omega00)

ec = ax.contourf(Lat1[:,::-1],Z1, omega1[::-1,::-1]*1.4,levels=np.arange(omega1.min(),0,6),cmap = 'gray',extend='both')

cb = fig.colorbar(ec,extend='both',label='hPa/s')
ax2 = cb.ax
ax2.set_yticks(np.arange(omega1.min(),0.001,6),[ -4, -3, -2, -1, 0])
levels = np.arange(-10,20,2)
#二维列表逆序
ax.streamplot(Lat1[:,::-1],Z1,v1[::-1,::-1],omega1[::-1,::-1]*50,color='k',linewidth=2,arrowsize=0.5,density=3)
ax.quiver(Lat1[:,::-1],Z1,v1[::-1,::-1]*25,-omega1[::-1,::-1]*100,scale=150000,width=0.0025, headwidth=8, headlength=3, headaxislength=1)
ax.semilogy([90,90], [1000,100], color='C3')
ax.set_ylim(1000, 100)
ax.set_yticks(np.linspace(100, 1000, 10),[1000,900,800,700,600,500,400,300,200,100][::-1])
ax.set_xlim(-30, 60)
ax.set_xlabel('Latitude[°]')
ax.set_ylabel('level[hPa]')
ax.set_xticks(range(-30,61,10),['30°','20°','10°S', 'EQ','10°','20°', '30°', '40°','50°','60°N'])
ax.set_xlim([-30,60])

plt.title("1月")

#002
ax = plt.subplot(212)
#绘制阴影区
v2 = v2/np.mean(v2)
omega2 = omega2/np.mean(omega00)
ec2 = ax.contourf(Lat2[:,::-1],Z2, omega2[::-1,::-1]*5,levels=np.arange(round(omega2.min(),0),0,19.6),cmap = 'gray',extend='both')
cb2 = fig.colorbar(ec2,extend='both',label='hPa/s')
ax2 = cb2.ax
ax2.set_yticks(np.arange(round(omega2.min(),0),0,19.6),[-4, -3, -2, -1, 0])
levels = np.arange(-10,20,2)
#二维列表逆序
ax.streamplot(Lat2[:,::-1],Z2,v2[::-1,::-1],omega2[::-1,::-1]*10,color='k',linewidth=2,arrowsize=0.5,density=3)
ax.quiver(Lat2[:,::-1],Z2,v2[::-1,::-1]*150,-omega2[::-1,::-1]*170,scale=180000,width=0.0025, headwidth=8, headlength=3, headaxislength=1)
ax.semilogy([90,90], [1000,100], color='C3')
ax.set_ylim(1000, 100)
ax.set_yticks(np.linspace(100, 1000, 10),[1000,900,800,700,600,500,400,300,200,100][::-1])
ax.set_xlim(-30, 60)
ax.set_xlabel('Latitude[°]')
ax.set_ylabel('level[hPa]')
ax.set_xticks(range(-30,61,10),['30°','20°','10°S', 'EQ','10°','20°', '30°', '40°','50°','60°N'])
ax.set_xlim([-30,60])
plt.title("7月")
plt.suptitle("图4-1-10 75°—110°E的平均经圈环流",fontsize='xx-large')

plt.savefig("图4-1-10（黑白）.png",dpi=800)

plt.show()

 这两幅图像的配色与数据处理能有更好的方法。