Matplotlib学习笔记：绘制矢量箭头图

meteorology/visual, programming/python
目录

气象上常见的矢量图有三种表现形式：箭头图、流线图、风羽图。 本文主要介绍矢量箭头图的绘制。

与填充图略微不同的是，矢量图需要准备两个分量的数据。 其余部分与填充图相同，都需要绘制底图。

准备数据

以读取 GRIB2 数据为例说明如何准备数据，与绘制填充图博文一样，使用 nuwe_pyeccodes 库读取 GRIB2 数据。

使用 GRAPES GFS 的 GRIB2 数据，序号 131 和 171 的场代表 850HPa 的 u、v 分量。

import nuwe_pyeccodes

# get message
grib_file = nuwe_pyeccodes.GribFileHandler()
grib_file.openFile("F:/data/gfs/gmf.gra.2019042500000.grb2")

# gh 850HPa
for i in range(0, 131):
    u_message = grib_file.next()

for i in range(131, 171):
    v_message = grib_file.next()

创建一个类 GridData 方便将 GRIB 消息转成 numpy 矩阵。

import numpy as np

class GridData(object):
    def __init__(self, grib_message):
        left_lon = grib_message.getDouble('longitudeOfFirstGridPointInDegrees')
        right_lon = grib_message.getDouble('longitudeOfLastGridPointInDegrees')
        lon_step = grib_message.getDouble('iDirectionIncrementInDegrees')
        nx = grib_message.getLong('Ni')
        lon_array = np.arange(left_lon, right_lon + lon_step / 2, lon_step)

        top_lat = grib_message.getDouble('latitudeOfFirstGridPointInDegrees')
        bottom_lat = grib_message.getDouble('latitudeOfLastGridPointInDegrees')
        lat_step = grib_message.getDouble('jDirectionIncrementInDegrees')
        ny = grib_message.getLong('Nj')
        lat_array = np.arange(top_lat, bottom_lat - lat_step / 2, -lat_step)

        self.lons, self.lats = np.meshgrid(lon_array, lat_array)
        values = grib_message.getDoubleArray('values')
        self.grid_values = values.reshape(ny, nx)

使用GribData生成Matplotlib需要的数据。

u_data = GridData(u_message)
v_data = GridData(v_message)

绘制箭头图

Matplotlib提供 quiver 绘制箭头图。 使用 cartopy 为数据添加投影，并绘制地图。

注意直接使用原始数据绘图会导致箭头过多，需要对数据进行采样。

import matplotlib.pyplot as plt
import cartopy.crs as ccrs

vector_crs = ccrs.PlateCarree()

f = plt.figure(figsize=(10, 5))

ax = f.add_subplot(
    111,
    projection=ccrs.PlateCarree(central_longitude=150)
)

ax.quiver(
    u_data.lons[::20, ::20], u_data.lats[::20, ::20],
    u_data.grid_values[::20, ::20], v_data.grid_values[::20, ::20],
    transform=vector_crs,
    color='blue',
)

ax.coastlines()
plt.show()

结果如下图所示。