用Python画ERA5气象要素分布图
2022-07-14更新：文章作者制作了一个基于 matplotlib 与 cartopy 的 ERA5 数据快速画图封装库: LiMeteo. LiMeteo文档
往期Python基础教程：Python教程(因为换了服务器图挂了……)
本文以ERA5月资料为例，讲解使用Python简单处理和绘制气象要素空间分布图
数据获取
从Climate Data Store下载ERA5月平均资料，所选配置为：
- 变量：位势高度、温度
- 高度层：500、700、850hPa
- 时间：2020年1-12月
- 区域：全球
- 格式：NC格式
数据请求完成后下载，并改名为“2020_data.nc”，为了方便，在于数据同一文件目录下新建Python脚本，命名为"draw.py"。
数据读取
使用IDE打开脚本，我使用的是(PyCharm)[https://www.jetbrains.com/pycharm/download]社区版，不建议使用盗版的PyCharm专业版，与免费的社区版相比，里面的绝大部分功能都不会用到，科学视图也不一定是有用的工具。
首先导入需要的库，读取nc文件本文使用netCDF4库，numpy库用于数据处理，matplotlib库用于绘图，basemap库用于画地图。
from netCDF4 import Dataset
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.basemap import Basemapfrom XXX import YYY 是指从某个库中导入其中YYY部分，as 是可以将导入的库的名称简写。
读取NC文件，使用我们导入的Dataset函数：
data = Dataset("2020_data.nc")双引号中是nc文件与py脚本的相对路径，因为两者在同一目录下，因此直接写nc文件名即可。这样就将数据储存在data变量中。
因为我们所选的数据，有不同的时间、空间位置和变量，因此是一个多维度的数据。想查看数据的维度属性，可以使用：
print(data.variables.keys())其中variables是读取数据不同维度变量的函数，返回的是字典格式的数据，这个函数是属于Dataset的，因为我们把Dataset赋值给了data，因此使用它时用.连接在data后面，keys()则是返回字典格式的keys函数。运行后输出：
dict_keys(['longitude', 'latitude', 'level', 'time', 'z', 't'])这里输出的就是下载的数据的各个变量，有经度、纬度、高度、时间、位势高度和温度，如果我们想查看其中某一变量的属性，例如高度，可以使用:
print(data.variables["z"])这里的中括号是因为variables函数返回的是字典格式，因此可以使用中括号检索，运行后返回：
<class 'netCDF4._netCDF4.Variable'>
int16 z(time, level, latitude, longitude)
scale_factor: 0.746896363282621
add_offset: 33853.14608306836
_FillValue: -32767
missing_value: -32767
units: m**2 s**-2
long_name: Geopotential
standard_name: geopotential
unlimited dimensions:
current shape = (12, 3, 721, 1440)
filling on这样就能看到各个变量的属性。若我们查看时间变量的属性和所有的值：
print(data.variables["time"]) # 输出时间维度属性
print(data.variables["time"][:]) # 输出时间的所有值这里的[:]为数组的切片引用，完整格式为[a:b]，即截取数组的第a到b-1个，若ab为空，即选取所有数组元素。运行输出：
<class 'netCDF4._netCDF4.Variable'>
int32 time(time)
units: hours since 1900-01-01 00:00:00.0
long_name: time
calendar: gregorian
unlimited dimensions:
current shape = (12,)
filling on, default _FillValue of -2147483647 used
[1051896 1052640 1053336 1054080 1054800 1055544 1056264 1057008 1057752
1058472 1059216 1059936]注意到时间的值都是一百多万的数字，这是因为ERA5资料的时间值是按照从1900年1月1日0时算起的以小时为单位的数值，即距离1900年1月1日0时过去了多少小时。关于时间的正确引用会在下文说明。
数据处理
这里以求年平均为例，因为下载的是月数据，因此把12个月的数据平均即可。首先说明如何读取数据，对于温度或位势高度这样一个气象要素来说，在某一层上是一个三维矩阵。即某一高度上，温度会随着经度的不同，纬度的不同，时间的不同都在变化。
如果我们想要获取指定时间范围，指定地理范围的数据，就需要计算各个纬度下所需值的序号。例如我只需要100-150E，10-50N这样一个范围的资料，我们首先用上面的方法，输出经度和纬度的所有值：
[0. 0.25 0.5 ... 359.25 359.5 359.75]
[90. 89.75 89.5 ... -89.5 -89.75 -90.]因为数组太大，中间就省略了，可以看出资料的分辨率的0.25度的。现在需要找出100N、150E、10、50N在里面的序号，以经度为例，100E对应的应该是第400个，150E对应的是第600个，同理纬度为第160个和第320个（纬度是从到小的，所以是从50N到10N）。若要验证是否正确，可以直接读取输出：
print(data.variables["longitude"][400], data.variables["longitude"][600])
print(data.variables["latitude"][160], data.variables["latitude"][320])运行输出：
100.0 150.0
50.0 10.0现在，提取这个区域的位势高度数据，在上面查看位势高度属性时，第二行有z(time, level, latitude, longitude)，这是说明高度场是一个四维矩阵，按顺序分别是时间、高度、纬度和经度，我们读取的时候，也要按照这个顺序来。读取数据与上面的读取所有值一样，只不过不是使用[:]，因为这是读取全部数据，我们需要的是部分区域的：
z = data.variables["z"][:, 0, 160:320, 400:600]
print(np.shape(z))这里中括号内，被逗号分成四个部分，其实对应的就是z(time, level, latitude, longitude)这几个纬度，第一个时间纬度我们输入了:，是因为要提取所有时间节点来做平均。第二个高度纬度我们输入了0，其实是代表了500hPa高度，选取不同高度层可以像经纬度一样通过输入对应层的序号来选取；第三和第四个就是纬度和经度，输入我们上面计算好的序号。
在时间维度上，如果我们只需要几个月的时间，也可以通过计算序号来提取。因为下载的资料是月数据，所以每一个时间节点间隔都是一个月，已知第一个时间节点为2020年1月，即1月的序号为0，因此第i个月的序号就是i-1。例如上面只想要6-8月的资料，把:替换为5:8即可。（注意：[a:b]提取的是第a至b-1个元素，第b个元素是取不到的，所以第b个元素也提取的话，需要[a:b+1]，上面的经纬度范围也是取不到150E和10N的，若需要，则更换为160:321, 400:601）np.shape()函数可以返回我们提取后的矩阵长度，运行输出:
(12, 160, 200)这个结果说明我们提取之后的数据，变成了三个维度，因为我们在高度层只选了一层，所以高度层的维度就消失了，或者说提取后的所有数据都是在500hPa上，所以不需再要高度层这个维度了。第一个12是代表有12个时间节点，剩下两个就是经纬度的格点数。
以资料格点的某一点来看，这个点上的数据其实就是这个地方12个月每个月的值，我们把12个值加以平均，就得到这个格点的年平均值。对每一个格点都这样处理的话，就能得到全部格点的均值。但是这样一个个处理太麻烦，我们可以直接使用更加方便的函数来处理：
avg_z = np.mean(z, axis=0)
print(np.shape(avg_z))这里的np.mean()函数，就是可以很方便地对数据矩阵进行平均处理的函数。第一个参数是z，也即是待处理的数据矩阵，参数axis是指平均的维度，axis=0是指以第0个维度进行平均。什么意思呢，这里举一个例子:
假设有这样一个2×5的矩阵，若不使用axis参数，则返回矩阵中所有元素的均值；axis=0时，以列为轴，即所有行的第i个元素的均值，在这个例子里就是[1与6的均值, 2与7的均值, ...]；axis=1是，以行为轴，计算每一行的均值，即[1、2、3、4、5的均值, 6、7、8、9、10的均值]。而我们提取后的数据，时间纬度是在第一个的，即是序号为0，所以我们要在时间上平均，就输入axis=0。上面的代码运行后输出：
(160, 200)这说明，在时间纬度平均后，时间纬度也消失了，我们得到的就是年平均的一个二维资料，两个维度分别是纬度和经度，这样，我们就可以以二维图片的形式将这个结果画出来了。
顺便一提，可以思考一下上面平均处理时，若axis为1或2，得到的结果是什么呢？
绘图
等值线
我们上面得到的，是一个160×200的矩阵，要将他们画到地图上，首先需要对他们定位，即告诉程序每一个格点的位置。每一个格点的位置，其实就是原来数据里面的那些经纬度，我们只要把他们按照截取的区域对应地读出来就行了。
lon = data.variables["longitude"][400:600]
lat = data.variables["latitude"][160:320]以相同的序号进行提取后，avg_z中的第i行第j个格点的经度就是lon的第j个，lat的第i个。
接着我们就可以尝试着把它们画出来了。
首先我们先尝试着画等值线图，使用的是plt.contour()函数，到目前为止，我们整个脚本是这样的：
from netCDF4 import Dataset
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.basemap import Basemap
data = Dataset("2020_data.nc") # 读取NC文件
z = data.variables["z"][:, 0, 160:320, 400:600] # 提取指定区域和高度的位势高度场
avg_z = np.mean(z, axis=0) # 求年均值
lon = data.variables["longitude"][400:600] # 格点经度坐标
lat = data.variables["latitude"][160:320] # 格点纬度坐标
plt.contour(lon, lat, avg_z) # 画等值线图
plt.show() # 查看画图结果
最后两行是画图的语句，plt.contour()需要3个必要的参数，分别是格点的x坐标，y坐标和格点的值，运行得到这样的结果：
地图底图
结果还不错，那么我们就加入地图底图吧。地图的绘制使用Basemap函数：
m = Basemap(llcrnrlat=10, llcrnrlon=100, urcrnrlat=50, urcrnrlon=150)
m.drawcoastlines()
plt.contour(lon, lat, avg_z)
plt.show()Basemap函数是创建了一个地图，地图的大小由里面的llcrnrlat, llcrnrlon, urcrnrlat, urcrnrlon，分别是地图左下角的纬度和经度以及右上角的纬度和经度，对应上面我们所需要的区域输入即可。第二行的m.drawcoastlines()是画海岸线。
可以注意到等值线的颜色都是不一样的，不同数值对应不同颜色，如果我们想让它都变成黑色，再加上数值表示，可以这样写：
m = Basemap(llcrnrlat=10, llcrnrlon=100, urcrnrlat=50, urcrnrlon=150)
m.drawcoastlines()
C = plt.contour(lon, lat, avg_z, colors='k')
plt.clabel(C)
plt.show()
其中colors参数指定了等值线的颜色，查看更多颜色值可以点我。这里使用了变量C来储存等值线，然后使用plt.clabel()函数画上等值线的数值，这个函数是专门给等值线画数值的，因此输入C。
之后我们再稍加美化一下：
plt.figure(figsize=(8, 6))
m = Basemap(llcrnrlat=10, llcrnrlon=100, urcrnrlat=50, urcrnrlon=150)
m.drawcoastlines(linewidth=0.7)
m.drawmeridians([100, 110, 120, 130, 140, 150], labels=[0, 0, 0, 1])
m.drawparallels([10, 20, 30, 40, 50], labels=[1, 0, 0, 0])
C = plt.contour(lon, lat, avg_z, np.arange(53000, 60000, 400), colors='k')
plt.clabel(C, fmt="%.0f")
plt.show()
plt.figure(figsize=(8, 6))设置了画板的大小为800×600像素，这个函数可用于创建和设置画板。- 在
m.drawcoastlines()中加入了linewidth=0.7参数，是设置海岸线线条宽度，使其较细。 m.drawmeridians()和m.drawparallels()是画经度和纬度网格，第一个参数是要画的值，第二个labels参数是在四个边显示对应的经纬度值，数组分别对应[上，右，下，左]。- 在
plt.coutour()中，第四个参数代表的是要画的等值线的值，是一个数组，这里用了np.arange()来生成了从53000到60000每隔400一个值的数组。 plt.clabel()中的fmt参数可以指定显示数字格式，"%a.bf"表示小数点前显示a位，小数点后显示b位。
等值线填色
填色用的函数是plt.contourf()，它的前四个参数与plt.contour都是一模一样的，填色后，可以使用plt.colorbar()来绘制色条。
plt.figure(figsize=(8, 6))
m = Basemap(llcrnrlat=10, llcrnrlon=100, urcrnrlat=50, urcrnrlon=150)
m.drawcoastlines(linewidth=0.7)
m.drawmeridians([100, 110, 120, 130, 140, 150], labels=[0, 0, 0, 1])
m.drawparallels([10, 20, 30, 40, 50], labels=[1, 0, 0, 0])
plt.contourf(lon, lat, avg_z, np.arange(53000, 60000, 400))
plt.colorbar()
C = plt.contour(lon, lat, avg_z, np.arange(53000, 60000, 400), colors='k')
plt.clabel(C, fmt="%.0f")
plt.show()
只要第四个参数与等值线的plt.contour()相同，那么填色区域就在等值线之间。现在这个颜色条太大了，颜色也不好看，可以美化一下：
plt.contourf(lon, lat, avg_z, np.arange(53000, 58000, 400), cmap="rainbow")
plt.colorbar(shrink=0.8)
这里为了使颜色能充分体现，将填色范围改成了53000~57600，cmap参数定义了填色使用的颜色，更多cmap预设点我，colorbar中的shrink为大小参数，使其变小一点。
至此，简易的绘图流程就结束啦，想要美化或绘制更复杂的图表和参数可参考Matplotlib官方文档（英文），地图底图可参考Basemap文档（英文），下面是完整的代码：
from netCDF4 import Dataset
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.basemap import Basemap
data = Dataset("2020_data.nc")
z = data.variables["z"][:, 0, 160:320, 400:600]
avg_z = np.mean(z, axis=0)
lon = data.variables["longitude"][400:600]
lat = data.variables["latitude"][160:320]
plt.figure(figsize=(8, 6))
m = Basemap(llcrnrlat=10, llcrnrlon=100, urcrnrlat=50, urcrnrlon=150)
m.drawcoastlines(linewidth=0.7)
m.drawmeridians([100, 110, 120, 130, 140, 150], labels=[0, 0, 0, 1])
m.drawparallels([10, 20, 30, 40, 50], labels=[1, 0, 0, 0])
plt.contourf(lon, lat, avg_z, np.arange(53000, 58000, 400), cmap="rainbow")
plt.colorbar(shrink=0.8)
C = plt.contour(lon, lat, avg_z, np.arange(53000, 58000, 400), colors='k')
plt.clabel(C, fmt="%.0f")
plt.show()