[Bug]: ValueError when using MetPy to calculate geostrophic winds

[starforge]

What happened?

I'm trying to use the MetPy function mpcalc.geostrophic_wind() to calculate geostrophic winds from a wrfout file.

I'm getting "ValueError: Must provide dx/dy arguments or input DataArray with latitude/longitude coordinates", along with a warning, "warnings.warn('More than one ' + axis + ' coordinate present for variable'".

I don't know what's causing the problem.

Minimal Complete Verifiable Example

import metpy.calc as mpcalc
import xarray as xr
import xwrf

# Open the NetCDF file
filename = "wrfout_d01_2016-10-04_12:00:00"
ds = xr.open_dataset(filename).xwrf.postprocess()

# Extract the geopotential height
z = ds['geopotential_height']

# Compute the geostrophic wind
geo_wind_u, geo_wind_v = mpcalc.geostrophic_wind(z)

Relevant log output

/mnt/iusers01/fatpou01/sees01/w34926hb/.conda/envs/metpy_env/lib/python3.9/site-packages/metpy/xarray.py:355: UserWarning: More than one latitude coordinate present for variable "geopotential_height".
  warnings.warn('More than one ' + axis + ' coordinate present for variable'
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/mnt/iusers01/fatpou01/sees01/w34926hb/.conda/envs/metpy_env/lib/python3.9/site-packages/metpy/xarray.py", line 1508, in wrapper
    raise ValueError('Must provide dx/dy arguments or input DataArray with '
ValueError: Must provide dx/dy arguments or input DataArray with latitude/longitude coordinates.

Environment

System Information
------------------
xWRF commit : None
python      : 3.9.13 | packaged by conda-forge | (main, May 27 2022, 16:58:50)
[GCC 10.3.0]
python-bits : 64
OS          : Linux
OS-release  : 3.10.0-1127.19.1.el7.x86_64
machine     : x86_64
processor   : x86_64
byteorder   : little
LC_ALL      : None
LANG        : en_GB.UTF-8
LOCALE      : ('en_GB', 'UTF-8')

Installed Python Packages
-------------------------
cf_xarray   : 0.7.5
dask        : 2022.11.0
donfig      : 0.7.0
matplotlib  : 3.6.2
metpy       : 1.3.1
netCDF4     : 1.6.2
numpy       : 1.23.5
pandas      : 1.5.1
pint        : 0.20.1
pooch       : v1.6.0
pyproj      : 3.4.0
xarray      : 2022.11.0
xgcm        : 0.8.0
xwrf        : 0.0.2

Anything else we need to know?

No response

[lpilz]

Hi @starforge, thanks for raising this issue. We definitely want to support interoperability with metpy, so this is important to us.

I tried reproducing your issue, however unfortunately I did not manage to do so. With my own wrfout files, everything works fine. Could you maybe post print(z), so that we can see whether the variable is malformed somehow? Mine looks like this:

xarray.DataArray 'geopotential_height' (Time: 169, z_stag: 42, y: 170, x: 155)>
[...]
Coordinates:
    XLAT     (y, x) float32 38.8 38.82 38.84 38.86 ... 60.74 60.7 60.66 60.63
    XLONG    (y, x) float32 -3.86 -3.693 -3.525 -3.357 ... 30.41 30.67 30.93
    XTIME    (Time) timedelta64[ns] 00:00:00 01:00:00 ... 7 days 00:00:00
  * z_stag   (z_stag) float32 1.0 0.9981 0.9957 0.9926 ... 0.01458 0.006846 0.0
  * Time     (Time) datetime64[ns] 2020-04-01 2020-04-01T01:00:00 ... 2020-04-08
  * y        (y) float64 -1.267e+06 -1.252e+06 ... 1.253e+06 1.268e+06
  * x        (x) float64 -1.086e+06 -1.071e+06 ... 1.209e+06 1.224e+06
Attributes:
    units:          m
    standard_name:  geopotential_height
    stagger:        Z
    grid_mapping:   wrf_projection

[starforge]

Hi @lpilz , Thanks very much for your reply! Here's what print(z) returns:

<xarray.DataArray 'geopotential_height' (Time: 1, z_stag: 80, y: 1366, x: 1366)>
array([...],
dtype=float32)
Coordinates:
XLAT (y, x) float32 25.65 25.66 25.67 25.68 ... 67.97 67.94 67.91 67.89
XLONG (y, x) float32 -44.18 -44.14 -44.11 -44.07 ... 53.11 53.16 53.22
XTIME (Time) timedelta64[ns] 3 days
CLAT (y, x) float32 25.65 25.66 25.67 25.68 ... 67.97 67.94 67.91 67.89
* z_stag (z_stag) float32 1.0 0.9938 0.9871 0.9797 ... 0.004977 0.002432 0.0
* Time (Time) datetime64[ns] 2016-10-04T12:00:00
* y (y) float64 -2.73e+06 -2.726e+06 -2.722e+06 ... 2.726e+06 2.73e+06
* x (x) float64 -2.73e+06 -2.726e+06 -2.722e+06 ... 2.726e+06 2.73e+06
Attributes:
units: m
standard_name: geopotential_height
stagger: Z
grid_mapping: wrf_projection

(Sorry, I couldn't figure out how to format it properly.) I notice that my z array has a coordinate called CLAT, which yours does not have.

[lpilz]

Hi @starforge,
thanks for the output and sorry for the delay. I don't see anything obviously wrong with your DataArray. The next step would be to look at the wrfout file. Could you share that, please?

PS: you can format code into a block by using triple ticks (```) before and after the relevant section.

[lpilz]

Closing because stale

[spacekace]

Closing because stale

Can this issue be reopened? I am trying to use xWRF to open a CONUS404 NetCDF file and compute some variables using MetPy. I use the same code snippet as in the original post and receive the same warning. However, the traceback information is different. I have tried a number of things to fix this issue, including using wrf-python with NetCDF4 instead of xWRF, passing dx/dy/lat arrays to the mpcalc function, etc., amongst other things. I am at a loss for opening CONUS404 (ideally as an xarray for functionality with some additional calculations) and using MetPy to calculate geostrophic wind, inertial advective wind, etc. Would prefer to use this functionality rather than having to start over and compute everything from scratch.

Thank you!!

Code:

import metpy.calc as mpcalc
import xarray as xr
import xwrf

# Open the NetCDF file
filename = "/glade/work/kshourd/29_June_2012/wrfout_d01_2012-06-29_12:00:00"
ds = xr.open_dataset(filename).xwrf.postprocess()

# Extract the geopotential height and wind variables
z = ds['geopotential_height']
print(z)

ua = ds['wind_east']
va = ds['wind_north']

# Smooth height data
#z = mpcalc.smooth_gaussian(z, 3)

# Compute the geostrophic wind
geo_wind_u, geo_wind_v = mpcalc.geostrophic_wind(z)

# Calculate ageostrophic wind components
ageo_wind_u = ua - geo_wind_u
ageo_wind_v = va - geo_wind_v

Log:

/glade/work/kshourd/conda-envs/wrfplot/lib/python3.9/site-packages/metpy/xarray.py:382: UserWarning: More than one latitude coordinate present for variable "geopotential_height".
  warnings.warn('More than one ' + axis + ' coordinate present for variable'
Traceback (most recent call last):
  File "/glade/work/kshourd/SCREAMING.py", line 23, in <module>
    geo_wind_u, geo_wind_v = mpcalc.geostrophic_wind(z)
  File "/glade/work/kshourd/conda-envs/wrfplot/lib/python3.9/site-packages/metpy/calc/tools.py", line 1178, in wrapper
    return func(*bound_args.args, **bound_args.kwargs)
  File "/glade/work/kshourd/conda-envs/wrfplot/lib/python3.9/site-packages/metpy/xarray.py", line 1329, in wrapper
    result = func(*bound_args.args, **bound_args.kwargs)
  File "/glade/work/kshourd/conda-envs/wrfplot/lib/python3.9/site-packages/metpy/units.py", line 320, in wrapper
    return func(*args, **kwargs)
  File "/glade/work/kshourd/conda-envs/wrfplot/lib/python3.9/site-packages/metpy/calc/kinematics.py", line 620, in geostrophic_wind
    f = coriolis_parameter(latitude)
  File "/glade/work/kshourd/conda-envs/wrfplot/lib/python3.9/site-packages/metpy/xarray.py", line 1329, in wrapper
    result = func(*bound_args.args, **bound_args.kwargs)
  File "/glade/work/kshourd/conda-envs/wrfplot/lib/python3.9/site-packages/metpy/calc/basic.py", line 609, in coriolis_parameter
    latitude = _check_radians(latitude, max_radians=np.pi / 2)
  File "/glade/work/kshourd/conda-envs/wrfplot/lib/python3.9/site-packages/metpy/calc/basic.py", line 1296, in _check_radians
    if np.any(np.greater(np.abs(value), max_radians)):
TypeError: bad operand type for abs(): 'NoneType'

Here is the result from printing z. Again, nothing obviously unusual:

Coordinates:
    XLAT     (y, x) float32 17.65 17.66 17.67 17.68 ... 51.74 51.73 51.71 51.69
    XLONG    (y, x) float32 -122.6 -122.5 -122.5 -122.5 ... -57.17 -57.12 -57.07
    XTIME    (Time) timedelta64[ns] 11960 days 12:00:00
    CLAT     (y, x) float32 17.65 17.66 17.67 17.68 ... 51.74 51.73 51.71 51.69
  * z_stag   (z_stag) float32 1.0 0.9938 0.9864 0.9779 ... 0.01243 0.0056 0.0
  * Time     (Time) datetime64[ns] 2012-06-29T12:00:00
  * x        (x) float64 -2.732e+06 -2.728e+06 ... 2.728e+06 2.732e+06
  * y        (y) float64 -2.028e+06 -2.024e+06 -2.02e+06 ... 2.024e+06 2.028e+06
Attributes:
    units:          m
    standard_name:  geopotential_height
    stagger:        Z
    grid_mapping:   wrf_projection

[lpilz]

It works with my data. Moving this to discussions (#152) until we figure out whether there is a problem with xwrf. If we do, I'm more than happy to reopen.

[spacekace]

So perhaps this is an issue with the CONUS404 WRF data and not the xwrf/metpy interface in general? Just trying to figure out how to work with this data...

…

On Thu, Dec 7, 2023, 1:48 AM Lukas Pilz ***@***.***> wrote: It works with my data. Moving this to discussions until we figure out whether there is a problem with xwrf. If we do, I'm more than happy to reopen. — Reply to this email directly, view it on GitHub <#120 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AGJ6SGGY333JCRLYIIULMCLYIGGFDAVCNFSM6AAAAAASSZ66C2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQNBVGAYTSNRTGU> . You are receiving this because you commented.Message ID: ***@***.***>