How to properly produce a cross section with model levels and vertical interpolation

[guidocioni]

I want to produce a cross section using data on model levels to properly capture the structure of a cold front. Note that I do have data on pressure levels but I want to use model hybrid levels to have a better vertical resolution.

The input dataset (obtained from the operational COSMO-D2 model, https://opendata.dwd.de/weather/nwp/cosmo-d2/) has some standard variables defined on the same vertical axis and a 2D variable containing the height of model levels in meters above sea level. What I want to obtain is a cross section on height levels, so not on model levels. I figured out that I have to do vertical interpolation as well as the usual cross section interpolation, so here is what I've come up with.

First compute the cross section interpolation for all the variables in the dataset.

start = (48, 5.5)
end = (55.5, 13.5)

cross = cross_section(dset, start, end, steps=300, interp_type='nearest')

This data is still on model levels so I cannot plot it on a contour plot with height as y axis. In oder to interpolate vertically from model levels to height levels the only function that (I believe) exists is interpolate_1d which unfortunately does not allow to pass an entire dataset but only arrays without the time dimension.
To decide which levels are to be used for the interpolation I decided to take a sample of levels where the terrain is "flat". This way the data "inside" the orography will be hopefully masked.
If HHL is the variable that contains the 2D height of model levels

vertlevels = cross.HHL[:, cross.HHL.argmin(dim='index').min().compute().item()].values.round()

returns

array([8020., 7542., 7085., 6646., 6227., 5827., 5445., 5081., 4734.,
       4404., 4091., 3794., 3512., 3246., 2994., 2756., 2532., 2322.,
       2124., 1939., 1765., 1604., 1453., 1313., 1183., 1063.,  953.,
        851.,  757.,  671.,  593.,  522.,  457.,  399.,  346.,  298.,
        256.,  217.,  182.,  151.,  123.,   97.,   73.,   51.,   30.,
         10.], dtype=float32)

Makes sense. Let's use this to do vertical interpolation of some variables

timestep = 2
data = cross.isel(time=timestep)
theta_interp, t_wind_interp, n_wind_interp = interpolate_1d(vertlevels,
                                                            data.HHL.values,
                                                            data.theta_e.values,
                                                            data.t_wind.values,
                                                            data.n_wind.values)

Let's plot these values by setting a background so that we can see where the array is masked

fig = plt.figure(1, figsize=(16., 9.))
ax = plt.axes()
ax.set_facecolor('brown')

cmap = get_colormap('temp_meteociel')

theta_contour = ax.contourf(cross['lon'],
                           vertlevels,
                           theta_interp,
                           levels=np.arange(-10, 40, 0.5),
                           cmap=cmap,
                           extend='both')

After some polishing this is what I obtain

I think the results are correct: I can see what I would expect and the data is masked over mountains.

But the question is...isn't there a better way to do this? Especially for the vertical interpolation method I believe there should be a function that just takes the whole cross section dataset (cross in this case) to do a vertical interpolation over defined levels. I believe this is what most of the people would like to see if they're using model levels. Right now interp_1d only accepts numpy arrays.

I hope I didn't overlooked something...