用Python画ERA5气象要素分布图

·
日常 no tag May 10, 2021

往期Python基础教程:Python教程(因为换了服务器图挂了……)

本文以ERA5月资料为例,讲解使用Python简单处理和绘制气象要素空间分布图

数据获取

从Climate Data Store下载ERA5月平均资料,所选配置为:
QQ截图20210510175450.png

  • 变量:位势高度、温度
  • 高度层:500、700、850hPa
  • 时间:2020年1-12月
  • 区域:全球
  • 格式:NC格式

QQ截图20210510180512.png

数据请求完成后下载,并改名为“2020_data.nc”,为了方便,在于数据同一文件目录下新建Python脚本,命名为"draw.py"。

QQ截图20210510180728.png

数据读取

使用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个月的数据平均即可。首先说明如何读取数据,对于温度或位势高度这样一个气象要素来说,在某一层上是一个三维矩阵。即某一高度上,温度会随着经度的不同,纬度的不同,时间的不同都在变化。

QQ截图20210510184741.jpg

如果我们想要获取指定时间范围,指定地理范围的数据,就需要计算各个纬度下所需值的序号。例如我只需要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个维度进行平均。什么意思呢,这里举一个例子:
QQ截图20210510200507.png
假设有这样一个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坐标和格点的值,运行得到这样的结果:
微信截图_20210510203952.png

地图底图

结果还不错,那么我们就加入地图底图吧。地图的绘制使用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()

微信截图_20210510204517.png
其中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()

微信截图_20210510210055.png

  • 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()

微信截图_20210510210849.png
只要第四个参数与等值线的plt.contour()相同,那么填色区域就在等值线之间。现在这个颜色条太大了,颜色也不好看,可以美化一下:

plt.contourf(lon, lat, avg_z, np.arange(53000, 58000, 400), cmap="rainbow")
plt.colorbar(shrink=0.8)

微信截图_20210510211223.png
这里为了使颜色能充分体现,将填色范围改成了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()
  • FF14集体模式下视野大小计算
  • 皆忆——一款辅助记忆小程序
取消回复

说点什么?
Title
数据获取
数据读取
数据处理
绘图
等值线
地图底图
等值线填色

© 2022 夜航船 · TOYOHAY Clouds. Using Typecho & Moricolor | 粤ICP备18131337号.