Python气象绘图——风场

一、作业要求

1）绘制850 hPa 2020年夏季平均（6-7-8）的uv风场矢量图

2）绘制850 hPa 2020年夏季平均（6-7-8）的流线图

二、数据介绍

数据来源于欧洲中期天气预报中心（ECMWF）的再分析数据（ERA5），使用的是风场矢量数据，由于是风矢量，我们使用的数据以纬线为横轴，经线为纵轴分解为纬向风（u）和经向风（v），数据文件也是分开的。

在Pycharm中读取数据，并通过print函数展示数据信息。

可以看到数据是2020年的月平均风场，分辨率为1°×1°，垂直分层10、20、...、850、925、1000hPa，我们本次使用的是850hPa的。

ECMWF再分析数据的下载网址如下：

https://cds.climate.copernicus.eu/cdsapp#!/search?type=datasethttps://cds.climate.copernicus.eu/cdsapp#!/search?type=dataset

三、代码

1、引入需要的库（图像中文显示与负号“-”显示问题解决）

import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker

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

2、数据处理

filename1 = r'E:\python\pysx\ERA5.u_component_of_wind.1-1.mon.2020.nc'
filename2 = r'E:\python\pysx\ERA5.v_component_of_wind.1-1.mon.2020.nc'
f1 = xr.open_dataset(filename1)
#print(f1) # 查看数据信息
sumu_850 = f1.u.loc[f1.time.dt.month.isin([6,7,8])].loc['2020-01-01':'2020-12-01',850] # 读取2020-01-01至2020-12-01的
# 月份为6、7、8月（即夏季）的层次为850hPa的纬向风场
f2 = xr.open_dataset(filename2)
sumv_850 = f2.v.loc[f1.time.dt.month.isin([6,7,8])].loc['2020-01-01':'2020-12-01',850] # 同理，读取经向风场
aversumu_850 = (sumu_850.values[0,:,:]+sumu_850.values[1,:,:]+sumu_850.values[2,:,:])/3. # 求取夏季平均，sumu_850为DataFrame
# 格式数据包含行列名称等内容，计算仅需对风场数据进行计算sumu_850.values即风场数据的矩阵，不包含行列名称等。
aversumv_850 = (sumv_850.values[0,:,:]+sumv_850.values[1,:,:]+sumv_850.values[2,:,:])/3. 
lon = f1.longitude.values # 读取经度数据
lat = f1.latitude.values # 读取纬度数据

3、绘图

两张图分开绘制。

（1）矢量图

# 检验数据读取绘图
fig = plt.figure(figsize=(12, 6))
leftlon, rightlon, lowerlat, upperlat = (-180, 180, -90, 90)
proj = ccrs.PlateCarree(central_longitude=180)  # 指定投影为经纬度投影。
ax1 = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection=proj)  # 绘制第一个子图
ax1.set_extent([leftlon, rightlon, lowerlat, upperlat], crs=proj)
# 设置经纬度标签
ax1.set_xticks(np.arange(leftlon, rightlon + 30, 30), crs=proj)
ax1.set_yticks(np.arange(lowerlat, upperlat + 20, 20), crs=proj)
lon_formatter = cticker.LongitudeFormatter()
lat_formatter = cticker.LatitudeFormatter()
ax1.xaxis.set_major_formatter(lon_formatter)
ax1.yaxis.set_major_formatter(lat_formatter)
coastline_50m = cfeature.COASTLINE.with_scale('50m')
ax1.add_feature(coastline_50m)
c = ax1.quiver(lon[::5], lat[::5], aversumu_850[::5,::5], aversumv_850[::5,::5],
               color='navy',angles='xy',scale=200,width=0.002,transform=ccrs.PlateCarree())
# lon[::5], lat[::5], aversumu_850[::5,::5], aversumv_850[::5,::5]经纬度与风场都是格点数据，每个数据对应一个经纬坐标，
# 由于风场数据分辨率很高，绘制的风场过于密集，经纬度与风场数据都需要间隔5取用，才能使得图像清晰美观。
ax1.set_title('2020年夏季850hPa风场矢量图',fontsize=15,loc='center')
plt.quiverkey(c,X=0.98,Y=1.02,U=5,angle=90,labelpos='W',label='Wind:5m/s',color='r',labelcolor='r')
# 这句代码是设置矢量图的标识矢量，使用方法就不赘述了

（2）流线图

# 检验数据读取绘图
fig = plt.figure(figsize=(12, 6))
leftlon, rightlon, lowerlat, upperlat = (-180, 180, -90, 90)
proj = ccrs.PlateCarree(central_longitude=180)  # 指定投影为经纬度投影（等距圆柱投影）。
ax2 = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection=proj)  # 绘制第一个子图
ax2.set_extent([leftlon, rightlon, lowerlat, upperlat], crs=proj)
# 设置经纬度标签
ax2.set_xticks(np.arange(leftlon, rightlon + 30, 30), crs=proj)
ax2.set_yticks(np.arange(lowerlat, upperlat + 20, 20), crs=proj)
lon_formatter = cticker.LongitudeFormatter()
lat_formatter = cticker.LatitudeFormatter()
ax2.xaxis.set_major_formatter(lon_formatter)
ax2.yaxis.set_major_formatter(lat_formatter)
# 添加网格线
gl=ax2.gridlines(draw_labels=True,
                 #xlocs=np.arange(leftlon, rightlon + 40, 40),
                 #ylocs=np.arange(lowerlat, upperlat + 30, 30),
                 color='gray',alpha=0.8,linestyle="--",linewidth=1)
gl.right_labels=False
gl.top_labels=False
gl.left_labels=False
gl.bottom_labels=False
coastline_50m = cfeature.COASTLINE.with_scale('50m')
ax2.add_feature(coastline_50m)
c = ax2.streamplot(lon, lat, aversumu_850, aversumv_850,
               color='navy',density=3.0,linewidth=1.2,transform=ccrs.PlateCarree())
ax2.set_title('2020年夏季850hPa风场流线图',fontsize=15,loc='center')

plt.show()

四、绘图结果

图中可以清晰分辨出北太平洋副热带高压的反气旋式风场，以及东亚地区夏季风，南亚地区夏季风等夏季天气系统，由此可以检验图像绘制的正确性。

五、结语

受限于本人专业知识水平以及代码功底限制，多有错误，欢迎批评指正！