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

目录

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

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

准备数据

以读取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()

结果如下图所示。