Errors when using metpy to calculate isentropic potential vorticity

[tik619]

I am trying to calculate isentropic vorticity using the metpy function "mpcalc.potential_vorticity_baroclinic", but I am encountering an error, which says:

Traceback (most recent call last):

  File ~/miniconda3/lib/python3.8/site-packages/spyder_kernels/py3compat.py:356 in compat_exec
    exec(code, globals, locals)

  File /mnt/d/paper/Recurving/pv_isentropic.py:62
    pv2 = mpcalc.potential_vorticity_baroclinic(

  File ~/miniconda3/lib/python3.8/site-packages/metpy/calc/tools.py:1178 in wrapper
    return func(*bound_args.args, **bound_args.kwargs)

  File ~/miniconda3/lib/python3.8/site-packages/metpy/xarray.py:1328 in wrapper
    result = func(*bound_args.args, **bound_args.kwargs)

  File ~/miniconda3/lib/python3.8/site-packages/metpy/units.py:325 in wrapper
    return func(*args, **kwargs)

  File ~/miniconda3/lib/python3.8/site-packages/metpy/calc/kinematics.py:1042 in potential_vorticity_baroclinic
    np.shape(potential_temperature)[vertical_dim] < 3

IndexError: tuple index out of range.

Can anyone help me how to sove this issue. For your convenience I am also giving the code which I have prepared:

# Import libraries
import numpy as np
import xarray as xr
import metpy.calc as mpcalc
from metpy.units import units
# Load data
d1 = xr.open_dataset('/mnt/h/Recurving_data/AS/RECURVING_RIGHT/ID3/temp.nc').sel(longitude=slice(50,80),latitude=slice(0,30))
d2 = xr.open_dataset('/mnt/h/Recurving_data/AS/RECURVING_RIGHT/ID3/u_wind.nc').sel(longitude=slice(50,80),latitude=slice(0,30))
d3 = xr.open_dataset('/mnt/h/Recurving_data/AS/RECURVING_RIGHT/ID3/v_wind.nc').sel(longitude=slice(50,80),latitude=slice(0,30))
d4 = xr.open_dataset('/mnt/h/Recurving_data/AS/RECURVING_RIGHT/ID3/rh.nc').sel(longitude=slice(50,80),latitude=slice(0,30))
time = d1['time'].values
lat = d1['latitude'].values
lon = d1['longitude'].values
lev = d1['plevel'].values * units.hPa # assuming pressure levels in hPa
# Define isentropic level and PV constant
isen_level = np.array([320]) * units.K
pv_constant = 1.5e-6 * units.m**2 / units.s / units.K / units.kg
# Create an empty list to store the results
pv_list = []
# Loop over the time dimension
for t in range(len(time)):
    # Extract the data for the current time step
    T = d1['TMP_prl'].isel(time=t).values * units.K # assuming temperature in K
    u = d2['UGRD_prl'].isel(time=t).values * units.m / units.s # assuming u-wind in m/s
    v = d3['VGRD_prl'].isel(time=t).values * units.m / units.s # assuming v-wind in m/s
    rh = d4['RH_prl'].isel(time=t).values # assuming relative humidity in fraction
    # Interpolate data to isentropic level
    isen_press, isen_u, isen_v, isen_T = mpcalc.isentropic_interpolation(
        isen_level, lev, T, u, v, rh, tmpk_out=True)
    # Calculate absolute vorticity
    f = mpcalc.coriolis_parameter(lat * units.degrees)
    dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)
    # vort = mpcalc.vorticity(isen_u[0], isen_v[0], dx=dx, dy=dy) + f
    # vort = mpcalc.vorticity(isen_u[0], isen_v[0], dx, dy) + f
    vort = mpcalc.vorticity(isen_u[0], isen_v[0], dx=dx, dy=dy) + f[:, None]
    # Calculate potential temperature
    isen_theta = mpcalc.potential_temperature(isen_press[0], isen_T[0]* units.kelvin)
   #######################################    
   # Calculate PV using metpy.potential_vorticity_baroclinic
    pv2 = mpcalc.potential_vorticity_baroclinic(
    isen_theta, isen_press[0],isen_u[0], isen_v[0], dx, dy,
    latitude=lat * units.degrees)
    # Mask the values below the PV constant
    pv_masked = np.ma.masked_less(pv2, pv_constant)
    # Append the result to the list
    pv_list.append(pv_masked)
# Convert the list to a numpy array
pv_array = np.array(pv_list)
# Save the result to a new netcdf file
 result = xr.Dataset({'PV': (['time', 'lat', 'lon'], pv_array)})
  result.to_netcdf('your_result.nc')

[Z-Richard]

Hi @tik619, it would be helpful if you can share a subset of data and wrap your code posted above with ``` so that we can more easily look into the issue. Thanks!

[tik619]

Please find the google drive link attached for the files which I used for the calculations, along with the code.
data

[Z-Richard]

Interesting. Based on your code, I think you are interpolating the pressure-level data to the 320K isentropic level, but potential_vorticity_baroclinic only provides functionality to calculate potential vorticity at each pressure level, e.g.,

$$PV = -g \left(\frac{\partial u}{\partial p}\frac{\partial \theta}{\partial y} - \frac{\partial v}{\partial p}\frac{\partial \theta}{\partial x} + \frac{\partial \theta}{\partial p}(\zeta + f) \right)$$

Note that the ${\partial p}$ indicates that the input should have an additional pressure dimension, whereas currently you are supplying a 2-D array to the potential_vorticity_baroclinic function. If you need potential vorticity at the isentropic coordinates, my suggestion is to compute them using pressure coordinates before interpolating them to the isentropic coordinates. If you want something different, there are many preliminary functions such as first_derivative and second_derivative in MetPy that might suit your needs. I hope this helps!

[tik619]

Actually i want the potential vorticity at 320k isentropic level. Can you help me with that by correcting the code.

…

On Fri, 11 Aug, 2023, 18:41 Richard Zhuang, ***@***.***> wrote: Interesting. Based on your code, I think you are interpolating the pressure-level data to the 320K isentropic level, but potential_vorticity_baroclinic only provides functionality to calculate potential vorticity at each pressure level, e.g., $$PV = -g \left(\frac{\partial u}{\partial p}\frac{\partial \theta}{\partial y} - \frac{\partial v}{\partial p}\frac{\partial \theta}{\partial x} + \frac{\partial \theta}{\partial p}(\zeta + f) \right)$$ Note that the ${\partial p}$ indicates that the input should have an additional pressure dimension, whereas currently you are supplying a 2-D array to the potential_vorticity_baroclinic function. If you need potential vorticity at the isentropic coordinates, my suggestion is to compute them using pressure coordinates before interpolating them to the isentropic coordinates. If you want something different, there are many preliminary functions such as first_derivative and second_derivative in MetPy that might suit your needs. I hope this helps! — Reply to this email directly, view it on GitHub <#3143 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ASYYLPN3IS5AXQT6X5CMXG3XUYVQVANCNFSM6AAAAAA3HI3MRA> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[Z-Richard]

Well, you are not far away from what you want. And I should correct my earlier comment: you can use potential_vorticity_baroclinic in isentropic coordinate. What you need to do is to just interpolate the data to more isentropic levels, say 318, 320, 322K instead of 320K so that the function can compute a derivative in the vertical direction.

[tik619]

Can you say how to do that, like by modifying my code?

…

On Fri, 11 Aug, 2023, 19:11 Richard Zhuang, ***@***.***> wrote: Well, you are not far away from what you want. And I should correct my earlier comment: you can use potential_vorticity_baroclinic in isentropic coordinate. What you need to do is to just interpolate the data to more isentropic levels, say 318, 320, 322K instead of 320K so that the function can compute a derivative in the vertical direction. — Reply to this email directly, view it on GitHub <#3143 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ASYYLPML2E3SIQJM2XKFYJLXUYZBTANCNFSM6AAAAAA3HI3MRA> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[sgdecker]

You can try modifying line 16 to:

isen_level = np.array([318, 320, 322]) * units.K