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

Question

I am trying to plot a CMC grib2 pressure forecast file using matplotlib to plot the pressure contours. The description of the grib2 grid can be found here: 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/ and starts with CMC_reg_PRMSL_MSL_0_ps10km followed by the date. It is a grib file containing pressure at mean sea level.

My problem is that I end up having some straight line contours that follow the lines of latitude on top of the actual pressure contours. I thought it might be because I am plotting in PlateCarree as opposed to Geodetic but the contour plot will not allow using Geodetic. The result of my plot is:

Code is as follows:

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)

Any suggestions would be appreciated.

UPDATE

For anyone without PyNIO the following can be used to reproduce using the dump files in the comments section.

Just remove all the references to NIO and replace the lats, lons, obs assignment with the following.

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

Answer 1

The problem

The problem is that the grid winds around the earth. Hence there will be points on the grid at -180° whose nearst neighbor sits at +180°, i.e. the grid wraps around the antimeridian. The following plots the grid index along both directions. One can see that the first grid row (black) appears on both sides of the plot.

Hence a contour line following the pacific westwards needs to then cross straight through the plot to continue towards japan on the other side of the plot. This will lead to the undesired lines

A solution

A solution is to mask the outer points of the PlateCarree out. Those occur in the middle of the grid. Cutting the grid at coordinates of longitude larger than 179° or smaller than -179°, as well as leaving the north pole out would look like

where the blue denotes the cut out points.

Applying this to the contour plot gives:

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

If you are sum up your longitude by +180 to avoid negative coordinates, your code should be running. A coordinate transformation should be legit from my point of view.