[Python] MetPy "interpolate_to_isosurface" results in "IndexError: Unlabeled multi-dimensional...

I need help on the usage of MetPy's interpolate_to_isosurface function (link).

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":

print(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?

And is there a better way to do this, or a better Python library for this purpose?

Continue reading...