Plotting height anomaly as a cross section

[rolypolyman]

I am wondering if anyone knows a good way to show height anomaly on a vertical cross section chart (or even a horizontal one). I am trying to generate a few educational graphics showing this.

GFS supposedly has "Geopotential_height_anomaly_isobaric" in the variables but it's missing from the ncss.variables list so I decided to just ballpark it against a standard atmosphere using the function mpcalc.pressure_to_height_std(). I came up with the code below. I am not sure if I wrote it properly, and it is kind of messy since there was some weirdness with the units that I struggled with, but it does seem to generate marginally useful output.

Is there a better way to obtain height anomaly, maybe via subtraction of fields from normals via a Siphon query? I see a lot of websites doing this but it's not clear to me how.

Code for plotting on a vertical cross section:

        # pressure_levels_hPa = cross['isobaric']
        hgt = cross['Geopotential_height_isobaric']

        # Convert hPa to Pa (1 hPa = 100 Pa) and attach units
        pressure_levels_Pa = (cross['isobaric'].values) * units.Pa
        
        # Create a new DataArray with units
        pressure_levels_Pa = xr.DataArray(pressure_levels_Pa, coords=cross['isobaric'].coords, dims=cross['isobaric'].dims)
        
        # Quantify the DataArray with MetPy
        pressure_levels_Pa = pressure_levels_Pa.metpy.quantify()
        
        # Calculate standard atmosphere geopotential heights in kilometers
        std_height_km = mpcalc.pressure_to_height_std(pressure_levels_Pa)
        
        # Convert to meters and extract values
        std_height_m_values = std_height_km.values * 1000  # Convert to meters
        
        # Replicate std_height_m_values across the second dimension to match the shape of hgt.values
        num_points = hgt.shape[1]  # Assuming hgt has shape (41, 100)
        std_height_m_expanded = np.tile(std_height_m_values[:, np.newaxis], (1, num_points))
        
        # Perform subtraction using expanded values
        departure_values = hgt.values - std_height_m_expanded
        
        # Convert the result back into a DataArray
        departure = xr.DataArray(departure_values, coords=hgt.coords, dims=hgt.dims)

        hgtanom_contour = ax.contour(distances, cross['isobaric'], departure,
                                 levels=np.arange(-10000, 10000, 40), colors='g', linewidths=1, zorder=2)

        halabels = plt.clabel(hgtanom_contour, hgtanom_contour.levels[1::2], fontsize=8, colors='g', inline=1,
                   inline_spacing=8, fmt='%i', rightside_up=True, use_clabeltext=True, zorder=2)

[dopplershift]

First things, here's how you can accomplish calculating the difference without so much fighting of things:

# This all is assuming we already have a 2D cross section of geopotential height

# This gets the internal data storage for this DataArray to use a Pint Quantity so unit-tracking
# and math "just works"
hgt = cross['Geopotential_height_isobaric'].metpy.quantify()

# Due to limitations in xarray, since "isobaric" is used as coordinates for the dataset, we can't use the "quantify()"
# utility from above, so use MetPy to get "isobaric" as a plain Pint Quantity (a numpy array with units attached)
pressure_levels_Pa = cross['isobaric'].metpy.unit_array

# Calculate standard atmosphere heights from pressure.
std_height = mpcalc.pressure_to_height_std(pressure_levels_Pa)

# Perform the subtraction, using NumPy broadcasting:
# https://numpy.org/doc/stable/user/basics.broadcasting.html
# Essentially this changes the size of std_height to be (41, 1), which numpy can "broadcast" against
# hgt which is (41, 100).
# Any needed unit conversion is handled automatically by pint and the result has unit information
departure = hgt - std_height[:, None]

hgtanom_contour = plt.contour(distances, cross['isobaric'], departure,
                         levels=np.arange(-10000, 10000, 40), colors='g', linewidths=1, zorder=2)

halabels = plt.clabel(hgtanom_contour, hgtanom_contour.levels[1::2], fontsize=8, colors='g', inline=1,
           inline_spacing=8, fmt='%i', rightside_up=True, use_clabeltext=True, zorder=2)

This code works with the features of pint/xarray rather than fighting it or needing a lot of manual switching.

WIth that said, calculating anomalies from a standard atmosphere isn't something I've ever seen used. Normally, anomalies are calculated against a recent (30-year) climatological norm. If you can get that, the calculation should be even simpler.