绘图精进 | ERA5 弯曲箭头风场图绘制

用户11172986

发布于 2025-01-13 12:29:32

11600

代码可运行

文章被收录于专栏：气python风雨气python风雨

运行总次数：0

代码可运行

python | ERA5 弯曲箭头风场图绘制

前言

前两日看到深雨露公众号发了个弯曲箭头的文章，接着云台书使也发了

那我也凑凑热闹吧

弯曲箭头从前是NCL语言独有的绘制方法，在2019后它不再维护，NOAA转身投入了PYTHON的怀抱

于是弯曲箭头不再是NCL的专属，我们PYTHON也有自己的curved vector了

小编将其称为寝取的curved vector

废话不多说，show u the code

摘要

项目使用pyngl 主要进行ERA5风场数据弯曲箭头的可视化

分为两部分，首先将官方示例分模块讲解

其次使用pyngl对常见的era5数据进行风场的弯曲箭头可视化

预览效果如下

Image Name

Image Name

个人信息

公众号：气python风雨

Image Name

关注我获取更多学习资料，第一时间收到我的Python学习资料，也可获取我的联系方式沟通合作

环境设置

安装依赖

由于和鲸的气象镜像3.9已配置好了pyngl库与pynio库,那么我们跳过这一步,

具体怎么配置网络上也有诸多帖子，不再赘述

!pip install xxx

导入必要的库

import Ngl
import Nio
import numpy as np
import os
from IPython.display import Image  # 用于在Notebook中显示图像

数据读取

加载纬向风和经向风数据

# 打开NetCDF文件
dirc = Ngl.pynglpath("data")
file = Nio.open_file(os.path.join(dirc, "cdf", "pop.nc"))

# 获取变量
urot = file.variables["urot"]
vrot = file.variables["vrot"]
lat2d = file.variables["lat2d"]
lon2d = file.variables["lon2d"]
print(urot)

Variable: urot
Type: float
Total Size: 491520 bytes
            122880 values
Number of Dimensions: 2
Dimensions and sizes: [nlat | 384] x [nlon | 320]
Coordinates: 
            nlat: not a coordinate variable
            nlon: not a coordinate variable
Number of Attributes: 8
         missing_value : 9.96921e+36
         _FillValue : 9.96921e+36
         cell_methods : time: mean
         coordinates : lat2d lon2d
         units : centimeter/s
         long_name : Zonal Velocity
         time : 365031
         z_t : 500.622

好的，是384 x 320的二维数据

数据处理

提取目标区域数据

u = Ngl.add_cyclic(urot[290:])
v = Ngl.add_cyclic(vrot[290:])
lon = Ngl.add_cyclic(lon2d[290:])
lat = Ngl.add_cyclic(lat2d[290:])

数据可视化

绘制弯曲箭头图

Ngl.vector_map 创建并绘制地图上的矢量图。

vmap = Ngl.vector_map(wks, u, v, res=None)

参数

wks 调用 Ngl.open_wks 返回的工作站标识符。

u, v 矢量的 U 和 V 分量。u 和 v 应为二维的 NumPy 数组或 NumPy 掩码数组（维度为 ny x nx）。

res=None （可选）一个 Resources 类的实例，包含 PyNGL 资源作为属性。

返回值

vmap 表示创建的矢量图的 PlotId。

描述

此函数在给定的工作站上创建并绘制地图上的矢量图，并推进帧。可以通过 res 变量设置绘图选项。

为了在地图上叠加矢量图，必须告诉矢量图它在地图上的位置（以纬度/经度度数为单位）。可以通过 vfXArray/vfYArray 资源或 vfXCStartV/vfXCEndV/vfYCStartV/vfYCEndV 资源来实现。

从版本 1.3.0 开始，如果 u 和/或 v 是掩码数组，则任何等于相应填充值的值将不会被绘制。如果 u 和/或 v 不是掩码数组且包含缺失值，则应将资源 vfMissingUValueV 和/或 vfMissingVValueV 设置为这些值。

请注意，PyNGL 会根据其他资源的设置自动为您设置一些资源。有关更多信息，请参阅默认设置列表。

关键点总结：

功能：绘制地图上的矢量图。
输入：
- wks：工作站标识符。
- u 和 v：矢量的 U 和 V 分量，必须是二维数组或掩码数组。
- res：可选参数，用于设置绘图资源。
输出：返回一个表示矢量图的 PlotId。
注意事项：
- 必须通过 vfXArray/vfYArray 或 vfXCStartV/vfXCEndV/vfYCStartV/vfYCEndV 指定矢量图的位置。
- 如果数据是掩码数组，缺失值将不会被绘制。
- 如果数据不是掩码数组但包含缺失值，需设置 vfMissingUValueV 和 vfMissingVValueV。

# 创建工作站（输出为PNG文件）
wks_type = "png"# 指定输出文件的类型为PNG格式
wks = Ngl.open_wks(wks_type, "curly_vectors")  # 创建一个工作站，输出文件名为 "curly_vectors.png"

# 设置资源
vcres = Ngl.Resources()  # 创建一个资源对象，用于设置绘图的属性

# 设置坐标数组
vcres.vfXArray = lon  # 设置经度坐标数组，用于指定矢量的X轴位置
vcres.vfYArray = lat  # 设置纬度坐标数组，用于指定矢量的Y轴位置

# 设置地图投影和填充颜色
vcres.mpProjection = "Orthographic"# 设置地图投影为正交投影（Orthographic）
vcres.mpFillOn = True# 启用地图填充
vcres.mpLandFillColor = "Tan1"# 设置陆地的填充颜色为 "Tan1"
vcres.mpOceanFillColor = "SkyBlue"# 设置海洋的填充颜色为 "SkyBlue"
vcres.mpInlandWaterFillColor = "SkyBlue"# 设置内陆水域的填充颜色为 "SkyBlue"
vcres.mpLimitMode = "LatLon"# 设置地图范围模式为经纬度模式
vcres.mpCenterLonF = -80.0# 设置地图中心的经度为 -80.0 度
vcres.mpCenterLatF = 55# 设置地图中心的纬度为 55 度
vcres.mpMinLatF = 60# 设置地图的最小纬度为 60 度
vcres.mpDataBaseVersion = "MediumRes"# 设置地图数据库版本为中等分辨率

# 设置向量资源
vcres.vcMinDistanceF = 0.013# 控制向量之间的最小距离，用于调整向量密度
vcres.vcGlyphStyle = "CurlyVector"# 设置向量样式为弯曲箭头（CurlyVector）
vcres.vcRefMagnitudeF = 30.0# 设置参考向量的大小为 30.0
vcres.vcRefLengthF = 0.08# 设置参考向量的长度为 0.08
vcres.tiMainString = "Curly Vectors"# 设置图的主标题为 "Curly Vectors"

# 绘制弯曲箭头
plot = Ngl.vector_map(wks, u, v, vcres)  # 在工作站上绘制矢量图，u 和 v 分别是经向和纬向风分量

# 结束绘图
Ngl.end()  # 结束绘图会话，释放资源并关闭工作站

# 在Notebook中显示图像
Image(filename="curly_vectors.png")  # 在Jupyter Notebook中显示生成的PNG图像

完整代码

import Ngl
import Nio
import numpy as np
import os
from IPython.display import Image  # 用于在Notebook中显示图像

# 打开NetCDF文件
dirc = Ngl.pynglpath("data")
file = Nio.open_file(os.path.join(dirc, "cdf", "pop.nc"))

# 获取变量
urot = file.variables["urot"]
vrot = file.variables["vrot"]
lat2d = file.variables["lat2d"]
lon2d = file.variables["lon2d"]
print(urot)
# 添加循环点
u = Ngl.add_cyclic(urot[290:])
v = Ngl.add_cyclic(vrot[290:])
lon = Ngl.add_cyclic(lon2d[290:])
lat = Ngl.add_cyclic(lat2d[290:])

# 创建工作站（输出为PNG文件）
wks_type = "png"
wks = Ngl.open_wks(wks_type, "curly_vectors")

# 设置资源
vcres = Ngl.Resources()

# 设置坐标数组
vcres.vfXArray = lon
vcres.vfYArray = lat

# 设置地图投影和填充颜色
vcres.mpProjection = "Orthographic"
vcres.mpFillOn = True
vcres.mpLandFillColor = "Tan1"
vcres.mpOceanFillColor = "SkyBlue"
vcres.mpInlandWaterFillColor = "SkyBlue"
vcres.mpLimitMode = "LatLon"
vcres.mpCenterLonF = -80.0
vcres.mpCenterLatF = 55
vcres.mpMinLatF = 60
vcres.mpDataBaseVersion = "MediumRes"

# 设置向量资源
vcres.vcMinDistanceF = 0.013# 控制向量密度
vcres.vcGlyphStyle = "CurlyVector"# 使用弯曲箭头
vcres.vcRefMagnitudeF = 30.0# 参考向量大小
vcres.vcRefLengthF = 0.08# 参考向量长度
vcres.tiMainString = "Curly Vectors"# 标题

# 绘制弯曲箭头
plot = Ngl.vector_map(wks, u, v, vcres)
# 结束绘图

Ngl.end()

# # 在Notebook中显示图像
Image(filename="curly_vectors.png")

Variable: urot
Type: float
Total Size: 491520 bytes
            122880 values
Number of Dimensions: 2
Dimensions and sizes: [nlat | 384] x [nlon | 320]
Coordinates: 
            nlat: not a coordinate variable
            nlon: not a coordinate variable
Number of Attributes: 8
         missing_value : 9.96921e+36
         _FillValue : 9.96921e+36
         cell_methods : time: mean
         coordinates : lat2d lon2d
         units : centimeter/s
         long_name : Zonal Velocity
         time : 365031
         z_t : 500.622

ERA5示例

在xarray发展起来的当下，我们可以选择熟悉的xarray库读取数据，而不是使用pynio

反正填入绘图函数的参数numpy也可

xarray读取

import xarray as xr
import numpy as np
import Ngl
import Nio
import numpy as np
import os
from IPython.display import Image
ds = xr.open_dataset("/home/mw/input/era58091/ERA5-2023-08_pl.nc")
u= ds.u[0,30]   ## 850hpa,也可用sel获取更优雅些
v= ds.v[0,30]   ## 850hpa

lat = ds["latitude"][:]  # 纬度
lon = ds["longitude"][:]  # 经度
lon_grid, lat_grid = np.meshgrid(lon,lat)