Unable to assign y and x coordinates to xarray.DataArray

Question

After working through the MetPy Cross Section Example, I unsuccessfully tried to generalize that example to NCEP NAM-12km GRIB2 files. By comparing the DataArray for my file to the example file (a netCDF file), I find xarray.open_dataset() is not generating x- and y-coordinates. Nor is it finding sufficient information in the input GRIB file to generate a metpy_crs coordinate. I tried converting the GRIB file to netCDF but this did not resolve the problem.

Even treating the input file as non-CF complaint, I find I am unable to assign x- and y-coordinates. By inspecting the netCDF file used in the Cross-Section example, it is clear that the x and y arrays represent a spatial grid with GeoX and GeoY coordinate axis types. What is not clear to me, however, is how to create a analogous arrays for my file.

The latitude and longitude arrays are not amenable to attempts to collapse them manually (as in selecting the mode of latitudes along a row, which was one my attempts to create an appropriate coordinate vector.)

Some assistance with solving this issue will be greatly appreciated. Thank you!

Following is the code I am running.

import sys
#  pygrib is installed to the "herbie" environment, so we add the path
#  to this so we can load those packages.  Pygrib is used to get the 
#  complete set of keys for a message so we can try to see what is
#  available.

sys.path.append('/miniconda3/envs/herbie/lib/python3.8/site-packages')
import pygrib

%matplotlib inline

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

import metpy.calc as mpcalc
from metpy.cbook import get_test_data
from metpy.interpolate import cross_section

# NAM 218 AWIPS Grid - CONUS 
# (12-km Resolution; full complement of pressure level fields and some surface-based fields)
# https://nomads.ncep.noaa.gov/pub/data/nccf/com/nam/prod/nam.20210925/

filename = '/NAM/nam.t00z.awphys00.tm00.grib2'

#   NCEP GRIB files have non-unique keys, so pick from recommended list to filter
#    --Temperature and Pressure levels
#   Non-CF Compliant so don't use .parse_cf()

data = xr.open_dataset(filename, engine="cfgrib",filter_by_keys={'typeOfLevel': 'isobaricInhPa','name':'Temperature'})

#  We are missing a mapping projection, so add it:
data1=data.metpy.assign_crs(
    grid_mapping_name='lambert_conformal_conic',
    latitude_of_projection_origin=35,
    longitude_of_central_meridian=-100,
    standard_parallel=(30,60),
    earth_radius=6371229.0)

print(data1)

#--------------------------  Print Output  ---------------------------------
#<xarray.Dataset>
#Dimensions:        (isobaricInhPa: 39, y: 428, x: 614)
#Coordinates:
#    time           datetime64[ns] 2021-09-25
#    step           timedelta64[ns] 00:00:00
#  * isobaricInhPa  (isobaricInhPa) float64 1e+03 975.0 950.0 ... 75.0 50.0
#    latitude       (y, x) float64 ...
#    longitude      (y, x) float64 ...
#    valid_time     datetime64[ns] 2021-09-25
#    metpy_crs      object Projection: lambert_conformal_conic
#Dimensions without coordinates: y, x
#Data variables:
#    t              (isobaricInhPa, y, x) float32 ...
#Attributes:
#    GRIB_edition:            2
#    GRIB_centre:             kwbc
#    GRIB_centreDescription:  US National Weather Service - NCEP 
#    GRIB_subCentre:          0
#    Conventions:             CF-1.7
#    institution:             US National Weather Service - NCEP 
#    history:                 2021-09-28T22:45 GRIB to CDM+CF via cfgrib-0.9.9...
#  Missing x- and y-coordinates
#----------------------------------------------------------------------

#  Try to add coordinates for y and x
data2=data1.metpy.assign_y_x()

#  Output Error:
# ValueError: Projected y and x coordinates cannot be collapsed to 1D within tolerance. Verify that your latitude and longitude coordinates correspond to your CRS coordinate.'''

Answer 1

So the problem here is pointed to by the error message you included:

Verify that your latitude and longitude coordinates correspond to your CRS coordinate.

With the provided CRS parameters, the calculated x/y coordinates don't seem to be 1D values, which makes them problematic for use, and indicate that the projection parameters in general aren't quite right.

I'm not sure where the parameters you are using in assign_crs came from, but they're not correct for the NAM 218 grid. The parameters you want are:

  grid_mapping_name = "lambert_conformal_conic",
  latitude_of_projection_origin = 25.0,
  longitude_of_central_meridian = 265.0,
  standard_parallel = 25.0,
  earth_radius = 6371229.0

The most robust way to get these is to read a message using pygrib. Assuming all GRIB messages in a file are from the same coordinate system (likely and is the case for those NWS files) you could do this:

import pygrib
from pyproj import CRS

grib = pygrib.open(filename)
msg = grib.message(1)

data1 = data.metpy.assign_crs(CRS(msg.projparams).to_cf())
data2 = data1.metpy.assign_y_x()

The message from pygrib provides the appropriate pyproj coordinates, which can readily be mapped to the requisite Cf parameters.

It would be great if cfgrib could provide these parameters natively, but that's not currently the case. There's an open issue to add this.