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

Я пытаюсь построить файл прогноза давления 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, за которым следует дата. Это гриб-файл, содержащий давление на среднем уровне моря.

Моя проблема заключается в том, что у меня в конечном итоге появляются контуры прямых линий, которые следуют линиям широты поверх контуров фактического давления. Я подумал, что это может быть из-за того, что я рисую в 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')

python matplotlib cartopy grib
person Campbell    schedule 18.11.2018    source источник
comment
К сожалению, я так не думаю (хотя я открыт для предложений о том, как этого можно достичь). Переменные представляют собой массивы 824x935 numpy, поэтому их немного сложно напрямую вставить в код. Поскольку я не уверен, почему он ведет себя так, я не знаю, как я мог бы создать меньший пример, чтобы воспроизвести проблему.   -  person Campbell    schedule 18.11.2018
comment
Это работает. Вот ссылки на файлы (по какой-то причине я получал ошибки, включая их в мое обновление вопроса: drive.google.com/open?id=1A7mZtqFcPCk92WRJ11TxOnHrsniZH_s- drive.google .com/open?id=15RC-IsoojBLIG9VMknTpFP1rPqV2qDN3 drive.google.com/ open?id=1An3oXDzev-4hGmapvYykduG293fhC1PH   -  person Campbell    schedule 18.11.2018
comment
Проблема в сетке, которая полностью деформирована в декартовых координатах и ​​накручивается вокруг северного полюса. Я построил индекс сетки в обоих измерениях здесь, чтобы визуализировать проблему. В принципе, эта проблема похожа на эту. Однако один ответ, похоже, не решает проблему; другой ответ (написанный мной) не может быть применен здесь из-за формы сетки.   -  person ImportanceOfBeingErnest    schedule 19.11.2018
comment
Я все еще очень новичок в данных с нумерацией и сеткой ... как вы думаете, можно ли отказаться от верхних 15-20 ° широты? Я мог обойтись только данными до 80° с.ш. Я просто не уверен, как найти, где сделать разрез и как это повлияет на контуры.   -  person Campbell    schedule 19.11.2018

Ответы (2)


arrow_upward
3
arrow_downward

Эта проблема

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

введите здесь описание изображения

person ImportanceOfBeingErnest    schedule 19.11.2018

arrow_upward
0
arrow_downward

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

person dl.meteo    schedule 04.03.2019