Argthtjtr

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')

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()