大气科学可视化绘图——利用python绘制hovmoller图教程（关于如何在静态图片中显示一些运动）

Hovmöller 图将所有值平均在一列经度或一行纬度中，并将这些值放在一个轴上。另一个轴表示时间。在上面的示例中，我们使用水平维度（X 轴）来表示时间，使用垂直维度（Y 轴）来表示纬度。通过这种方式，我们可以展示雨水在一年中如何向北或向南移动。

但是当阅读了很多文献的时候，我们常常发现，一般是以经度为X轴，时间为Y轴，按照纬度平均，按照时间顺序排列。

而在大气中使用hovmoller图的时候，我们一般有两种值，一种是选择降水的值进行绘制，一种是选择u和v的值进行绘制。这里我们举例的就拿最常见的降水来展示。

此处我们选取的是2004年6月19日至6月30日的滤波后的GPCP降水数据。

此处导入需要的模块

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.gridspec as gridspec
import metpy.calc as mpcalc
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from matplotlib.ticker import AutoMinorLocator, MultipleLocator, FuncFormatter
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
import matplotlib.pyplot as plt
import xarray as xr
import time
import numpy as np
import math
import xarray as xr

下一步读取我们需要的数据并且将需要计算的数据计算出来：

# 降水数据
predata=xr.open_dataset(path).sel(lat=slice(0,15),lon=slice(100,240),time=slice('2004-06-19','2004-06-30'))
pre=predata.pre.values
lon=predata.lon.values
lat=predata.lat.values
time=predata.time.values.astype('datetime64[ns]').astype("datetime64[D]")
#此处将时间格式从毫秒改为天，因为我们的值是daily，可以根据自己数据的情况进行修改
weights=np.cos(np.deg2rad(lat))
avg_data=(pre*weights[None,:,None]).sum(axis=1)/np.sum(weights)
# 计算权重并在整个纬度维度上取加权平均值

很简单，接下来就可以进行绘图准备：

fig = plt.figure(figsize=(12,18))#创建画布
gs = gridspec.GridSpec(nrows=2, ncols=1, height_ratios=[1, 5], hspace=0.05)#设置绘图的比例和图与图之间的距离

我们在这里先绘制我们计算的经纬度范围的地图，这样我们可以直观的在地图上对应看到一个移动的趋势：（此处设置了一个画地图的函数，在其他图中也可以这样用，会比较方便）

# 起始图
def make_map(ax, title):
    # set_extent  set crs
    ax.set_extent(box, crs=ccrs.PlateCarree())
    land = cfeature.NaturalEarthFeature('physical',
                                        'land',
                                        scale,
                                        edgecolor='grey', 
                                        facecolor='grey'
                                        ,zorder=2
                                        )
    ax.add_feature(land)  # set land color
    ax.coastlines(scale)  # set coastline resolution
    # set coordinate axis
    ax.set_xticks(np.arange(box[0],box[1]+xstep, xstep),crs=ccrs.PlateCarree())
    ax.set_yticks(np.arange(box[2], box[3]+ystep, ystep),crs=ccrs.PlateCarree())
    ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label =False))#经度0不加标识
    ax.yaxis.set_major_formatter(LatitudeFormatter())
    # plt.tick_params(labelsize=25)

    ax.set_title(title, fontsize=10, loc='left',pad=12)
    ax.yaxis.set_minor_locator(AutoMinorLocator(5))
    ax.xaxis.set_minor_locator(AutoMinorLocator(10))
    ax.tick_params(which='minor', 
                    direction='inout', length=4,width=0.59,
                    right=True, top=True)
    ax.tick_params(which='major', 
                    direction='out', 
                    length=8,
                    width=0.99, 
                    pad=3, 
                    labelsize=14,
                    bottom=True, left=True, right=True, top=True)
    return ax
# # prepare 
box = [100, 180+60, 0, 15]  
scale = '50m'            
xstep, ystep = 20, 5 
# cmap=plt.get_cmap('Reds')#'RdYlBu_r'
levels=np.arange(0,31)
# 地势最高的地块（制作小地图）
ax1 = fig.add_subplot(gs[0, 0], projection=ccrs.PlateCarree(central_longitude=180))

make_map(ax1,'map')

# 添加地缘边界以供地图参考
ax1.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax1.add_feature(cfeature.LAKES.with_scale('50m'), color='black', linewidths=0.5)
# 设置一些标题
plt.title('Hovmoller Diagram', loc='left', fontsize=14)
plt.title('GPCP Reanalysis', loc='right', fontsize=14)

接下来绘制我们的重头戏hovmoller图：

# Hovmoller图的底图
ax2 = fig.add_subplot(gs[1, 0])
ax2.invert_yaxis()  # 颠倒时间顺序，先做最早的
# 所选变量的图在纬度上取平均值并稍作平滑
clevs = np.arange(-4,4.8,0.8)
cf = ax2.contourf(lon, time, mpcalc.smooth_n_point(
    avg_data, 9, 2), clevs, cmap=plt.cm.RdBu_r, extend='both')
cs = ax2.contour(lon, time, mpcalc.smooth_n_point(
    avg_data, 9, 2), clevs, colors='k',linewidths=1.8)
cbar = plt.colorbar(cf, orientation='horizontal', pad=0.04, aspect=50, extendrect=True)
cbar.set_label('$mm/day$', fontsize=14)

# 做一些刻度线和刻度线标签
ax2.set_xticks(np.arange(box[0],box[1]+xstep, xstep),crs=ccrs.PlateCarree(),fontsize=16)
ax2.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label =False))
# ax2.set_xticks(np.arange(100,box[1]+xstep, xstep),crs=ccrs.PlateCarree(central_longitude=180), fontsize=14)
ax2.set_yticks(time, fontsize=16)
ax2.set_xlabel('longitude', fontsize=16)
ax2.set_ylabel('time', fontsize=16)
ax2.tick_params(axis='x',direction='in', length=6, width=1,colors='k')
y_major_locator=MultipleLocator(1)
# 设置一些标题
plt.title('850-hPa precipitation', loc='left', fontsize=16)
plt.title('Time Range: 2004-06-19 ——2004-06-30',loc='right', fontsize=16)

这样我们的图就画完啦！是不是非常简单

 还可以根据数据进行一定的自我设置哦，比如如果需要分析风速和降水的关系，我们可以在降水的底色图上加上风速的等值线，这样就可以进行一个对比。但是这里的av_data0一定要重新计算哦，他的经纬度和time都会因为数据的不同需要进行自定义的修改！

cs = ax2.contour(lon0, time0, mpcalc.smooth_n_point(
    avg_data0【换成想要的等值线的数据】, 9, 2), clevs, colors='k',linewidths=1.8)

初次见面，请多关照！希望能解决你的一点小烦恼哦！

一个在努力学习python的小菜鸟！

水平有限，欢迎指正！！！

欢迎评论、收藏、点赞、转发、关注。

关注我不后悔，记录学习进步的过程~~