Python气象绘图——垂直剖面图

1.画布创建

首先创建的画布不需要设置投影类型，使用默认的画布即可。

fig=plt.figure(figsize=(w,h),dpi=dpi)
axe=plt.subplot(1,1,1)

2.设置图题

axe.set_title(str, fontsize=size,y=y)

3.设置绘图边界

axe.set_xlim(l, r)
axe.set_ylim(b, t)
axe.invert_yaxis()  # 翻转纵坐标

这里需要注意两点：

a.默认画布的范围设置函数有所不同

b.由于我们所做垂直剖面图，一般纵坐标都是气压（离地高度也能够添加，后面做介绍）。而近地面的气压大，高空气压小，因此纵坐标由下至上的变化过程是由大到小。而set_ylim函数要求t>b，因此我们需要翻转纵坐标

4.网格绘制

axe.grid(color=grid_color, linestyle=grid_type, linewidth=grid_linewidth, axis=axis)

不是投影画布了，因此也不需要gridlines啦，可以省点力气

这里要提一下axis，可以设置的参数有'both'， 'x'，'y'。如果不设置axis的话默认就是both，x与y分别表示只显示x轴或者y轴的网格。

5.坐标轴设置

首先设置横坐标：

plt.xticks(position_list, name_list, fontsize=ticks_size, color=ticks_color)
axe.set_xlabel(str, fontproperties=Simsun, fontsize=ylabelsize)
axe.get_xaxis().set_visible('True')

plt.xticks函数用于自定义横坐标的刻度，position_list表示插入自定义刻度的位置；name_list与position_list对应，表示对应位置下插入的自定义刻度，string与float均可以。

set_xlabel设置的是横坐标的标签

axe.get_xaxis().set_visible('True')是获取了axe的横坐标刻度+标签，并设置可视化，True表示显示，False表示隐藏

然后设置纵坐标：

plt.yticks(fontsize=ticks_size, color=ticks_color)
axe.get_yaxis().set_visible('True')
axe.set_yticks(np.arange(b, t, i))
axe.set_ylabel(str, fontproperties=Simsun, fontsize=ylabelsize)

axe.set_yticks是另一种设置坐标刻度的方式，不支持自定义方式，好处就是用起来简单快捷

最后对刻度进行统一格式修改：

axe.tick_params(length=tick_length)

axe.tick_params有很多功能，这里只用到了设置刻度线的长短

labels = axe.get_xticklabels() + axe.get_yticklabels() 
[label.set_fontproperties(FontProperties(fname="./font/Times.ttf", size=8)) for label in labels]

然后是将刻度的字体进行修改一下，修改为Time New Roman字体，同时修改字体大小

6.气象信息绘制

把画布绘制完以后，就可以进行气象信息的绘制啦。由于是垂直剖面，因此在数据处理上需要花一些功夫，同时在绘制的时候也要注意一些技巧。这一部分干讲会比较难懂，先上一个图，然后逐步分析。

6.1数据处理

图中填色的内容是臭氧浓度，箭头表示剖切面的风速，右上角的红色线条表示剖切面。因此我们首先需要读取数据，然后对数据进行插值与处理，最后才是根据数据绘制图形。

在本文中以wrf模式的输出结果为例，进行示范。如果是同样使用wrf模式的，可以参考下面的数据处理过程。如果使用的不是wrf模式，可以跳过直接看后面的数据绘制过程。（由于只是为了绘图展示用，因此请不要过度在意图中数值的准确与否）

• 首先声明一下需要引入的库：

import netCDF4 as nc
from wrf import to_np,getvar,CoordPair,vertcross,interplevel

• 然后确定插值的起始经纬度与压力范围：

big_p,small_p,interval_p=1000,700,50
start_lon,end_lon,interval_lon=121,122,0.2
start_lat,end_lat=30.8,31.2

• 然后需要用到CoordPair库，将始末点包装为CoordPair对象，用于插值用：

startpoint = CoordPair(lat=start_lat, lon=start_lon)
endpoint = CoordPair(lat=end_lat, lon=end_lon)

• 然后照例读取wrf输出文件：

ncfile=nc.Dataset('D:\wrf_simulation\\2meic\\wrfout_d03_2016-07-21_12_2meic')

• 根据自己的需要读取数据：

o3=getvar(ncfile,'o3',timeidx=126)
o3 = o3 * 1000 / 22.4 * 48 * 273.15 / (getvar(ncfile, 'tk', timeidx=126)) * (getvar(ncfile, 'pressure', timeidx=126)) / 1013.25
lat=getvar(ncfile,'lat')
lon=getvar(ncfile,'lon')
height=getvar(ncfile,'height')
hgt=getvar(ncfile,'HGT')
height2earth=height-hgt
ua=getvar(ncfile,'ua',timeidx=126)
va=getvar(ncfile,'va',timeidx=126)
wa=getvar(ncfile,'wa',timeidx=126)
wa = wa * 10
p = getvar(ncfile, 'pressure', timeidx=126)

由于wrf-chem中o3的浓度单位是ppmv，需要将其转换为ug/m3的单位。而后分别读取纬度、经度；读取海拔高度与地形高度用于计算离地高度（之后有用）；读取风速数据，其中由于垂直风向很小，因此放大10倍用于观测方便；最后读取气压

• 读取完数值以后就开始进行插值了：

o3_vert = vertcross(o3, p, wrfin=ncfile, start_point=startpoint, end_point=endpoint, latlon=True)
ua_vert = vertcross(ua, p, wrfin=ncfile, start_point=startpoint, end_point=endpoint, latlon=True)
va_vert = vertcross(va, p, wrfin=ncfile, start_point=startpoint, end_point=endpoint, latlon=True)
wa_vert = vertcross(wa, p, wrfin=ncfile, start_point=startpoint, end_point=endpoint, latlon=True)

这里使用到的是vertcross函数，首先输入读取的数值与压力，然后由于本次插值输入的始末点是CoordPair对象的经纬度坐标，因此需要告诉vertcross指定文件是之前读取的ncfile，然后是始末点，最后的latlon=True可以使得最后的插值结果更详细一些，False的话结果更简单

输出的数据结果如下图，是一个xarray.DataArray对象。除了数据本身以外还包含坐标信息，而这个坐标包含两维，xy_loc表示横向的数据，即沿着剖切线方向的各个点的经纬度；vertical表示垂直方向的坐标，单位是hPa，由大到小。由于数据本身并不规整，并且xarray数据本身插值即存在一定不严格的情况，导致插值数据的起始经纬度与用户设置的会有一些出入，属于正常现象。只要初始经纬度都在数据所包括的范围内，就不会有太大问题。

• 由于风速具有u和v两个分量的，而剖切线方向并不一定会垂直于经纬度方向，所以需要将水平方向的风速投影到剖切线上：

#计算风速
ws_vert = np.sqrt(ua_vert ** 2 + va_vert ** 2)
#计算风向夹角
wdir_vert = np.arctan2(va_vert, ua_vert) * 180 / np.pi
line_angel = np.arctan2(end_lat - start_lat, end_lon - start_lon) * 180 / np.pi
vl_angel = wdir_vert - line_angel
#vl_angel = np.cos(vl_angel / 180 * np.pi)
vl_angel = np.cos(np.deg2rad(vl_angel))
ws_vert = ws_vert * vl_angel

关于这几个参数的关系，可以看下面这个简单的示意图。

最终要满足的就是让风速*vl_angel从而投影到剖切线方向。ws_vert是一个正值，不需要在这里考虑风速的方向，因为最后在剖切线上的正负值是由vl_angel来负责的。需要注意的是在使用np.cos的时候要将弧度转换为角度。而np.arctan2是个很好用的函数，可以自动输出360°的arctan值。

• 之前说了，由于插值的时候插值数据的网格点并不会规整分布，因此如果要精确绘图的话需要读取插值点的经纬度以及气压数据：

首先读取一下插值点的气压数据，这个不用处理，很简单。由于插值的臭氧，风速始末点都一致，因此每个数据点的经纬度与气压都是相同的，随便找一个插值数据读取即可。to_np是让数据从xarray转换为ndarray方便后续操作的，va_vert.coords可以直接读取xarray对象的Coordinates属性参数。

plist = to_np(va_vert.coords['vertical'])

然后根据气压值，去计算对应的离地高度，为了后续绘图用。这里用到的是interplevel，wrf-python的垂直方向插值函数，计算对应气压的离地高度。这里有一点需要注意，因为绘图的时候气压还是按照初始设定的气压起始以及间隔绘制的，因此高度也要对应的插值。

同时由于每个格点的离地高度多少会有点区别，但是区别又不大，因此采用了np.max函数取同一层下的最大值。最后用了个简单地方法快速取个整，如果要保留小数点的话建议用np.around函数

hlist = []
for i in range(big_p, small_p - interval_p, -interval_p):
    hlist.append(float(np.max(interplevel(height2earth, p, i)).values))
hlist = np.array([int(i) for i in hlist]) #简单地取整

接下去要读取插值数据的经纬度了，由于插值数据的xy_loc是一个难以直接读取的数据结构，因此干脆直接转换成字符串形式读取一下：

for i in range(len(ua_vert.coords['xy_loc'])):
     s = str(ua_vert.coords['xy_loc'][i].values)
     print(s)

他的经纬度数据还是有规律可循的，基本上每一个lat和lon后面就分别跟着经纬度，于是用字符串的find函数就能够解决，然后要将读取的字符串再变回float：

lonlist, latlist= [], []
for i in range(len(ua_vert.coords['xy_loc'])):
    s = str(ua_vert.coords['xy_loc'][i].values)
    lonlist.append(float(s[s.find('lon=') + 4:s.find('lon=') + 12]))
    latlist.append(float(s[s.find('lat=') + 4:s.find('lat=') + 12]))

最后处理我们保存的经纬度列表，以及之后要插入刻度的字符串：

str_lonlist, float_lonlist = [], []
a = np.mgrid[0:len(lonlist) - 1:complex(str(int(round((end_lon + interval_lon - start_lon) / interval_lon))) + 'j')] #生成位置列表
a = np.around(a, decimals=0) #四舍五入取整
#处理插入x轴的刻度标签
for i in range(int((end_lon + interval_lon - start_lon) / interval_lon)):
    float_lonlist.append(lonlist[int(a[i])])
    lo, la = round(lonlist[int(a[i])], 2), round(latlist[int(a[i])], 2)
    str_lonlist.append(str(lo) + '°E' + "\n" + str(la) + '°N')

由于我们的刻度肯定不是每个格点都插入的，要按照我们预期的初始经纬度和间隔插入，因此需要一个位置列表：a。而str_lonlist是用于保存插入的刻度的，做了个简单的处理让刻度可以显示两行，上面经度下面纬度。别的不多做解释了，都是python的基础处理方式。

6.2数据绘制

呼~数据终于处理好了

图像绘制就简单很多了：

o3_level=np.arange(40,410,20)
contourf = axe.contourf(lonlist, plist, o3_vert, levels=o3_level, cmap=cmaps.NCV_jaisnd,extend='neither')
axe.quiver(lonlist, plist, ws_vert, wa_vert, pivot='mid',
           width=0.001, scale=120, color='k', headwidth=4,
           alpha=1)

关于quiver和contourf可以看下面两篇文章，不做赘述了：爱垂钓的猫：Python气象绘图笔记（四）——填色与colorbar爱垂钓的猫：Python气象绘图笔记（六）——矢量箭头

6.3colorbar的绘制

同样可以参考这篇文章爱垂钓的猫：Python气象绘图笔记（四）——填色与colorbar

这里就贴上代码：

fig.subplots_adjust(right=0.78)
rect=[0.84,0.07,0.01,0.57]# 分别代表，水平位置，垂直位置，水平宽度，垂直宽度
cbar_ax=fig.add_axes(rect)
cb=fig.colorbar(contourf,drawedges=True,ticks=o3_level,cax=cbar_ax,orientation='vertical',spacing='uniform')
cb.set_label('$\mathrm{O_3}$浓度',fontsize=12)
cb.ax.tick_params(length=0)
labels=cb.ax.get_xticklabels()+cb.ax.get_yticklabels()
[label.set_fontproperties(FontProperties(fname="./font/Times.ttf",size=10))forlabelinlabels]

6.4小地图的绘制

其实右上角的小地图不是什么很难的东西，无非就是在画布上创建一个新的子图：

axe1 = fig.add_axes([0.77,0.69,0.2,0.2], projection=ccrs.PlateCarree())

不过这个子图我们给定了x的位置，y的位置，x方向的长度与y方向的长度。同时由于是地理绘图因此加上了投影类型。

下面就加上地理图的细节部分，贴上代码：

axe1.add_feature(cfeat.COASTLINE.with_scale("10m"), linewidth=0.7, color='k')
axe1.add_feature(cfeat.RIVERS.with_scale("10m"))
axe1.add_feature(cfeat.OCEAN.with_scale("10m"))
axe1.add_feature(cfeat.LAKES.with_scale("10m"))
axe1.set_extent([120, 123, 29.5, 32.5], crs=ccrs.PlateCarree())

gl = axe1.gridlines(crs=ccrs.PlateCarree(), draw_labels=False, linewidth=0.5, color='gray',linestyle=':')
gl.xlocator = mticker.FixedLocator(np.arange(120, 123.1, 1))
gl.ylocator = mticker.FixedLocator(np.arange(29.5, 32.6, 1))
axe1.set_xticks(np.arange(120, 123.1, 1), crs=ccrs.PlateCarree())
axe1.set_yticks(np.arange(29.5, 32.6, 1), crs=ccrs.PlateCarree())
plt.tick_params(top=True,bottom=True,left=True,right=True) #tick刻度显示
plt.tick_params(labeltop=True,labelleft=True,labelright=True,labelbottom=True) #tick标签显示
axe1.xaxis.set_major_formatter(LongitudeFormatter())
axe1.yaxis.set_major_formatter(LatitudeFormatter())
labels = axe1.get_xticklabels() + axe1.get_yticklabels()
[label.set_fontproperties(FontProperties(fname="./font/Times.ttf",size=6)) for label in labels]

关于代码的细节部分可以参考这篇文章，不多赘述了：

https://zhuanlan.zhihu.com/p/372821243爱垂钓的猫：Python气象绘图笔记（三）——地理图

最后在小地图上绘制剖切线，剖切点以及字体：

axe1.plot([lonlist[0], lonlist[-1]], [latlist[0], latlist[-1]], color='red', linewidth=1.5,linestyle='-')
axe1.plot(lonlist[0], latlist[0], marker='o', color='red', markersize=2.5)
axe1.plot(lonlist[-1], latlist[-1], marker='o', color='red', markersize=2.5)
axe1.text(lonlist[0], latlist[0], 'A', fontproperties=FontProperties(fname="./font/Times.ttf", size=8))
axe1.text(lonlist[-1], latlist[-1], 'A'+"'", fontproperties=FontProperties(fname="./font/Times.ttf", size=8))

plot与text是matplotlib的基础函数，参考资料很多也有很多写的很好的，我就不多说了，有疑问可以问我。

7.翻转纵坐标轴绘制

可能有人需要在主图的右侧绘制另一个纵坐标轴，这样可以标识离地高度。如下图所示：

那么首先，需要创建一个翻转y轴，共享x轴：

axe_sub = axe.twinx()

然后设置翻转y轴的范围：

axe_sub.set_ylim(hlist[0],hlist[-1])

接着，设置y轴添加ticks的位置。这一步是必需的，因为离地高度与气压的关系并不是线性的，但是对matplotlib直接添加我们事先准备好的hlist的话，会导致刻度错位。所以要先创建一个刻度的位置：

axe_sub.set_yticks(np.mgrid[hlist[0]:hlist[-1]:complex(str(int(hlist.shape[0])) + 'j')])

这个时候刻度的数值是这样的：

然后我们要将刻度位置的数值替换为之前计算好的数值：

round_hlist=np.around(hlist,-1)
axe_sub.set_yticklabels(round_hlist)

最后和之前一样，添加一下y轴标签，修改下参数：

axe_sub.set_ylabel('离地高度$\mathrm{(m))}$',fontsize=12)
axe_sub.tick_params(labelcolor='black',length=3)
labels=axe_sub.get_xticklabels()+axe_sub.get_yticklabels()
[label.set_fontproperties(FontProperties(fname="./font/Times.ttf",size=8))forlabelinlabels]

终于大功告成啦！最后贴上完整图和全部代码：

importnetCDF4asnc
fromwrfimportto_np,getvar,CoordPair,vertcross,interplevel
importcartopy.crsasccrs
importcartopy.featureascfeat
importmatplotlib.pyplotasplt
importcmaps
importmatplotlibasmpl
importmatplotlib.tickerasmticker
fromcartopy.mpl.tickerimportLongitudeFormatter,LatitudeFormatter
frommatplotlib.font_managerimportFontProperties
importnumpyasnp

Simsun=FontProperties(fname="./font/SimSun.ttf")
Times=FontProperties(fname="./font/Times.ttf")
mpl.rcParams['axes.unicode_minus']=False
config={
"font.family":'serif',
"mathtext.fontset":'stix',
"font.serif":['SimSun'],
}
mpl.rcParams.update(config)

#初始数据
big_p,small_p,interval_p=1000,700,50
start_lon,end_lon,interval_lon=121,122,0.2
start_lat,end_lat=30.8,31.2
#插值初始点与末点的坐标数据
startpoint=CoordPair(lat=start_lat,lon=start_lon)
endpoint=CoordPair(lat=end_lat,lon=end_lon)
ncfile=nc.Dataset('D:\wrf_simulation\\2meic\\wrfout_d03_2016-07-21_12_2meic')
#读取及初步处理数据
o3=getvar(ncfile,'o3',timeidx=126)
o3=o3*1000/22.4*48*273.15/(getvar(ncfile,'tk',timeidx=126))*(getvar(ncfile,'pressure',timeidx=126))/1013.25
lat=getvar(ncfile,'lat')
lon=getvar(ncfile,'lon')
height=getvar(ncfile,'height')
hgt=getvar(ncfile,'HGT')
height2earth=height-hgt
ua=getvar(ncfile,'ua',timeidx=126)
va=getvar(ncfile,'va',timeidx=126)
wa=getvar(ncfile,'wa',timeidx=126)
wa=wa*10
p=getvar(ncfile,'pressure',timeidx=126)
#插值数据
o3_vert=vertcross(o3,p,wrfin=ncfile,start_point=startpoint,end_point=endpoint,latlon=True)
ua_vert=vertcross(ua,p,wrfin=ncfile,start_point=startpoint,end_point=endpoint,latlon=True)
va_vert=vertcross(va,p,wrfin=ncfile,start_point=startpoint,end_point=endpoint,latlon=True)
wa_vert=vertcross(wa,p,wrfin=ncfile,start_point=startpoint,end_point=endpoint,latlon=True)
#计算风速
ws_vert=np.sqrt(ua_vert**2+va_vert**2)
#计算风向夹角
wdir_vert=np.arctan2(va_vert,ua_vert)*180/np.pi
line_angel=np.arctan2(end_lat-start_lat,end_lon-start_lon)*180/np.pi
vl_angel=wdir_vert-line_angel
#vl_angel = np.cos(vl_angel / 180 * np.pi)
vl_angel=np.cos(np.deg2rad(vl_angel))
ws_vert=ws_vert*vl_angel

lonlist,latlist,hlist=[],[],[]
plist=to_np(va_vert.coords['vertical'])
#生成插值数据的经纬度列表
foriinrange(len(ua_vert.coords['xy_loc'])):
s=str(ua_vert.coords['xy_loc'][i].values)
lonlist.append(float(s[s.find('lon=')+4:s.find('lon=')+12]))
latlist.append(float(s[s.find('lat=')+4:s.find('lat=')+12]))
#插值数据对应气压的离地高度
foriinrange(big_p,small_p-interval_p,-interval_p):
hlist.append(float(np.max(interplevel(height2earth,p,i)).values))
hlist=np.array([int(i)foriinhlist])#简单地取整
str_lonlist,float_lonlist=[],[]
a=np.mgrid[0:len(lonlist)-1:complex(str(int(round((end_lon+interval_lon-start_lon)/interval_lon)))+'j')]#生成位置列表
a=np.around(a,decimals=0)#四舍五入取整
#处理插入x轴的刻度标签
foriinrange(int((end_lon+interval_lon-start_lon)/interval_lon)):
float_lonlist.append(lonlist[int(a[i])])
lo,la=round(lonlist[int(a[i])],2),round(latlist[int(a[i])],2)
str_lonlist.append(str(lo)+'°E'+"\n"+str(la)+'°N')

fig=plt.figure(figsize=(5,5),dpi=150)
axe=plt.subplot(1,1,1)
axe.set_title('垂直剖面图，臭氧浓度与风向',fontsize=12,y=1.05)
axe.set_xlim(start_lon,end_lon)
axe.set_ylim(small_p,big_p)
axe.invert_yaxis()# 翻转纵坐标
axe.grid(color='gray',linestyle=':',linewidth=0.7,axis='y')
plt.xticks(float_lonlist,str_lonlist,fontsize=8,color='k')
axe.set_xlabel('经纬度坐标',fontproperties=Simsun,fontsize=12)
axe.get_xaxis().set_visible('True')
plt.yticks(fontsize=8,color='k')
axe.get_yaxis().set_visible('True')
axe.set_yticks(np.arange(small_p,big_p,interval_p))
axe.set_ylabel('气压$\mathrm{(hPa))}$',fontproperties=Simsun,fontsize=12)
axe.tick_params(length=2)
labels=axe.get_xticklabels()+axe.get_yticklabels()
[label.set_fontproperties(FontProperties(fname="./font/Times.ttf",size=8))forlabelinlabels]

axe.grid(color='gray',linestyle=":",linewidth=0.7,axis='y')

o3_level=np.arange(40,410,20)
contourf=axe.contourf(lonlist,plist,o3_vert,levels=o3_level,cmap=cmaps.NCV_jaisnd,extend='neither')
axe.quiver(lonlist,plist,ws_vert,wa_vert,pivot='mid',
width=0.001,scale=120,color='k',headwidth=4,
alpha=1)

fig.subplots_adjust(right=0.78)
rect=[0.87,0.07,0.01,0.57]# 分别代表，水平位置，垂直位置，水平宽度，垂直宽度
cbar_ax=fig.add_axes(rect)
cb=fig.colorbar(contourf,drawedges=True,ticks=o3_level,cax=cbar_ax,orientation='vertical',spacing='uniform')
cb.set_label('$\mathrm{O_3}$浓度',fontsize=12)
cb.ax.tick_params(length=0)
labels=cb.ax.get_xticklabels()+cb.ax.get_yticklabels()
[label.set_fontproperties(FontProperties(fname="./font/Times.ttf",size=10))forlabelinlabels]

axe1=fig.add_axes([0.8,0.69,0.2,0.2],projection=ccrs.PlateCarree())
axe1.add_feature(cfeat.COASTLINE.with_scale("10m"),linewidth=0.7,color='k')
axe1.add_feature(cfeat.RIVERS.with_scale("10m"))
axe1.add_feature(cfeat.OCEAN.with_scale("10m"))
axe1.add_feature(cfeat.LAKES.with_scale("10m"))
axe1.set_extent([120,123,29.5,32.5],crs=ccrs.PlateCarree())

gl=axe1.gridlines(crs=ccrs.PlateCarree(),draw_labels=False,linewidth=0.5,color='gray',linestyle=':')
gl.xlocator=mticker.FixedLocator(np.arange(120,123.1,1))
gl.ylocator=mticker.FixedLocator(np.arange(29.5,32.6,1))
axe1.set_xticks(np.arange(120,123.1,1),crs=ccrs.PlateCarree())
axe1.set_yticks(np.arange(29.5,32.6,1),crs=ccrs.PlateCarree())
plt.tick_params(top=True,bottom=True,left=True,right=True)#tick刻度显示
plt.tick_params(labeltop=True,labelleft=True,labelright=True,labelbottom=True)#tick标签显示
axe1.xaxis.set_major_formatter(LongitudeFormatter())
axe1.yaxis.set_major_formatter(LatitudeFormatter())
labels=axe1.get_xticklabels()+axe1.get_yticklabels()
[label.set_fontproperties(FontProperties(fname="./font/Times.ttf",size=6))forlabelinlabels]

axe1.plot([lonlist[0],lonlist[-1]],[latlist[0],latlist[-1]],color='red',linewidth=1.5,linestyle='-')
axe1.plot(lonlist[0],latlist[0],marker='o',color='red',markersize=2.5)
axe1.plot(lonlist[-1],latlist[-1],marker='o',color='red',markersize=2.5)
axe1.text(lonlist[0],latlist[0],'A',fontproperties=FontProperties(fname="./font/Times.ttf",size=8))
axe1.text(lonlist[-1],latlist[-1],'A'+"'",fontproperties=FontProperties(fname="./font/Times.ttf",size=8))

axe_sub=axe.twinx()
axe_sub.set_ylim(hlist[0],hlist[-1])
axe_sub.set_yticks(np.mgrid[hlist[0]:hlist[-1]:complex(str(int(hlist.shape[0]))+'j')])
round_hlist=np.around(hlist,-1)
axe_sub.set_yticklabels(round_hlist)
axe_sub.set_ylabel('离地高度$\mathrm{(m))}$',fontsize=12)
axe_sub.tick_params(labelcolor='black',length=3)
labels=axe_sub.get_xticklabels()+axe_sub.get_yticklabels()
[label.set_fontproperties(FontProperties(fname="./font/Times.ttf",size=8))forlabelinlabels]

plt.show()

原文链接：https://zhuanlan.zhihu.com/p/380626604