It only takes one line to interpolate to cross-section using metpy.

It's just the cleanup before (parse_cf) and after (getting coordinates to show on x-axis) that might be tricky!

#python #atmos #datavi

from io import BytesIO
import metpy.interpolate
import requests
import xarray as xr
import matplotlib.pyplot as plt

# download dataset to memory and read it with xarray
url = "https://downloads.psl.noaa.gov/Datasets/ncep.reanalysis/Monthlies/pressure/air.mon.1981-2010.ltm.nc"
with requests.get(url) as response:
    with BytesIO(response.content) as fobj:
        ds = xr.open_dataset(fobj, decode_times=False)[["air"]].isel(time=0).load()

# clean up bad units and parse cf metadata
ds["time"].attrs.pop("units")
ds = ds.metpy.parse_cf()

# create the cross section (lat, lon)
start = (40, 235)
end = (80, 260)
ds_cross = metpy.interpolate.cross_section(ds, start, end)
ds_cross["coords"] = "index", [
    f"({lat:.1f}N, {lon:.1f}E)"
    for lat, lon in zip(ds_cross["lat"].values, ds_cross["lon"].values)
]

# plot the cross section
fig, ax = plt.subplots(1, 1, figsize=(14, 6))
ax.contourf(ds_cross["coords"], ds_cross["level"], ds_cross["air"], cmap="coolwarm")
ax.set_ylim(ax.get_ylim()[::-1])
ax.set_xlabel("Coordinates")
ax.set_ylabel("Pressure (mb)")
ax.set_xticks(ax.get_xticks()[::16])
plt.show()