MetPy
MetPy copied to clipboard
Unidata
Necessary horizontal dimension coordinates are missing, causing v1.4 (and later) kinematics calculations to no longer work with latitude/longitude coordinates alone
in metpy 1.3.1 this runs fine:
modelDataArray = xr.open_dataset("winds.grib2", engine="cfgrib").sel(isobaricInhPa=500)
uwind = modelDataArray.u
uwind = uwind.metpy.quantify()
uwind = uwind.metpy.convert_units("m/s")
vwind = modelDataArray.v
vwind = vwind.metpy.quantify()
vwind = vwind.metpy.convert_units("m/s")
vortData = mpcalc.vorticity(uwind, vwind)
on 1.4 this produces the following traceback:
File "/home/stgardner4/hdwx-operational/hdwx-modelplotter/modelPlot.py", line 583, in vort500Plot
vortData = mpcalc.vorticity(uwind, vwind)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/stgardner4/micromamba/envs/HDWX/lib/python3.11/site-packages/metpy/calc/tools.py", line 1094, in wrapper
grid_deltas = grid_prototype.metpy.grid_deltas
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/stgardner4/micromamba/envs/HDWX/lib/python3.11/site-packages/metpy/xarray.py", line 544, in grid_deltas
raise AttributeError(
AttributeError: Grid deltas cannot be calculated since horizontal dimension coordinates cannot be found.. Did you mean: 'time_deltas'?
Can you share the data files you're working with or at least a preview of the dataarrays? I'd be curious to see why these worked before but not now.
Either way, you'll need identifiable x, y (or lat, lon) coordinates for these functions to work without specifying any other information. You can get around this by either renameing your coordinates to match CF Standard Names (or some "close enough" names that MetPy can understand), or by manually telling MetPy explicitly how your coordinates map via .metpy.assign_coordinates. Or you can supply dx, dy and skip xarray handling, though I'd recommend against that here.
preview:
<xarray.Dataset>
Dimensions: (isobaricInhPa: 3, y: 311, x: 564)
Coordinates:
time datetime64[ns] ...
step timedelta64[ns] ...
* isobaricInhPa (isobaricInhPa) float64 850.0 500.0 250.0
latitude (y, x) float64 ...
longitude (y, x) float64 ...
valid_time datetime64[ns] ...
Dimensions without coordinates: y, x
Data variables:
u (isobaricInhPa, y, x) float32 ...
v (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: 2023-03-17T21:24 GRIB to CDM+CF via cfgrib-0.9.1
(HDWX-dev) stgardner4@wx4stg:~$ curl -L "https://nomads.ncep.noaa.gov/cgi-bin/filter_hrrr_2d.pl?dir=%2Fhrrr.20230605%2Fconus&file=hrrr.t00z.wrfsfcf00.grib2&var_UGRD=on&var_VGRD=on&lev_850_mb=on&lev_500_mb=on&lev_250_mb=on" -o winds.grib2
% Total % Received % Xferd Average Speed Time Time Time Current
Dload Upload Total Spent Left Speed
100 4119k 0 4119k 0 0 14.0M 0 --:--:-- --:--:-- --:--:-- 14.0M
(HDWX-dev) stgardner4@wx4stg:~$ python3
Python 3.11.3 | packaged by conda-forge | (main, Apr 6 2023, 08:57:19) [GCC 11.3.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import xarray as xr
>>> import metpy
>>> from metpy import calc as mpcalc
>>> modelDataArray = xr.open_dataset("winds.grib2", engine="cfgrib").sel(isobaricInhPa=500)
>>> uwind = modelDataArray.u
>>> uwind = uwind.metpy.quantify()
>>> uwind = uwind.metpy.convert_units("m/s")
>>> vwind = modelDataArray.v
>>> vwind = vwind.metpy.quantify()
>>> vwind = vwind.metpy.convert_units("m/s")
>>> uwind
<xarray.DataArray 'u' (y: 1059, x: 1799)>
<Quantity([[11.928843 11.928843 11.928843 ... 5.9913425 5.9288425
5.8663425 ]
[11.991343 11.991343 11.991343 ... 6.1788425 6.1163425
6.0538425 ]
[11.991343 11.991343 11.991343 ... 6.3038425 6.3038425
6.2413425 ]
...
[12.553843 12.616343 12.741343 ... 0.74134254 0.74134254
0.74134254]
[12.616343 12.741343 12.803843 ... 0.80384254 0.80384254
0.80384254]
[12.678843 12.803843 12.928843 ... 0.92884254 0.92884254
0.86634254]], 'meter / second')>
Coordinates:
time datetime64[ns] ...
step timedelta64[ns] ...
isobaricInhPa float64 500.0
latitude (y, x) float64 ...
longitude (y, x) float64 ...
valid_time datetime64[ns] ...
Dimensions without coordinates: y, x
Attributes: (12/32)
GRIB_paramId: 131
GRIB_dataType: fc
GRIB_numberOfPoints: 1905141
GRIB_typeOfLevel: isobaricInhPa
GRIB_stepUnits: 1
GRIB_stepType: instant
... ...
GRIB_missingValue: 3.4028234663852886e+38
GRIB_name: U component of wind
GRIB_shortName: u
GRIB_units: m s**-1
long_name: U component of wind
standard_name: eastward_wind
>>> vwind
<xarray.DataArray 'v' (y: 1059, x: 1799)>
<Quantity([[ 3.8059063 3.8059063 3.8059063 ... 10.430906 10.493406
10.618406 ]
[ 3.8684063 3.8684063 3.8684063 ... 10.430906 10.493406
10.618406 ]
[ 3.8684063 3.8684063 3.8684063 ... 10.430906 10.493406
10.618406 ]
...
[-13.569094 -13.631594 -13.756594 ... 4.7434063 4.6809063
4.6809063]
[-13.694094 -13.756594 -13.881594 ... 4.6809063 4.6184063
4.6184063]
[-13.819094 -13.881594 -14.006594 ... 4.6184063 4.5559063
4.5559063]], 'meter / second')>
Coordinates:
time datetime64[ns] ...
step timedelta64[ns] ...
isobaricInhPa float64 500.0
latitude (y, x) float64 ...
longitude (y, x) float64 ...
valid_time datetime64[ns] ...
Dimensions without coordinates: y, x
Attributes: (12/32)
GRIB_paramId: 132
GRIB_dataType: fc
GRIB_numberOfPoints: 1905141
GRIB_typeOfLevel: isobaricInhPa
GRIB_stepUnits: 1
GRIB_stepType: instant
... ...
GRIB_missingValue: 3.4028234663852886e+38
GRIB_name: V component of wind
GRIB_shortName: v
GRIB_units: m s**-1
long_name: V component of wind
standard_name: northward_wind
>>> avort = mpcalc.absolute_vorticity(uwind, vwind)
<stdin>:1: UserWarning: More than one time coordinate present for variable "u".
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/home/stgardner4/metpy/src/metpy/calc/tools.py", line 1094, in wrapper
grid_deltas = grid_prototype.metpy.grid_deltas
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/stgardner4/metpy/src/metpy/xarray.py", line 543, in grid_deltas
raise AttributeError(
AttributeError: Grid deltas cannot be calculated since horizontal dimension coordinates cannot be found.. Did you mean: 'time_deltas'?
by re-ordering the if-elif in tools.py so that lines 1097 through 1102 are preferred to 1093 through 1096, like so:
if longitude is not None and latitude is not None and crs is not None:
# TODO: de-duplicate .metpy.grid_deltas code
geod = None if crs is None else crs.get_geod()
bound_args.arguments['dx'], bound_args.arguments['dy'] = (
nominal_lat_lon_grid_deltas(longitude, latitude, geod)
)
elif grid_prototype is not None:
grid_deltas = grid_prototype.metpy.grid_deltas
bound_args.arguments['dx'] = grid_deltas['dx']
bound_args.arguments['dy'] = grid_deltas['dy']
...
and then specifying latitude, longitude, and crs in my code: (run following example in previous comment)
>>> from pyproj import CRS
>>> from metpy.units import units
>>> avort = mpcalc.absolute_vorticity(uwind, vwind, latitude=uwind.latitude * units.degree_north, longitude=uwind.longitude * units.degree_east, crs=CRS('+proj=latlon'))
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/home/stgardner4/metpy/src/metpy/calc/tools.py", line 1097, in wrapper
nominal_lat_lon_grid_deltas(longitude, latitude, geod)
File "/home/stgardner4/metpy/src/metpy/xarray.py", line 1328, in wrapper
result = func(*bound_args.args, **bound_args.kwargs)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/stgardner4/metpy/src/metpy/calc/tools.py", line 875, in nominal_lat_lon_grid_deltas
raise ValueError(
ValueError: Cannot calculate nominal grid spacing from longitude and latitude arguments that are not one dimensional.
As noted in the initial error message
AttributeError: Grid deltas cannot be calculated since horizontal dimension coordinates cannot be found.
MetPy's calculations involving geospatial grid deltas require horizontal dimension coordinates (i.e., y and x coordinates, not just latitude and longitude). Since your dataset from cfgrib includes latitude and longitude auxiliary coordinates, but not a CF-compliant grid mapping or the necessary y and x coordinates, we can add them as follows (which is a bit messier than we'd like it to be at this point, see https://github.com/Unidata/MetPy/issues/2473 for thoughts on making this easier):
grid = modelDataArray.metpy.assign_crs({
"semi_major_axis": 6371200.0,
"semi_minor_axis": 6371200.0,
"grid_mapping_name": "lambert_conformal_conic",
"standard_parallel": [
modelDataArray["u"].attrs["GRIB_Latin1InDegrees"],
modelDataArray["u"].attrs["GRIB_Latin2InDegrees"]
],
"latitude_of_projection_origin": modelDataArray["u"].attrs["GRIB_LaDInDegrees"],
"longitude_of_central_meridian": modelDataArray["u"].attrs["GRIB_LoVInDegrees"],
}).metpy.assign_y_x()
print(grid)
<xarray.Dataset>
Dimensions: (y: 311, x: 564)
Coordinates:
time datetime64[ns] ...
step timedelta64[ns] ...
isobaricInhPa float64 500.0
latitude (y, x) float64 16.0 16.02 16.05 16.08 ... 48.88 48.85 48.82
longitude (y, x) float64 231.2 231.3 231.4 231.6 ... 305.6 305.7 305.9
valid_time datetime64[ns] ...
metpy_crs object Projection: lambert_conformal_conic
* y (y) float64 -5.523e+05 -5.401e+05 ... 3.215e+06 3.227e+06
* x (x) float64 -3.617e+06 -3.604e+06 ... 3.235e+06 3.247e+06
Data variables:
u (y, x) float32 -1.281 -1.261 -1.251 ... -4.651 -4.181 -3.801
v (y, x) float32 -4.467 -4.467 -4.477 ... 8.623 7.923 7.573
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: 2023-06-05T05:40 GRIB to CDM+CF via cfgrib-0.9.1...
With this, you should be able to compute vorticity without concern:
vortData = mpcalc.vorticity(grid.u, grid.v)
print(vortData)
<xarray.DataArray (y: 311, x: 564)>
<Quantity([[-1.99231253e-05 -1.53559366e-05 -1.32796818e-05 ... -4.19774231e-05
-4.69511698e-05 -5.19246961e-05]
[-1.90933695e-05 -1.70162675e-05 -1.32799313e-05 ... -4.36305248e-05
-4.69465399e-05 -5.02620898e-05]
[-1.74338684e-05 -1.82610650e-05 -1.82583878e-05 ... -4.40410323e-05
-4.36295288e-05 -4.48737385e-05]
...
[-2.91201929e-05 -2.10776388e-05 -9.89749391e-06 ... -1.09145971e-06
-3.31177852e-06 8.01820970e-06]
[-2.24088316e-05 -2.01775769e-05 -1.34678251e-05 ... -6.78989148e-05
-6.74049967e-05 -3.53058467e-05]
[-1.25554145e-05 -1.47978150e-05 -1.43550163e-05 ... -1.14015576e-04
-1.06286680e-04 -3.89474551e-05]], '1 / second')>
Coordinates:
time datetime64[ns] ...
step timedelta64[ns] ...
isobaricInhPa float64 500.0
latitude (y, x) float64 16.0 16.02 16.05 16.08 ... 48.88 48.85 48.82
longitude (y, x) float64 231.2 231.3 231.4 231.6 ... 305.6 305.7 305.9
valid_time datetime64[ns] ...
metpy_crs object Projection: lambert_conformal_conic
* y (y) float64 -5.523e+05 -5.401e+05 ... 3.215e+06 3.227e+06
* x (x) float64 -3.617e+06 -3.604e+06 ... 3.235e+06 3.247e+06
Given the missing metadata, I'm surprised MetPy 1.3 didn't choke on this too. Must have been subtly more forgiving...
