My goal is to compute frontogenesis at a given pressure level using the mpcalc.frontogenesis method.
I have read in the relevant data from NOMADS:
t7 = data['temp'].sel(lev=700.0,lat=lats,lon=lons)
u7 = data['u'].sel(lev=700.0,lat=lats,lon=lons).squeeze()*1.94384449
v7 = data['v'].sel(lev=700.0,lat=lats,lon=lons).squeeze()*1.94384449
h7 = data['gph'].sel(lev=700.0,lat=lats,lon=lons).squeeze()
Computed the relevant coordinate/grid spacing values:
x, y = t7.metpy.coordinates('x', 'y')
dx, dy = mpcalc.lat_lon_grid_deltas(lons, lats)
lat, lon = xr.broadcast(y, x)
Tagged each array with units:
t7u = t7.values * units.degC
u7u = u7.values * units.knots
v7u = v7.values * units.knots
Then I go to compute FGEN:
h7_fgen = mpcalc.frontogenesis(t7u,u7u,v7u,dx,dy,dim_order='xy')
which produces the following error:
ValueError: operands could not be broadcast together with shapes (359,219) (358,220)
The shapes of the temperature, u, and v arrays respectively are:
(220, 360)
(220, 360)
(220, 360)
So while dx and dy will understandably have lengths n-1, I'm not sure why this isn't working.
If pulling data directly from the NOMADS GFS source, it looks like your data should be in yx dimension order (which is also the order assumed by MetPy's lat_lon_grid_deltas function). And so, I would presume that something got messed up in the script with regards to dimension ordering. Resolving this, and making additional use of some of MetPy's xarray integration features to clean up the implementation, we have:
# Imports
import metpy.calc as mpcalc
from metpy.units import units
import xarray as xr
# Open data and parse for CRS
data = xr.open_dataset('http://nomads.ncep.noaa.gov/dods/gfs_0p25_1hr/gfs20201230/gfs_0p25_1hr_12z').metpy.parse_cf()
# Subset out the desired data
# Note that this leaves in time, but that can also be selected out if only
# one time is needed and you want to have smaller data arrays to compute on
t7 = data['tmpprs'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310))
u7 = data['ugrdprs'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310)).squeeze()
v7 = data['vgrdprs'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310)).squeeze()
h7 = data['hgtprs'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310)).squeeze()
# Clean up the xarray metadata for MetPy's use, since provided data is missing
# the units attribute
t7.attrs['units'] = 'K'
u7.attrs['units'] = v7.attrs['units'] = 'm/s'
h7.attrs['units'] = 'm'
# Get grid deltas
dx, dy = mpcalc.grid_deltas_from_dataarray(t7)
# Compute frontogenesis
h7_fgen = mpcalc.frontogenesis(t7, u7, v7, dx, dy)
In MetPy v1.0, the grid arguments are automatically extracted from DataArray inputs, and so, this can be simplified to
# Imports
import metpy.calc as mpcalc
from metpy.units import units
import xarray as xr
# Open data and parse for CRS
data = xr.open_dataset('http://nomads.ncep.noaa.gov/dods/gfs_0p25_1hr/gfs20201230/gfs_0p25_1hr_12z').metpy.parse_cf()
# Subset out the desired data
t7 = data['tmpprs'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310))
u7 = data['ugrdprs'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310)).squeeze()
v7 = data['vgrdprs'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310)).squeeze()
h7 = data['hgtprs'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310)).squeeze()
# Clean up the xarray metadata for MetPy's use, since provided data is missing
# the units attribute
t7.attrs['units'] = 'K'
u7.attrs['units'] = v7.attrs['units'] = 'm/s'
h7.attrs['units'] = 'm'
# Compute frontogenesis
h7_fgen = mpcalc.frontogenesis(t7, u7, v7)