MetPy "interpolate_to_isosurface" results in "IndexError: Unlabeled multi-dimensional array cannot be used for indexing: pressure_level"

[beuraieon-a]

Hello, newbie Python programmer here... I need help on the usage of MetPy's interpolate_to_isosurface function.

My goal was to interpolate a gridded meteorological dataset, such that its vertical coordinate is transformed from pressure-based (isobaric surfaces in hPa) to altitude-based (in kilometers above mean sea level). I used a NetCDF file containing gridded meteorological data from ERA5 that I preprocessed such that I end up with the following variables of concern: WVF_u, WVF_v and WVF_w (zonal, meridional and vertical components of isobaric water vapor flux), and height (height above sea level of isobaric surface). All of these variables have the dimensions pressure_level (vertical coordinate in terms of pressure level, i.e. pressure of isobaric surface, in hPa), latitude and longitude. I want to do the interpolation for WVF_u, WVF_v and WVF_w; while height serves to provide the necessary altitude data for the interpolation. The NetCDF file is here.

So far, my unfinished script looks like this:

import numpy as np
import xarray as xr
import metpy.calc as mpcalc
import metpy.interpolate as mpinterp
from metpy.units import units

###################
#### FUNCTIONS ####
###################

# FUNCTION 1: Data extraction
def extract_data(nc_file, variable_name, time):
    dataset = xr.open_dataset(nc_file)
    
    extracted_data = dataset[variable_name].sel(valid_time=time)
    dataset.close()

    return extracted_data

#####################################
#### MAIN PROGRAM IMPLEMENTATION ####
#####################################

# Setting the data extraction parameters
file_path = 'C:/Users/Brian/Downloads/wvf_isobaric.nc'
time = '2024-07-23T06:00:00'

# Extracting the required input data
uwind = extract_data(file_path, 'u', time)
vwind = extract_data(file_path, 'v', time)
omega = extract_data(file_path, 'w', time)
spechum = extract_data(file_path, 'q', time)
temp = extract_data(file_path, 't', time)
geopot = extract_data(file_path, 'z', time)

# Calculating altitude of pressure level
height = (mpcalc.geopotential_to_height(geopot * units('m^2/s^2'))).metpy.convert_units('km')
del geopot

# Calculating vertical component of wind (vertical velocity)
mixrat = mpcalc.mixing_ratio_from_specific_humidity(spechum * units('kg/kg'))
wwind = mpcalc.vertical_velocity(omega * units('Pa/s'),
                                 omega['pressure_level'] * units.hPa,
                                 temp * units.kelvin,
                                 mixrat)
del mixrat, omega, temp

# Calculating the components of water vapor flux
WVF_u = (spechum * units('kg/kg')) * (uwind * units('m/s'))
WVF_v = (spechum * units('kg/kg')) * (vwind * units('m/s'))
WVF_w = (spechum * units('kg/kg')) * wwind
del uwind, vwind, wwind, spechum

#####  Interpolating the dataset to an   #####
##### altitude-based vertical coordinate #####

# Defining the target altitude levels (in km)
target_altitudes = np.arange(0, 16.5, 0.5)  # 0 to 16 km in 0.5 km increments

# Getting the latitude-longitude coordinate values
lats = WVF_u['latitude'].values
lons = WVF_u['longitude'].values

# Preparing the lists to store interpolated slices
interp_u_list = []
interp_v_list = []
interp_w_list = []

# Looping through each target altitude level
# NOTE: MetPy's "interpolate_to_isosurface" function assumes that the
# highest vertical level (lowest pressure) is the zeroth index. The input dataset
# is structured such that the lowest vertical level (highest pressure) is zeroth index.
# Thus, to circumvent this, bottom_up_search is set to "False" (top-down) so that the
# function accordingly takes first the lowest pressure value (which is assigned at the
# highest index, i.e. last value).

for alt in target_altitudes:
    interp_u = mpinterp.interpolate_to_isosurface(height,
                                                  WVF_u,
                                                  alt * units.km,
                                                  bottom_up_search = False)
    interp_v = mpinterp.interpolate_to_isosurface(height,
                                                  WVF_v,
                                                  alt * units.km,
                                                  bottom_up_search = False)
    interp_w = mpinterp.interpolate_to_isosurface(height,
                                                  WVF_w,
                                                  alt * units.km,
                                                  bottom_up_search = False)
    
    interp_u_list.append(interp_u)
    interp_v_list.append(interp_v)
    interp_w_list.append(interp_w)
    
    break

# End of program

But this results in an IndexError: Unlabeled multi-dimensional array cannot be used for indexing: pressure_level, with this error message:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
C:\ProgramData\anaconda3\Lib\site-packages\spyder_kernels\customize\code_runner.py in ?(self, code, filename, ns_globals, ns_locals, post_mortem, exec_fun, capture_last_expression, global_warning)
    523                 error_type, error, tb = sys.exc_info()
    524                 self._post_mortem_excepthook(error_type, error, tb)
    525             else:
    526                 # We ignore the call to exec
--> 527                 self.shell.showtraceback(tb_offset=1)
    528         finally:
    529             __tracebackhide__ = "__pdb_exit__"

C:\ProgramData\anaconda3\Lib\site-packages\spyder_kernels\customize\utils.py in ?(code_ast, globals, locals, exec_fun, filename)
    207         if filename is None:
    208             filename = "<stdin>"
    209         exec_fun(compile(code_ast, filename, "exec"), globals, None)
    210     finally:
--> 211         if use_locals_hack:
    212             # Cleanup code
    213             globals.pop("_spyderpdb_code", None)
    214             if len(globals["_spyderpdb_locals"]) > 1:

c:\users\brian\desktop\wvf_3d.py in ?()
     74 # function accordingly takes first the lowest pressure value (which is assigned at the
     75 # highest index, i.e. last value).
     76 
     77 for alt in target_altitudes:
---> 78     interp_u = mpinterp.interpolate_to_isosurface(height,
     79                                                   WVF_u,
     80                                                   alt * units.km,
     81                                                   bottom_up_search = False)

C:\ProgramData\anaconda3\Lib\site-packages\metpy\pandas.py in ?(*args, **kwargs)
     17         # and this is more explicitly handling pandas
     18         args = tuple(a.values if isinstance(a, pd.Series) else a for a in args)
     19         kwargs = {name: (v.values if isinstance(v, pd.Series) else v)
     20                   for name, v in kwargs.items()}
---> 21         return func(*args, **kwargs)

C:\ProgramData\anaconda3\Lib\site-packages\metpy\interpolate\grid.py in ?(level_var, interp_var, level, bottom_up_search)
    348     above, below, good = find_bounding_indices(level_var, [level], axis=0,
    349                                                from_below=bottom_up_search)
    350 
    351     # Linear interpolation of variable to interpolated surface value
--> 352     interp_level = (((level - level_var[above]) / (level_var[below] - level_var[above]))
    353                     * (interp_var[below] - interp_var[above])) + interp_var[above]
    354 
    355     # Handle missing values and instances where no values for surface exist above and below

C:\ProgramData\anaconda3\Lib\site-packages\xarray\core\dataarray.py in ?(self, key)
    909         if isinstance(key, str):
    910             return self._getitem_coord(key)
    911         else:
    912             # xarray-style array indexing
--> 913             return self.isel(indexers=self._item_key_to_dict(key))

C:\ProgramData\anaconda3\Lib\site-packages\xarray\core\dataarray.py in ?(self, indexers, drop, missing_dims, **indexers_kwargs)
   1543 
   1544         indexers = either_dict_or_kwargs(indexers, indexers_kwargs, "isel")
   1545 
   1546         if any(is_fancy_indexer(idx) for idx in indexers.values()):
-> 1547             ds = self._to_temp_dataset()._isel_fancy(
   1548                 indexers, drop=drop, missing_dims=missing_dims
   1549             )
   1550             return self._from_temp_dataset(ds)

C:\ProgramData\anaconda3\Lib\site-packages\xarray\core\dataset.py in ?(self, indexers, drop, missing_dims)
   2776         *,
   2777         drop: bool,
   2778         missing_dims: ErrorOptionsWithWarn = "raise",
   2779     ) -> Self:
-> 2780         valid_indexers = dict(self._validate_indexers(indexers, missing_dims))
   2781 
   2782         variables: dict[Hashable, Variable] = {}
   2783         indexes, index_variables = isel_indexes(self.xindexes, valid_indexers)

C:\ProgramData\anaconda3\Lib\site-packages\xarray\core\dataset.py in ?(self, indexers, missing_dims)
   2556                     elif isinstance(index, CFTimeIndex):
   2557                         v = _parse_array_of_cftime_strings(v, index.date_type)
   2558 
   2559                 if v.ndim > 1:
-> 2560                     raise IndexError(
   2561                         "Unlabeled multi-dimensional array cannot be "
   2562                         f"used for indexing: {k}"
   2563                     )

IndexError: Unlabeled multi-dimensional array cannot be used for indexing: pressure_level

And this confuses me because, upon checking, the variables WVF_u, WVF_v, WVF_w and height are proper xarray DataArrays and that their dimensions and coordinates are defined. Example, for height:

<xarray.DataArray (pressure_level: 27, latitude: 85, longitude: 101)> Size: 927kB
<Quantity([[[ 3.95379309e-03  9.63872392e-03  5.75104542e-03 ...  8.68324190e-02
    8.91778395e-02  9.12300721e-02]
  [ 3.89006082e-03  3.91555391e-03 -5.83946879e-04 ...  8.54812562e-02
    8.80306214e-02  9.04652551e-02]
  [-2.22823862e-03 -4.26766928e-03 -5.36386296e-03 ...  8.45762268e-02
    8.70108679e-02  8.96367133e-02]
  ...
  [ 8.23455527e-02  8.17464516e-02  8.12748149e-02 ...  8.09689015e-02
    8.11218545e-02  8.10963660e-02]
  [ 8.20268765e-02  8.15424994e-02  8.10581222e-02 ...  8.04335326e-02
    8.04080293e-02  8.02805647e-02]
  [ 8.17719474e-02  8.12620744e-02  8.07267055e-02 ...  7.98726752e-02
    7.97961950e-02  7.95285106e-02]]

 [[ 2.22951382e-01  2.28725925e-01  2.25564584e-01 ...  3.11686873e-01
    3.13751996e-01  3.15613180e-01]
  [ 2.23091587e-01  2.24583030e-01  2.21459940e-01 ...  3.10310125e-01
    3.12591970e-01  3.14746350e-01]
  [ 2.19955742e-01  2.20006734e-01  2.19649807e-01 ...  3.09405059e-01
    3.11648637e-01  3.13956022e-01]
...
  [ 1.53792210e+01  1.53791199e+01  1.53784533e+01 ...  1.54257402e+01
    1.54255867e+01  1.54249716e+01]
  [ 1.53785553e+01  1.53781958e+01  1.53773260e+01 ...  1.54258442e+01
    1.54257402e+01  1.54252796e+01]
  [ 1.53781958e+01  1.53775826e+01  1.53765059e+01 ...  1.54261503e+01
    1.54261503e+01  1.54258947e+01]]

 [[ 1.68463440e+01  1.68462410e+01  1.68456783e+01 ...  1.67774582e+01
    1.67737160e+01  1.67702827e+01]
  [ 1.68426018e+01  1.68426018e+01  1.68419361e+01 ...  1.67784328e+01
    1.67744846e+01  1.67705383e+01]
  [ 1.68390675e+01  1.68389664e+01  1.68383999e+01 ...  1.67790985e+01
    1.67752018e+01  1.67710495e+01]
  ...
  [ 1.66797657e+01  1.66799202e+01  1.66800251e+01 ...  1.67107258e+01
    1.67112370e+01  1.67112370e+01]
  [ 1.66804848e+01  1.66806908e+01  1.66808949e+01 ...  1.67096996e+01
    1.67103653e+01  1.67105713e+01]
  [ 1.66809464e+01  1.66812019e+01  1.66814079e+01 ...  1.67086754e+01
    1.67094440e+01  1.67096996e+01]]], 'kilometer')>
Coordinates:
    number          int64 8B ...
    valid_time      datetime64[ns] 8B 2024-07-23T06:00:00
  * pressure_level  (pressure_level) float64 216B 1e+03 975.0 ... 125.0 100.0
  * latitude        (latitude) float64 680B 24.0 23.75 23.5 ... 3.5 3.25 3.0
  * longitude       (longitude) float64 808B 110.0 110.2 110.5 ... 134.8 135.0
    expver          <U4 16B ...

And for WVF_u:

<xarray.DataArray (pressure_level: 27, latitude: 85, longitude: 101)> Size: 927kB
<Quantity([[[-3.93375121e-02 -4.25711088e-02 -3.57796550e-02 ... -8.78068060e-02
   -8.19757804e-02 -7.58858547e-02]
  [-4.63820845e-02 -4.17888388e-02 -3.43502276e-02 ... -9.67195183e-02
   -9.02168304e-02 -8.35420415e-02]
  [-4.68467176e-02 -3.87023538e-02 -3.25297602e-02 ... -1.04566343e-01
   -9.89820808e-02 -9.43385214e-02]
  ...
  [ 7.97641203e-02  7.84882456e-02  7.21117482e-02 ...  2.79520098e-02
    2.69949734e-02  2.68616080e-02]
  [ 7.34354407e-02  7.12440759e-02  6.47160709e-02 ...  2.07603071e-02
    1.94839053e-02  1.92326922e-02]
  [ 6.57767281e-02  6.37183115e-02  5.74287325e-02 ...  1.13987373e-02
    9.09739826e-03  8.59958306e-03]]

 [[-3.75317186e-02 -4.18534800e-02 -3.64151001e-02 ... -8.95251855e-02
   -8.27948377e-02 -7.80771002e-02]
  [-4.64028344e-02 -4.54809740e-02 -3.90531942e-02 ... -9.86189842e-02
   -9.14381742e-02 -8.48975256e-02]
  [-5.65981083e-02 -5.27375229e-02 -4.21363339e-02 ... -1.05799690e-01
   -9.99277532e-02 -9.55018327e-02]
...
  [-1.12326372e-04 -1.09806926e-04 -1.07791471e-04 ... -4.47919301e-05
   -4.39612559e-05 -4.44887992e-05]
  [-1.07405052e-04 -1.04622130e-04 -1.02835133e-04 ... -4.15492104e-05
   -4.12685513e-05 -4.10974972e-05]
  [-1.15246665e-04 -1.09728840e-04 -1.04672712e-04 ... -3.92314178e-05
   -3.89648667e-05 -3.83068909e-05]]

 [[-1.03202095e-04 -1.04577986e-04 -1.07462663e-04 ... -2.44413532e-05
   -2.77988711e-05 -3.07358896e-05]
  [-1.14893250e-04 -1.17476266e-04 -1.21443351e-04 ... -1.94820186e-05
   -2.34835534e-05 -2.73336445e-05]
  [-1.31988840e-04 -1.33388428e-04 -1.34841292e-04 ... -1.36382887e-05
   -1.82957356e-05 -2.30648238e-05]
  ...
  [-3.23990025e-05 -3.31889169e-05 -3.28399801e-05 ... -4.17629008e-06
   -3.46313391e-06 -2.46502850e-06]
  [-3.24851280e-05 -3.33583776e-05 -3.30968905e-05 ... -4.55718919e-06
   -3.85740532e-06 -2.91683318e-06]
  [-3.24450157e-05 -3.27122798e-05 -3.19729370e-05 ... -5.28054079e-06
   -4.79992195e-06 -4.10016901e-06]]], 'meter / second')>
Coordinates:
    number          int64 8B 0
    valid_time      datetime64[ns] 8B 2024-07-23T06:00:00
  * pressure_level  (pressure_level) float64 216B 1e+03 975.0 ... 125.0 100.0
  * latitude        (latitude) float64 680B 24.0 23.75 23.5 ... 3.5 3.25 3.0
  * longitude       (longitude) float64 808B 110.0 110.2 110.5 ... 134.8 135.0
    expver          <U4 16B '0001'

Even a simple script such as the following would result in the same IndexError: Unlabeled multi-dimensional array cannot be used for indexing: pressure_level:

import numpy as np
import xarray as xr
import metpy.calc as mpcalc
import metpy.interpolate as mpinterp
from metpy.units import units

dataset = xr.open_dataset('C:/Users/Brian/Downloads/wvf_isobaric.nc')

x = (1.0 * (10 ** -6)) * units('(K*(m^2))/(kg*s)')

interp = mpinterp.interpolate_to_isosurface(dataset['pv'].isel(valid_time=0) * units('(K*(m^2))/(kg*s)'),
                                            dataset['t'].isel(valid_time=0) * units.kelvin,
                                            x,
                                            bottom_up_search = False)

Is there something wrong with the NetCDF file, or with how I manipulated the dataset? Or I am using the interpolate_to_isosurface function improperly? Or is there a better way to do this?