Parse CF Grid Mapping

[blaylockbk]

Let me know if I'm missing something, but it doesn't appear that cfgrib parses grid mappings details as defined in the CF convention (see Appendix F: Grid Mappings)

I suspect that the CF "grid_mapping" attribute could be derived from GRIB_gridType and related attributes.

[blaylockbk]

Just jotting down some ideas to parse the CF grid_mapping from the variable attributes loaded with cfgrib...

The CF grid_mapping parameters can be parsed for lambert and regular_ll grid types. If we assume the grid_mapping is the same for all variables in an xarray Dataset, the attributes from one variable can be used to build the grid_mapping.

In the example below, ds is an xarray Dataset that opened a GRIB file with cfgrib...

# Attach CF grid mapping
# ----------------------
# See http://cfconventions.org/Data/cf-conventions/cf-conventions-1.8/cf-conventions.html#appendix-grid-mappings

# Create the dummy variable to store the grid mapping information
grid_mapping = "gribfile_projection"
ds[grid_mapping] = None
ds[grid_mapping].attrs["long_name"] = f"gribfile grid projection"

# To get the grid_mapping parameters, we need attributes from a reference variable, the first will do
attrs = ds[list(ds)[0]].attrs
if attrs["GRIB_gridType"] == "lambert":
    ds[grid_mapping].attrs["grid_mapping_name"] = "lambert_conformal_conic"
    ds[grid_mapping].attrs["standard_parallel"] = (
        attrs.get("GRIB_Latin1InDegrees"),
        attrs.get("GRIB_Latin2InDegrees"),
    )
    ds[grid_mapping].attrs["longitude_of_central_meridian"] = attrs.get( "GRIB_LoVInDegrees")
    ds[grid_mapping].attrs["latitude_of_projection_origin"] = attrs.get("GRIB_LaDInDegrees")
    ds[grid_mapping].attrs["false_easting"] = 0
    ds[grid_mapping].attrs["false_northing"] = 0

elif attrs["GRIB_gridType"] == "regular_ll":
    ds[grid_mapping].attrs["grid_mapping_name"] = "latitude_longitude"
    # Nothing else needed for this as far as I can tell

elif attrs["GRIB_gridType"] == "polar_stereographic":
    # The GRIB file I tested with polar_stereographic coordinates didn't have enough info
    ds[grid_mapping].attrs["grid_mapping_name"] = "polar_stereographic"
    ds[grid_mapping].attrs["straight_vertical_longitude_from_pole"] = "INFO NOT IN GRIB FILE"
    ds[grid_mapping].attrs["latitude_of_projection_origin"] = "INFO NOT IN GRIB FILE"
    ds[grid_mapping].attrs["standard_parallel"] = "INFO NOT IN GRIB FILE"
    ds[grid_mapping].attrs["false_easting"] = 0
    ds[grid_mapping].attrs["false_northing"] = 0

else:
    print("GRIB_gridType not recognized")

# Set this grid_mapping for all variables in the Dataset
for var in list(ds):
    ds[var].attrs["grid_mapping"] = grid_mapping

One benefit of including the grid_mapping CF data is that the projection information can be understood by tools like MetPy's parse_cf to create a cartopy map projection (see here).

[dopplershift]

I'm not sure how to do the GRIB-handling part with cfgrib/eccodes, but pygrib can handle giving you the appropriate pyproj parameters already, and pyproj already handles mapping to CF:

import pygrib
from pyproj import CRS

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

cf_params = CRS(msg.projparams).to_cf()

[blaylockbk]

Thanks @dopplershift! That drastically simplifies what I was trying to do. I love it when my thirty lines of code can be done in two.