使用数据: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()
这两幅图像的配色与数据处理能有更好的方法。