用Python画ERA5气象要素分布图
往期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 Basemap
from 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()