convert metpy units GFS precipitation_rate

Question

I followed the example here (https://unidata.github.io/python-gallery/examples/Precipitation_Map.html) but when trying to access the units I get a dimensionality error. I'm assuming this is to do with pint and the way it parses units and the fact precipitation is a rate. Any help would be appreciated

%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import metpy
import datetime
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from xarray.backends import NetCDF4DataStore
import xarray as xr
import numpy as np
from metpy.units import masked_array, units
from siphon.catalog import TDSCatalog



best_gfs = TDSCatalog('http://thredds.ucar.edu/thredds/catalog/grib/NCEP/GFS/'
                      'Global_0p25deg/catalog.xml?dataset=grib/NCEP/GFS/Global_0p25deg/Best')
best_gfs.datasets


best_ds = list(best_gfs.datasets.values())[0]
ncss = best_ds.subset()

query = ncss.query()
query.lonlat_box(-66.243114,-25.762908,-28.708933, -2.191886).time_range(datetime.datetime.utcnow(), datetime.datetime.utcnow() + datetime.timedelta(days=10))
query.accept('netcdf4')
query.variables('Precipitation_rate_surface')
data = ncss.get_data(query)
data = xr.open_dataset(NetCDF4DataStore(data))

lon_2d, lat_2d = np.meshgrid(data['lon'], data['lat'])
precip = data['Precipitation_rate_surface']

precip.metpy.units

Answer 1

This is caused by the fact that the units string for that variable is kg.m-2.s-1, which is a UDUnits-compatible string, but does not work with the default unit parser in Pint, which is what MetPy uses for unit support.

This is fixed in MetPy 1.0. You can install the second release candidate for MetPy 1.0 with conda:

conda install -c conda-forge/label/metpy_rc metpy=1.0

or with pip:

pip install --pre metpy=1.0

A workaround for MetPy < 1.0 would be to overwrite the units metadata before calling precip.metpy.units:

precip.attrs['units'] = 'kg m^-2 s^-1'