天气学原理插图复现（五）——1000hPa平均风场

使用数据：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
from matplotlib.colors import ListedColormap

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

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

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] # 改月份
    u0 = u[index_month][index_level,:,:]
    v0 = v[index_month][index_level,:,:]
    u00, lon = add_cyclic_point(u0, coord=lon0)
    v00, lon = add_cyclic_point(v0, coord=lon0)
    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,u00,v00,wspd

Lon1,Lat1,u1,v1,wspd1=pa(1)
Lon2,Lat2,u2,v2,wspd2=pa(7)


fig = plt.figure(figsize=[12,10])
#001
ax = plt.subplot(211,projection=ccrs.PlateCarree(central_longitude=180))
ax.streamplot(Lon1,Lat1,u1,v1,transform=ccrs.PlateCarree(),color='k',linewidth=0.75,arrowsize=0.75,density=2)
levels = np.arange(0,24,2)
ac=ax.contourf(Lon1,Lat1, -wspd1,cmap = 'jet',levels=levels,transform=ccrs.PlateCarree())
ax.contourf(Lon1,Lat1, wspd1,cmap = 'jet',levels=levels,transform=ccrs.PlateCarree())
fig.colorbar(ac)
ax.set_extent([-180,180,-60,60])
ax.set_xticks([-180,-150,-120,-90,-60,-30,0,30,60,90,120,150,180])
ax.set_yticks([-60,-30,0,30,60])
ax.xaxis.set_major_formatter(LongitudeFormatter())
ax.yaxis.set_major_formatter(LatitudeFormatter()) 
ax.add_feature(cfeature.LAND,color='white')
ax.add_feature(cfeature.COASTLINE)
ax.tick_params(axis='both',which='major',labelsize=10,direction='out',length=10,width=0.8)
ax.minorticks_on()
ax.tick_params(axis='both',which='minor',direction='out',width=0.8,top=True,right=True)
plt.title("1月")

#002
ax = plt.subplot(212,projection=ccrs.PlateCarree(central_longitude=180))
ax.streamplot(Lon2,Lat2,u2,v2,transform=ccrs.PlateCarree(),color='k',linewidth=0.75,arrowsize=0.75,density=2)
levels = np.arange(0,24,2)
ac=ax.contourf(Lon2,Lat2, -wspd2,cmap = 'jet',levels=levels,transform=ccrs.PlateCarree())
ax.contourf(Lon2,Lat2, wspd2,cmap = 'jet',levels=levels,transform=ccrs.PlateCarree())
fig.colorbar(ac)
ax.set_extent([-180,180,-60,60])
ax.set_xticks([-180,-150,-120,-90,-60,-30,0,30,60,90,120,150,180])#指定要显示的经纬度 
ax.set_yticks([-60,-30,0,30,60])
ax.xaxis.set_major_formatter(LongitudeFormatter())#刻度格式转换为经纬度样式 
ax.yaxis.set_major_formatter(LatitudeFormatter()) 
ax.add_feature(cfeature.LAND,color='white')
ax.add_feature(cfeature.COASTLINE)
ax.tick_params(axis='both',which='major',labelsize=10,direction='out',length=10,width=0.8)
ax.minorticks_on()
ax.tick_params(axis='both',which='minor',direction='out',width=0.8,top=True,right=True)
plt.title("7月")

plt.suptitle("图4-1-6   1000hPa平均风场",fontsize='xx-large')

plt.savefig("图4-1-6.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
from matplotlib.colors import ListedColormap

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

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

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] # 改月份
    u0 = u[index_month][index_level,:,:]
    v0 = v[index_month][index_level,:,:]
    u00, lon = add_cyclic_point(u0, coord=lon0)
    v00, lon = add_cyclic_point(v0, coord=lon0)
    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,u00,v00,wspd

Lon1,Lat1,u1,v1,wspd1=pa(1)
Lon2,Lat2,u2,v2,wspd2=pa(7)


fig = plt.figure(figsize=[12,10])
#001
ax = plt.subplot(211,projection=ccrs.PlateCarree(central_longitude=180))
ax.streamplot(Lon1,Lat1,u1,v1,transform=ccrs.PlateCarree(),color='k',linewidth=0.75,arrowsize=0.75,density=2)
## 绘制阴影区
levels = np.arange(0,60,30)
ax.contourf(Lon1,Lat1, -wspd1,cmap = 'gray_r',levels=levels,transform=ccrs.PlateCarree())
ax.set_extent([-180,180,-60,60])
ax.set_xticks([-180,-150,-120,-90,-60,-30,0,30,60,90,120,150,180])
ax.set_yticks([-60,-30,0,30,60])
ax.xaxis.set_major_formatter(LongitudeFormatter()) 
ax.yaxis.set_major_formatter(LatitudeFormatter()) 
ax.add_feature(cfeature.LAND,color='white')
ax.add_feature(cfeature.COASTLINE)
ax.tick_params(axis='both',which='major',labelsize=10,direction='out',length=10,width=0.8)
ax.minorticks_on()
ax.tick_params(axis='both',which='minor',direction='out',width=0.8,top=True,right=True)
plt.title("1月")
#002
ax = plt.subplot(212,projection=ccrs.PlateCarree(central_longitude=180))
ax.streamplot(Lon2,Lat2,u2,v2,transform=ccrs.PlateCarree(),color='k',linewidth=0.75,arrowsize=0.75,density=2)
## 绘制阴影区
levels = np.arange(0,60,30)
ax.contourf(Lon2,Lat2, -wspd2,cmap = 'gray_r',levels=levels,transform=ccrs.PlateCarree())
ax.set_extent([-180,180,-60,60])
ax.set_xticks([-180,-150,-120,-90,-60,-30,0,30,60,90,120,150,180])
ax.set_yticks([-60,-30,0,30,60])
ax.xaxis.set_major_formatter(LongitudeFormatter())
ax.yaxis.set_major_formatter(LatitudeFormatter()) 
ax.add_feature(cfeature.LAND,color='white')
ax.add_feature(cfeature.COASTLINE)
ax.tick_params(axis='both',which='major',labelsize=10,direction='out',length=10,width=0.8)
ax.minorticks_on()
ax.tick_params(axis='both',which='minor',direction='out',width=0.8,top=True,right=True)
plt.title("7月")

plt.suptitle("图4-1-6   1000hPa平均风场",fontsize='xx-large')

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