天气学原理插图复现（三）——北半球500hPa平均西风风速分布

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

使用库：Matplotlib、NumPy、netCDF4、Cartopy、math

彩色图像

import matplotlib.pyplot as plt
import netCDF4 as nc
import numpy as np
import cartopy.mpl.ticker as cmt
import cartopy.crs as ccrs
import matplotlib.path as mpath
import cartopy.feature as cfeature
from cartopy.util import add_cyclic_point
from cartopy.mpl.ticker import (LongitudeFormatter, LatitudeFormatter)
import matplotlib as mpl
import matplotlib.ticker as mticker
import math

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

fu = nc.Dataset(r"uwnd.mon.ltm.1991-2020.nc",'r')
fv = nc.Dataset(r"vwnd.mon.ltm.1991-2020.nc",'r')

para_level = 500
levels = fu.variables['level'][:]
lat = fu.variables['lat'][:]
lon0 = fu.variables['lon'][:]
time = nc.num2date(fu.variables['time'][:],fu.variables['time'].units).data
months = np.array([i.month for i in time])
u = fu.variables['uwnd'][:]
v = fv.variables['vwnd'][:]


def uv(u,v):
    i=len(u)
    j=len(u[0])
    k=len(u[0][0])
    l=len(u[0][0][0])
    u=np.array(u)
    v=np.array(v)
    for a in range(i):
        u[a]=np.array(u[a])
        v[a]=np.array(v[a])
        for b in range(j):
            u[a][b]=np.array(u[a][b])
            v[a][b]=np.array(v[a][b])
            for c in range(k):
                u[a][b][c]=np.array(u[a][b][c])
                v[a][b][c]=np.array(v[a][b][c])
    uv=(u**2+v**2)**0.5
    wdir = 180.0 + np.arctan2(u, v)*180/np.pi
    for a in range(i):
        for b in range(j):
            for c in range(k):
                for d in range(l):
                    if 0 < wdir[a][b][c][d] < 180:
                        uv[a][b][c][d]=-uv[a][b][c][d]
    uv.tolist()
    for a in range(i):
        uv[a].tolist()
        for b in range(j):
            uv[a][b].tolist()
            for c in range(k):
                uv[a][b][c].tolist()
    return uv

def pa(para_month):
    index_level = np.where(levels==para_level)[0][0] # 改等压面
    index_month = np.where(months==para_month)[0][0] # 改月份        
    wspd = uv(u,v)
    wspd = wspd[index_month][index_level,:,:]
    wspd, lon = add_cyclic_point(wspd, coord=lon0)
    Lon, Lat = np.meshgrid(lon,lat)
    return Lon,Lat,wspd


Lon1,Lat1,wspd1=pa(1)
Lon2,Lat2,wspd2=pa(7)
fig = plt.figure(figsize=(16,6))

proj =ccrs.NorthPolarStereo(central_longitude=180)
leftlon, rightlon, lowerlat, upperlat = (-180,180,0,90)
img_extent = [leftlon, rightlon, lowerlat, upperlat]

#001
ax = plt.subplot(121,projection=proj)
ax.set_extent([-180,180,0,90], ccrs.PlateCarree())
theta = np.linspace(0, 2*np.pi, 100)
center, radius = [0.5, 0.5], 0.5
verts = np.vstack([np.sin(theta), np.cos(theta)]).T
circle = mpath.Path(verts * radius + center)
ax.set_boundary(circle, transform=ax.transAxes)
levels = np.arange(-20,45,5)
ec = ax.contourf(Lon1,Lat1, wspd1,cmap = 'rainbow',levels=levels,transform=ccrs.PlateCarree())
ac = ax.contour(Lon1,Lat1, wspd1,levels=levels,colors='black',linewidths=1,linestyles='-',transform=ccrs.PlateCarree())
ax.clabel(ac)
g1=ax.gridlines( draw_labels=True, linewidth=1, color='gray',linestyle='--')
g1.xlocator = mticker.FixedLocator(np.linspace(-180,180,13))
g1.ylocator = mticker.FixedLocator(np.linspace(0, 90,4))
g1.xlabel_style={'size':8}
g1.ylabel_style={'size':8}
ax.add_feature(cfeature.COASTLINE)
plt.title("1月")

#002
ax = plt.subplot(122,projection=proj)
ax.set_extent([-180,180,0,90], ccrs.PlateCarree())
theta = np.linspace(0, 2*np.pi, 100)
center, radius = [0.5, 0.5], 0.5
verts = np.vstack([np.sin(theta), np.cos(theta)]).T
circle = mpath.Path(verts * radius + center)
ax.set_boundary(circle, transform=ax.transAxes)
levels = np.arange(-20,45,5)
ec = ax.contourf(Lon2,Lat2, wspd2,cmap = 'rainbow',levels=levels,transform=ccrs.PlateCarree())
ac = ax.contour(Lon2,Lat2, wspd2,levels=levels,colors='black',linewidths=1,linestyles='-',transform=ccrs.PlateCarree())
fig.colorbar(ec)
ax.clabel(ac)
g1=ax.gridlines( draw_labels=True, linewidth=1, color='gray',linestyle='--')
g1.xlocator = mticker.FixedLocator(np.linspace(-180,180,13))
g1.ylocator = mticker.FixedLocator(np.linspace(0, 90,4))
g1.xlabel_style={'size':8}
g1.ylabel_style={'size':8}
ax.add_feature(cfeature.COASTLINE)
plt.title("7月")

plt.suptitle("图4-1-4   北半球500hPa平均西风风速分布（单位：m/s）",fontsize='xx-large')

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

黑白图像

 

import matplotlib.pyplot as plt
import netCDF4 as nc
import numpy as np
import cartopy.mpl.ticker as cmt
import cartopy.crs as ccrs
import matplotlib.path as mpath
import cartopy.feature as cfeature
from cartopy.util import add_cyclic_point
from cartopy.mpl.ticker import (LongitudeFormatter, LatitudeFormatter)
import matplotlib as mpl
import matplotlib.ticker as mticker
import math

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

fu = nc.Dataset(r"uwnd.mon.ltm.1991-2020.nc",'r')
fv = nc.Dataset(r"vwnd.mon.ltm.1991-2020.nc",'r')

para_level = 500
levels = fu.variables['level'][:]
lat = fu.variables['lat'][:]
lon0 = fu.variables['lon'][:]
time = nc.num2date(fu.variables['time'][:],fu.variables['time'].units).data
months = np.array([i.month for i in time])
u = fu.variables['uwnd'][:]
v = fv.variables['vwnd'][:]


def uv(u,v):
    i=len(u)
    j=len(u[0])
    k=len(u[0][0])
    l=len(u[0][0][0])
    u=np.array(u)
    v=np.array(v)
    for a in range(i):
        u[a]=np.array(u[a])
        v[a]=np.array(v[a])
        for b in range(j):
            u[a][b]=np.array(u[a][b])
            v[a][b]=np.array(v[a][b])
            for c in range(k):
                u[a][b][c]=np.array(u[a][b][c])
                v[a][b][c]=np.array(v[a][b][c])
    uv=(u**2+v**2)**0.5
    wdir = 180.0 + np.arctan2(u, v)*180/np.pi
    for a in range(i):
        for b in range(j):
            for c in range(k):
                for d in range(l):
                    if 0 < wdir[a][b][c][d] < 180:
                        uv[a][b][c][d]=-uv[a][b][c][d]
    uv.tolist()
    for a in range(i):
        uv[a].tolist()
        for b in range(j):
            uv[a][b].tolist()
            for c in range(k):
                uv[a][b][c].tolist()
    return uv

def pa(para_month):
    index_level = np.where(levels==para_level)[0][0]
    index_month = np.where(months==para_month)[0][0]   
    wspd = uv(u,v)
    wspd = wspd[index_month][index_level,:,:]
    wspd, lon = add_cyclic_point(wspd, coord=lon0)
    Lon, Lat = np.meshgrid(lon,lat)
    return Lon,Lat,wspd


Lon1,Lat1,wspd1=pa(1)
Lon2,Lat2,wspd2=pa(7)
fig = plt.figure(figsize=(16,6))

proj =ccrs.NorthPolarStereo(central_longitude=180)
leftlon, rightlon, lowerlat, upperlat = (-180,180,0,90)
img_extent = [leftlon, rightlon, lowerlat, upperlat]

#001
ax = plt.subplot(121,projection=proj)
ax.set_extent([-180,180,0,90], ccrs.PlateCarree())
theta = np.linspace(0, 2*np.pi, 100)
center, radius = [0.5, 0.5], 0.5
verts = np.vstack([np.sin(theta), np.cos(theta)]).T
circle = mpath.Path(verts * radius + center)
ax.set_boundary(circle, transform=ax.transAxes)
levels1 = np.arange(0,40,5)
ac1 = ax.contour(Lon1,Lat1, wspd1,levels=levels1,colors='black',linewidths=1,linestyles='-',transform=ccrs.PlateCarree())
levels2 = np.arange(-10,0,5)
ac2 = ax.contour(Lon2,Lat2, wspd1,levels=levels2,colors='black',linewidths=1,linestyles='dotted',transform=ccrs.PlateCarree())
ax.clabel(ac1)
ax.clabel(ac2)
g1=ax.gridlines( draw_labels=True, linewidth=1, color='gray',linestyle='--')
g1.xlocator = mticker.FixedLocator(np.linspace(-180,180,13))
g1.ylocator = mticker.FixedLocator(np.linspace(0, 90,4))
g1.xlabel_style={'size':8}
g1.ylabel_style={'size':8}
ax.add_feature(cfeature.COASTLINE)
plt.title("1月")

#002
ax = plt.subplot(122,projection=proj)
ax.set_extent([-180,180,0,90], ccrs.PlateCarree())
theta = np.linspace(0, 2*np.pi, 100)
center, radius = [0.5, 0.5], 0.5
verts = np.vstack([np.sin(theta), np.cos(theta)]).T
circle = mpath.Path(verts * radius + center)
ax.set_boundary(circle, transform=ax.transAxes)
levels1 = np.arange(0,40,5)
ac3 = ax.contour(Lon2,Lat2, wspd2,levels=levels1,colors='black',linewidths=1,linestyles='-',transform=ccrs.PlateCarree())
levels2 = np.arange(-10,0,5)
ac4 = ax.contour(Lon2,Lat2, wspd2,levels=levels2,colors='black',linewidths=1,linestyles='dotted',transform=ccrs.PlateCarree())
ax.clabel(ac3)
ax.clabel(ac4)
g1=ax.gridlines( draw_labels=True, linewidth=1, color='gray',linestyle='--')
g1.xlocator = mticker.FixedLocator(np.linspace(-180,180,13))
g1.ylocator = mticker.FixedLocator(np.linspace(0, 90,4))
g1.xlabel_style={'size':8}
g1.ylabel_style={'size':8}
ax.add_feature(cfeature.COASTLINE)
plt.title("7月")

plt.suptitle("图4-1-4   北半球500hPa平均西风风速分布（单位：m/s）",fontsize='xx-large')

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

 风速处理肯定有更好的方法，我写的太长了。