Necessary horizontal dimension coordinates are missing, causing v1.4 (and later) kinematics calculations to no longer work with latitude/longitude coordinates alone #2960
Labels
Area: Calc
Pertains to calculations
Area: Xarray
Pertains to xarray integration
Type: Question
Issues does not require work, but answers a user question
Comments
|
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 |
|
yes my apologies, it's a grib2 file from ncep/nomads grib filter... |
|
preview: |
|
bug is still present in 1.5.0 |
(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 MetPy's calculations involving geospatial grid deltas require horizontal dimension coordinates (i.e., 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)With this, you should be able to compute vorticity without concern: vortData = mpcalc.vorticity(grid.u, grid.v)
print(vortData) |
jthielen
added
Type: Question
jthielen
changed the title
wx4stg
added a commit
to wx4stg/hdwx-modelplotter
that referenced
this issue
Jun 7, 2023
|
Thanks for the help @jthielen |
|
Given the missing metadata, I'm surprised MetPy 1.3 didn't choke on this too. Must have been subtly more forgiving... |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
in metpy 1.3.1 this runs fine:
on 1.4 this produces the following traceback:
The text was updated successfully, but these errors were encountered: