Как построить контуры из полярного стереографического файла grib2 в Python

Question

Я пытаюсь построить файл прогноза давления CMC grib2, используя matplotlib для построения контуров давления. Описание сетки 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 в отличие от геодезического, но контурный график не позволяет использовать геодезические. Результатом моего сюжета является:

Код выглядит следующим образом:

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)

Любые предложения будут оценены.

UPDATE

Для тех, у кого нет PyNIO, для воспроизведения можно использовать следующие файлы дампа в разделе комментариев.

Просто удалите все ссылки на NIO и замените lats, lons, obs на следующее.

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

Answer 1

Если вы суммируете свою долготу на +180, чтобы избежать отрицательных координат, ваш код должен работать.Преобразование координат должно быть законным с моей точки зрения.

Answer 2

Проблема

Проблема в том, что сетка обвивается вокруг земли.Следовательно, на сетке будут точки при -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()