如何用python画金字塔形状_新手上路 | 用Python画一张高空天气图

由于之前没有找到微信这边怎么插代码，所以都是在知乎发的，现在找到了就搬运过来，之后会微信和知乎专栏一起同步发。

知乎专栏：天气瓶milky - 今天又在写代码

摸了一个暑假的鱼终于在八月快结束的时候学了点东西出来，写个文章算是总结一下这个月学到的东西吧，本来是想发在vx公众号上的，但是代码块的插入实在是太烦了，那就试试知乎的文章功能吧。

急匆匆看完了隔壁某大型交友视频网站的教程就开始动手了，过程中参考了不少大佬们的文章，感谢云台书使大佬不厌其烦地指导(下次玩安娜激素都给你)。

用了Anaconda来管理和安装需要的包，写代码则用的JupyterLab，看起来比较舒服。

使用ERA5 hourly data on pressure levels from 1979 to present的再分析资料，精度为0.25° * 0.25°，绘制了2020年5月21日02时(世界时)的位势高度、温度和风场。

成品图如下(中心气压的标志、槽线分析没有画。。。)：

在写代码的时候碰到了以下一些纠结了较长时间的问题，会在下面贴代码的时候顺便讲一下我自己的头铁解决方法，希望各位大佬多多指教。

一些在准备工作中出现的问题：

1. 1. 清华源换原始源；

2. 2. Cartopy自带的地图精度文件下载时间太长，出不了图；

3. 3. JupyterLab的默认文件夹和默认浏览器设置。

一些在画图过程中出现的问题：

1. 1. 直接使用set_extent会令等高度线的数值标签被切在画布外，如何让标签呈现在画布内；

2. 2. 风羽图相关：

   1.     (1)精度过高，风羽过密，如何降低风羽的密度；

   2.     (2)风羽的单位更改。

清华源换原始源

刚开始装的清华源，结果直接在Anaconda Prompt里使用conda install 命令安装Cartopy失败了，师兄说是清华源拉胯导致的。不推荐去Cartopy官网下载安装，万一装错了挺麻烦的。

在Anaconda Prompt中输入如下命令即可换回原始源：

$ conda config --remove-key channels

下载Cartopy中直接调用的地图信息

第一次画图的时候会自动下载到目录下，建议还是在第一次绘图的时候等大概五到十分钟让他自己下载，实在是太慢了没动静再考虑下面的手动安装，具体可以参考以下链接：

https://blog.csdn.net/x5675602/article/details/86429040

注意：该自带的地图信息中国国界存在问题，一定要使用正确的国界文件嗷。

默认设置修改

可参考：

https://blog.csdn.net/housz77/article/details/106125831

https://blog.csdn.net/luckywizard/article/details/102655407

正式开始画图了

import所需要的包

import matplotlib.pyplot as plt
import xarray as xr    #读取nc文件
import numpy as np
import cartopy.crs as ccrs   #投影方式
import cartopy.feature as cfeat   #使用shp文件
from cartopy.io.shapereader import Reader  #读取自定的shp文件
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter  #将x，y轴换成经纬度

为了防止无法显示中文，务必在开始加上这一行：

plt.rcParams['font.sans-serif']=['SimHei']   #黑体，也可以替换为其他字体

利用xarray打开nc文件，注意文件的绝对路径和相对路径。第二行可以查看该文件的信息(xarray的信息栏看起来比netCDF4的顺眼多了，建议还是用xarray)

obj = xr.open_dataset('2122h.nc')
obj

文件信息要仔细看，尤其是单位、开始的经纬度，后面出了什么问题很大可能是和这两者相关(比如换算的时候，我就被位势米到底怎么换给坑了...)

因为只是试试水，所以我只将时次0(2020.5.21 UTC02)和层次1(500hPa)的要素拿出来放在z、u、v、t里，同时进行了一下单位的换算和单位重置。

z = obj['z'][0][1]/98
z.attrs['units'] = 'dagpm'
u = obj['u'][0][1]
v = obj['v'][0][1]
t = obj['t'][0][1]-273.15
t.attrs['units'] = 'deg C'

“z”一下，可以看到现在z的属性情况，数据和单位已经改好了~

接着是读取经纬度，并且利用numpy.meshgrid()生成网格点坐标矩阵。这里是为了解决等压线的标签在所选区域外的准备工作，可以根据实际的情况去掉，直接用set_extent()来指定画图范围。

lat = obj['latitude'][:]
lon = obj['longitude'][:]
lons, lats = np.meshgrid(lon,lat)

画国界、省界，以及为底图上色。颜色可以利用对照表找自己喜欢的颜色，甚至直接安排颜色代码也是可以的(设计狗狂喜？)。

proj = ccrs.PlateCarree()
fig = plt.figure(figsize=(9,6),dpi=150)
ax = fig.subplots(1, 1, subplot_kw={'projection':proj})
country_shp = Reader('country1.shp')
pro_shp = Reader('行政边界_省级.shp')
fea_country = cfeat.ShapelyFeature(country_shp.geometries(), proj,
                              edgecolor='darkolivegreen', facecolor='ivory')
ax.add_feature(fea_country, linewidth=0.2, alpha=0.2)
fea_pro = cfeat.ShapelyFeature(pro_shp.geometries(), proj,
                              edgecolor='darkolivegreen', facecolor='none')
ax.add_feature(fea_pro, linewidth=0.3, alpha=0.5)
ax.add_feature(cfeat.OCEAN, alpha=0.5)

将x轴和y轴设置为经纬度格式，指定标出的经纬度。

ax.set_xticks([70,85,100,115,130])
ax.set_yticks([10,20,30,40,50,60])
ax.xaxis.set_major_formatter(LongitudeFormatter())
ax.yaxis.set_major_formatter(LatitudeFormatter())

画等压线：

z_levels = np.arange(500, 580+10, 4)
ac = ax.contour(lons[120:320,280:560], lats[120:320,280:560], z[120:320,280:560],
                levels=z_levels, extent='both',
                colors='mediumblue', linewidths=0.8)
plt.clabel(ac, inline=True, fontsize=8, fmt='%.0f')

利用levels可以规定等值线的间隔。

lons和lats是先前生成的网格点矩阵，逗号前为纬度，逗号后的是经度，为什么是这几个看起来莫名其妙的数据呢。。。我希望画的范围是10°N~60°N，70°E~135°E，该nc文件的精度是0.25° * 0.25°，因此一个经度/纬度实际上占了4个点，从对nc文件的描述中可以看到，其经度从0度开始向东走，因此要找到70°E只要4×70=280就可以了；但是对纬度来说，它是从90°N开始的，90°N是第1个，要推到60°N只能(90-60)×4，所以是120...

(我刚开始没发现，是一个个范围试出来的，为此还定位到了不知道哪个大洲，开学了得买份世界地图供着才行)

至于为什么不用set_extent()，其实我试过了，效果很棒，但是作为一个天气图，它的气压标注的不及格的，以900hPa这个图为例，可以看到左上角的标签都被切到看不见了，师兄给的解决办法是只取需要范围的数据画，这样标签就回到画布上了

2020年5月21日02时(世界时)900hPa 等高线

这样，通过网格点矩阵的定位，我们只取到了10°N~60°N，70°E~135°E这个范围中的z，就可以只画这个范围中的等高线了。

等温线也是同理，不过小于零的等温线是虚线，我们可以加一个np.where()来判断应该画什么线。(有warning？能画出来就行啦)

t_levels = np.arange(-40, 30, 4)
at = ax.contour(lons[120:320,280:560], lats[120:320,280:560], t[120:320,280:560],
                levels=t_levels, extent='both',
                colors='red', linewidths=0.5, linestyle=np.where(t>=0,'-','--'))
plt.clabel(at, inline=True, fontsize=5, fmt='%.0f')

最后是风羽图，如果我们和上面一样，直接使用

lons[120:320,280:560], lats[120:320,280:560]

这样的形式，那会出现一个直接起飞的结果——数据精度太高，画面变黑了...

所以要对每一个要素都取相同的间隔，再画图...(我偷懒没取相同间隔，又卡了很久)

另外，可以使用barb_increments来修改风矢杆长短杆线和三角分别代表的风速大小。^[5]

ax.barbs(lons[120:320:12,280:560:12], lats[120:320:12,280:560:12],
         u[120:320:12,280:560:12],v[120:320:12,280:560:12],
         barb_increments={'half':2,'full':4,'flag':20}, zorder=5, length=5, linewidth=0.2)

再加上标题，保存一哈，一张高空天气图就搞定了！

ax.set_title('2020年5月21日02时(世界时) 500hPa 温度(Temp) 高度(dagpm)',fontsize=12)
fig.show
fig.savefig('20052102_500hPa_天气图')

Python确实比GrADs好玩很多hhh，上手也很快，GrADs可定义的颜色确实有点少。在学了在学了，呜呜。

感谢阅读！第一次写这种教程，画个天气图其实也不是很难的东西，不过一步一步自己解决还挺有意思的，希望自己以后能画更多好看的图出来。

敬请各位大佬们批评指正！_(:3