如何在 Python 中从极地立体 grib2 文件中绘制轮廓答案

【问题标题】：How to plot contours from a polar stereographic grib2 file in Python如何在 Python 中从极地立体 grib2 文件中绘制轮廓
【发布时间】：2019-04-21 01:33:34
【问题描述】：

我正在尝试使用 matplotlib 绘制 CMC grib2 压力预测文件来绘制压力等值线。 grib2 网格的描述可以在这里找到：https://weather.gc.ca/grib/grib2_reg_10km_e.html。 grib2 文件位于此目录中：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 的任何人，可以使用 cmets 部分中的转储文件来重现以下内容。

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

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

【问题讨论】：

不幸的是，我不这么认为（尽管我对如何实现这一点持开放态度）。变量是 824x935 numpy 数组，因此直接放入代码有点困难。由于我不确定它为什么会这样，所以我不知道如何创建一个更小的示例来重现问题。
那行得通。以下是文件的链接（由于某种原因，我在问题更新中遇到了错误，包括它们：drive.google.com/open?id=1A7mZtqFcPCk92WRJ11TxOnHrsniZH_s-drive.google.com/open?id=15RC-IsoojBLIG9VMknTpFP1rPqV2qDN3 drive.google.com/open?id=1An3oXDzev-4hGmapvYykduG293fhC1PH
问题是网格，它在笛卡尔坐标中完全变形并且绕着北极旋转。我在两个维度上绘制了网格索引here 以可视化问题。原则上，这个问题似乎类似于this one。但是，一个答案似乎并不能解决问题；由于网格形状，其他答案（由我编写）无法在此处应用。
我对 numpy 和网格数据还是很陌生...您认为有可能丢弃纬度的上 15-20° 吗？我可以仅使用高达 80°N 的数据。我只是不确定如何找到在哪里进行切割以及对轮廓的影响。

标签： python matplotlib cartopy grib

【解决方案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()

【讨论】：

【解决方案2】：

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

【讨论】：