用Python绘制《天气学原理和方法》插图
作者：Vector
最近天气学原理需要绘制课本插图来做 翻转课堂，因此整理了课本第四章几个典型图片的画法和代码，共需要的人使用。
数据
本文使用数据来自：https://psl.noaa.gov/data/gridded/data.ncep.reanalysis.derived.surface.html
图4.2 北半球冬季平均经向风
难点：各种库的下载与安装，主要难点在xarray上。
因为不需要设置地图，此图还是比较简单的。结果和代码如下。
fig4.2
代码语言：javascript复制import matplotlib.pyplot as plt
import xarray as xr
import numpy as np
import matplotlib.ticker as ticker
dir_loc = r"D:meteo_datancarvwnd.mon.mean.nc"
def DrawVWind():
# 读取数据
vwind = xr.open_dataset(dir_loc)["vwnd"]
# 选择使用的年份
year_label = (vwind['time'].dt.year >= 1990) & (vwind['time'].dt.year <= 2019)
# 选择年份做季节平均
vwindM = vwind[year_label].groupby("time.season").mean(dim="time")
# 选择DJF，做 纬向平均
vwindM2 = vwindM.loc[dict(season="DJF")].mean("lon")
# 选择合适位置的数据
DrawData = vwindM2.loc[:100, 90:-12.5]
print(DrawData)
Levels = DrawData["level"]
lat = DrawData['lat']
# 开始画图
fig, ax = plt.subplots(figsize=(15, 13))
# 彩色底图
MC = ax.contourf(lat, Levels, DrawData, cmap="RdBu_r", levels=np.arange(-3.5, 3.6, 0.5))
# 反转x,y轴
ax.invert_xaxis()
ax.invert_yaxis()
# 画等值线
newLevel = np.concatenate([np.copy(MC.levels), np.array([-0.2, 0.2])])
IsoLine = ax.contour(MC, levels=newLevel, colors="k", origin="lower")
# 在线上标字
ax.clabel(IsoLine, fmt='%2.1f', fontsize=12)
# 设置横纵坐标
ax.set_ylabel("Press(hPa)", fontsize=15)
ax.set_xlabel("Latitude(°)", fontsize=15)
ax.xaxis.set_major_locator(ticker.MultipleLocator(10))
# ax.set_yscale("symlog")
ax.yaxis.set_major_locator(ticker.MultipleLocator(200))
ax.tick_params(labelsize=15)
# 设置colorbar
cbar = fig.colorbar(MC, orientation='horizontal',shrink=0.8)
cbar.set_label("DJF vwind over NH(m/s)", fontsize=13)
plt.savefig("fig4.2.png", dpi=300)
plt.show()
if __name__ == '__main__':
DrawVWind()
图4.5 500hPa平均等高线
难点：地图投影、如何把地图变圆、cartopy下载
cartopy下载可能遇到很多问题，慢慢来就好。
参考cartopy官网，设置原形边框可以：
代码语言：javascript复制 # 画圆所需圆框
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)
使用的投影为NorthPolarStereo
最后的结果为
完整代码：
代码语言：javascript复制import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
import matplotlib.path as mpath
import cartopy.feature
from cartopy.util import add_cyclic_point
from cartopy.mpl.ticker import (LongitudeFormatter, LatitudeFormatter)
fLoc = r"D:meteo_datancarhgt.mon.mean.nc"
def Draw500HgtJulyG():
# 打开文件
Hgt = xr.open_dataset(fLoc)["hgt"]
# 选定数据
time_label = (Hgt['time'].dt.year >= 1990) & (Hgt['time'].dt.year <= 2019) &
(Hgt["time"].dt.month == 7)
Hgt500Mean = Hgt[time_label].loc[dict(level=500.0)].mean("time").loc[:0, :]
# 设置画布
fig = plt.figure(figsize=(16, 10))
# 选定经纬度
lat = np.array(Hgt500Mean['lat'])
lon = np.array(Hgt500Mean["lon"])
# 创建一个投影
proj = ccrs.NorthPolarStereo(central_longitude=90)
# 创建一个画纸， 并指明投影类型
ax = plt.axes(projection=proj)
ax.coastlines() # 画海岸线
ax.add_feature(cartopy.feature.LAND)
# 画圆所需圆框
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)
# 经纬网格
gl = ax.gridlines(draw_labels=True, color='gray', alpha=0.5, linestyle='--',
xlocs=np.arange(-180, 181, 30), ylocs=[],
x_inline=False)
gl.xformatter = LongitudeFormatter()
gl.yformatter = LatitudeFormatter()
# 消除环形白带
cycleHgt500, cycle_lon = add_cyclic_point(np.array(Hgt500Mean), coord=lon)
cycle_LON, cycle_LAT = np.meshgrid(cycle_lon, lat)
# 彩色底图
m0 = ax.contourf(cycle_LON, cycle_LAT, cycleHgt500 / 10, levels=np.arange(540, 609, 8)
, cmap='RdBu',
transform=ccrs.PlateCarree()) # 注意选择投影
# 等值线
m1 = ax.contour(m0, levels=m0.levels, colors="k", transform=ccrs.PlateCarree())
# 等值线上标字
ax.clabel(m1, fmt='%2.0f', fontsize=13, colors="k")
# 设置colorbar
cbar = fig.colorbar(m0, shrink=0.8)
cbar.set_label("500 hPa geo-height over NH in July(10m)", fontsize=13)
plt.savefig("fig4.5.png", dpi=300)
plt.show()
if __name__ == '__main__':
Draw500HgtJulyG()
图4.7 7月平均海平面气压
难点：地图投影
主要难在地图投影上和消除白线上。消除白线可以参考：气象数据处理--Cartopy - Vector的文章 - 知乎 https://zhuanlan.zhihu.com/p/347944786
地图投影可以选择Miller
fig4.7
这里其实有一个问题，格陵兰岛根据NOAA的卫星资料在夏季是低压，可能资料有一点问题。
完整代码：
代码语言：javascript复制import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
import matplotlib.path as mpath
import cartopy.feature
from cartopy.util import add_cyclic_point
from cartopy.mpl.ticker import (LongitudeFormatter, LatitudeFormatter)
Floc = r"D:meteo_datancarSeaLevelP.mon.mean.nc"
def DrawSeaLPJuly():
Press = xr.open_dataset(Floc)["slp"]
time_label = (Press['time'].dt.year >= 1980) & (Press['time'].dt.year <= 2019) &
(Press["time"].dt.month == 7)
PressT = Press[time_label].mean(dim="time").loc[90:-50, :]
lat = PressT["lat"]
lon = PressT["lon"]
# 消除中心白线
cyclePressT, cycle_lon = add_cyclic_point(np.array(PressT), coord=lon)
cycle_LON, cycle_LAT = np.meshgrid(cycle_lon, lat)
# 开始画图
fig = plt.figure(figsize=(15, 7))
proj = ccrs.Miller(central_longitude=120)
ax = plt.axes(projection=proj)
ax.coastlines() # 画海岸线
ax.add_feature(cartopy.feature.LAND)
m0 = ax.contourf(cycle_LON, cycle_LAT, cyclePressT
, cmap='RdBu',
transform=ccrs.PlateCarree(),
levels=np.arange(995, 1031, 5)) # 注意选择投影
ax.set_extent([-180, 179.999, 90, -50], crs=ccrs.PlateCarree())
m1 = ax.contour(m0, transform=ccrs.PlateCarree(), colors="k")
ax.clabel(m1, fmt="%2.0f", fontsize=10, colors="k")
gl = ax.gridlines(draw_labels=True, color='gray', alpha=0.5, linestyle=':',
x_inline=False)
gl.xformatter = LongitudeFormatter()
gl.yformatter = LatitudeFormatter()
cbar = fig.colorbar(m0)
cbar.set_label("Mean sea level pressure in July(hPa)", fontsize=13)
print(PressT)
plt.savefig("fig4.7.png",dpi=300)
plt.show()
if __name__ == '__main__':
DrawSeaLPJuly()
图4.23 7月低纬地区流线
难点：streamplot应用及各种因系统产生的bug
流线需要使用streamplot，官方文档https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.streamplot.html。简单应用参考：python电场线画法 - Vector的文章 - 知乎 https://zhuanlan.zhihu.com/p/222450226
因为使用streamplot的时候会用到scipy进行插值，如果你的windows user‘s name是中文会报错，可以参考https://blog.csdn.net/wls666/article/details/103334152给出的步骤将你的用户名改为英文。
fig4.23
完整代码：
代码语言：javascript复制import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
import matplotlib.path as mpath
import cartopy.feature
from cartopy.util import add_cyclic_point
from cartopy.mpl.ticker import (LongitudeFormatter, LatitudeFormatter)
u_loc = r'D:meteo_datancarsuruwnd.mon.mean.nc'
v_loc = r"D:meteo_datancarsurvwnd.mon.mean.nc"
Floc = r"D:meteo_datancarSeaLevelP.mon.mean.nc"
def DrawStreamLine():
# 打开文件
Press = xr.open_dataset(Floc)["slp"]
Uwind = xr.open_dataset(u_loc)["uwnd"]
Vwind = xr.open_dataset(v_loc)["vwnd"]
# 设置时间
time_label = (Uwind['time'].dt.year >= 1980) & (Uwind['time'].dt.year <= 2019) &
(Uwind["time"].dt.month == 7)
# 选定合适的位置
UwindT = Uwind[time_label].mean(dim="time").loc[60:-60, :]
VwindT = Vwind[time_label].mean(dim="time").loc[60:-60, :]
PressT = Press[time_label].mean(dim="time").loc[60:-60, :]
# 设置画布
fig = plt.figure()
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=-140))
# 选定范围
ax.set_extent([-180.0, 179.99, -60.0, 60.0], crs=ccrs.PlateCarree())
ax.coastlines() # 画海岸线
# 设置横纵坐标与补齐数据
lon, lat = UwindT['lon'], UwindT['lat']
cycleUwindT, cycle_lon = add_cyclic_point(np.array(UwindT), coord=lon)
cycleVwindT, _ = add_cyclic_point(np.array(VwindT), coord=lon)
cyclePressT, _ = add_cyclic_point(np.array(PressT), coord=lon)
cycle_LON, cycle_LAT = np.meshgrid(cycle_lon, lat)
# 画流线
ax.streamplot(cycle_LON, cycle_LAT, cycleUwindT, cycleVwindT, transform=ccrs.PlateCarree(), color="k",
linewidth=0.8)
# 画海平面气压场
m0 = ax.contourf(cycle_LON, cycle_LAT, cyclePressT
, cmap='RdBu',
transform=ccrs.PlateCarree())
# 经纬线
gl = ax.gridlines(draw_labels=True, color='gray', alpha=0.5, linestyle=':')
gl.xformatter = LongitudeFormatter()
gl.yformatter = LatitudeFormatter()
# 设置colorbar
cbar = fig.colorbar(m0, orientation='horizontal', shrink=0.8)
cbar.set_label("Mean sea level pressure in July(hPa)")
plt.savefig("fig4.23.png", dpi=300)
plt.show()
if __name__ == '__main__':
DrawStreamLine()
