如何在 Python 中绘制极地立体 grib2 文件的等高线

Question

我正在尝试使用 matplotlib 绘制 CMC grib2 压力预测文件来绘制压力等值线。可以在此处找到 grib2 网格的描述：https://weather.gc.ca/grib/grib2_reg_10km_e.html. The grib2 file is found in this directory: http://dd.weather.gc.ca/model_gem_regional/10km/grib2/00/000/ 并以 CMC_reg_PRMSL_MSL_0_ps10km 开头，后跟日期。这是一个包含平均海平面压力的 grib 文件。

我的问题是我最终得到了一些直线等高线，这些等高线在实际压力等高线之上遵循纬度线。我认为这可能是因为我在 PlateCarree 中绘图而不是 Geodetic，但等高线图不允许使用 Geodetic。我的情节的结果是：

代码如下：

import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import datetime as dt
import cartopy
import cartopy.crs as ccrs
import Nio

gr = Nio.open_file('./data/CMC_reg_PRMSL_MSL_0_ps10km_2018111800_P000.grib2', 'r')
print(gr)
names = gr.variables.keys()
print("Variable Names:", names)
dims = gr.dimensions
print("Dimensions: ", dims)
attr = gr.attributes.keys()
print("Attributes: ", attr)

obs = gr.variables['PRMSL_P0_L101_GST0'][:]
lats = gr.variables["gridlat_0"][:]
lons = gr.variables["gridlon_0"][:]

fig = plt.figure(figsize=(15, 2))
intervals = range(95000, 105000, 400)
ax=plt.axes([0.,0.,1.,1.],projection=ccrs.PlateCarree())
obsobj = plt.contour(lons, lats, obs, intervals, cmap='jet',transform=ccrs.PlateCarree())
states_provinces = cartopy.feature.NaturalEarthFeature(
        category='cultural',
        name='admin_1_states_provinces_lines',
        scale='50m',
        facecolor='none')
ax.add_feature(cartopy.feature.BORDERS)
ax.coastlines(resolution='10m')  
ax.add_feature(states_provinces,edgecolor='gray')
obsobj.clabel()
colbar =plt.colorbar(obsobj)

如有任何建议，我们将不胜感激。

更新

对于没有 PyNIO 的任何人，以下内容可用于使用评论部分中的转储文件进行重现。

只需删除所有对 NIO 的引用，并将 lats、lons、obs 赋值替换为以下内容。

lats = np.load('lats.dump')
lons = np.load('lons.dump')
obs = np.load('obs.dump')

Answer 1

问题

问题是网格缠绕在地球周围。因此，网格上的点在 -180° 处，其最近的邻居位于 +180° 处，即网格环绕反子午线。下面绘制了沿两个方向的网格索引。可以看到第一个网格行（黑色）出现在图的两边。

因此，沿着太平洋向西的等高线需要直接穿过该地块，然后继续向地块另一侧的日本前进。这将导致不需要的行

一个解决方案

一种解决方案是将 PlateCarree 的外部点屏蔽掉。那些发生在网格的中间。在大于 179° 或小于 -179° 的经度坐标处切割网格，以及将北极留在外面看起来像

其中蓝色表示剪切点。

将其应用于等高线图得到：

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

lats = np.load('data/lats.dump')
lons = np.load('data/lons.dump')
obs = np.load('data/obs.dump')

intervals = range(95000, 105000, 400)

fig, ax = plt.subplots(figsize=(15,4), subplot_kw=dict(projection=ccrs.PlateCarree()))
fig.subplots_adjust(left=0.03, right=0.97, top=0.8, bottom=0.2)

mask = (lons > 179) | (lons < -179) | (lats > 89)
maskedobs = np.ma.array(obs, mask=mask)

pc = ax.contour(lons, lats, maskedobs, intervals, cmap='jet', transform=ccrs.PlateCarree())

ax.add_feature(cartopy.feature.BORDERS)
ax.coastlines(resolution='10m')  

colbar =plt.colorbar(pc)

plt.show()

Answer 2

如果您将经度加起来 +180 以避免负坐标，您的代码应该是运行。从我的角度来看，坐标转换应该是合法的。

如何在 Python 中绘制极地立体 grib2 文件的等高线

How to plot contours from a polar stereographic grib2 file in Python

python

matplotlib

cartopy

grib

问题

一个解决方案