Height scale for SkewT

[dopplershift]

Should be a standalone graphical element, secondary scale that pulls from the standard atmosphere conversion function we have. Should be do-able modelling off of some of matplotlib's examples with multiple scales.

[guytcc]

Borrowing liberally from SkewT python 3 port, I added a height scale to a plot with the code below (the class is just a holder). Most would be easy to port into Metpy.SkewT with some finagling.

Two obvious problems would need to be worked out:

1. The axis for height doesn't show up. I've tried a number of settings with no luck.
2. Better(smart?) positioning of the axes and wind barbs.

from metpy.units import units

class Sondes(object):
    def __init__(self, sounding, dt):
        self.file = sounding
        self.dt = dt
        self.dt_str = dt.strftime('%Y-%m-%d %HZ')

        self.p = sounding.variables['pressure'][:]
        self.T = sounding.variables['temperature'][:]
        self.Td = sounding.variables['dewpoint'][:]
        self.u = sounding.variables['u_wind'][:]
        self.v = sounding.variables['v_wind'][:]

        # Calculate the parcel profile.
        self.H = mpcalc.pressure_to_height_std(self.p)
        self.zero_levs = mpcalc.find_intersections(self.H, self.T, np.zeros_like(self.T))

    def plot(self, fig=None, xlims=None, ylims=None, plot_barbs=False, bold_temp=None, add_height=False):
        if fig is None:
            fig = plt.figure(figsize=(7, 7))
        skew = SkewT(fig)
        _ = skew.ax.set_title(self.dt_str)
        _ = skew.ax.set_ylim(ylims)

        skew.plot(self.p, self.T, 'r', linewidth=2)
        skew.plot(self.p, self.Td, 'g', linewidth=2)

        if bold_temp is not None:
            for temp in bold_temp:
                skew.ax.axvline(temp, color='c', linestyle='--', linewidth=2)

        if plot_barbs:
            _ = skew.plot_barbs(self.p, self.u, self.v)

        # This order of plotting matters
        _ = skew.ax.set_xlim(xlims)
            
        if add_height:
            majorLocatorKM = MultipleLocator(1)
            majorLocatorKFT = MultipleLocator(5)
            minorLocator = MultipleLocator(1)

            plims = skew.ax.get_ylim() * units('hectopascal')
            zmin = mpcalc.pressure_to_height_std(plims[0]).magnitude
            zmax = mpcalc.pressure_to_height_std(plims[1]).magnitude
            zminf = zmin * 3.2808
            zmaxf = zmax * 3.2808

            kmhax = fig.add_axes([0.775, 0.1, 1e-6, 0.8], frameon=True)
            kmhax.xaxis.set_ticks([], [])
            kmhax.spines['left'].set_color('k')
            kmhax.spines['right'].set_visible(False)
            kmhax.tick_params(axis='y', colors='k', labelsize=8)
            kmhax.set_ylim(zmin, zmax)
            kmhax.set_title("km/kft", fontsize=10)
            kmhax.get_yaxis().set_tick_params(which="both", direction='out')
            kmhax.yaxis.set_major_locator(majorLocatorKM)
            kmhax.yaxis.set_minor_locator(minorLocator)

            fthax = kmhax.twinx()
            fthax.xaxis.set_ticks([], [])
            fthax.tick_params(axis='y', colors='k', labelsize=8)
            fthax.set_ylim(zminf, zmaxf)
            fthax.get_yaxis().set_tick_params(which="both", direction='out')
            fthax.yaxis.set_major_locator(majorLocatorKFT)
            fthax.yaxis.set_minor_locator(minorLocator)
        return skew

[dopplershift]

Thanks for digging that up! A few comments:

1. I'd really like to have it separated from the plotting call
2. In an ideal world, it'd be dynamically updated from the y-axis in pressure coords. I'm not sure if a
   matplotlib spine is enough, or if maybe we do need the separate Axes, but we use matplotlib's event
   handling to update the range on the height scale
3. I don't think the scaling for heights is correct. Pressure is plotted with log scaling, but I don't see any
   log scaling associated with the heights.

[guytcc]

I may misunderstand, but I just want to note that the plot in the class above is not equivalent to the plot in Metpy.SkewT, but is a handler that basically does most things from the individual examples. I'm sure matplotlib event handling could be used for dynamic updates, but I'd have to look into that, since the values for the heights are computed from the limits of the plot maybe this would be easier. (?)

I'm open to another method. Just went the "easy" route of least resistance (development).

The scaline for the heights is a linear scale. I made no attempt to have a logarithmic axis as we often don't think of altitude the same way we do atmospheric pressure.

[dopplershift]

So what I'd want to merge is a patch that adds a method to the SkewT class that allows the user to turn on the height scale bar.

As far as the scaling is concerned, I think there's some confusion. You start by taking the pressure values and converting them to heights--no problem here. When we make the skewT, pressure is logarithmically spaced (not changing the values) along the y-axis. If you want the heights to properly align with the pressure values, the height axis has to be logarithmically scaled as well.

[dopplershift]

Nevermind, I mis-read the code before. I see what's going on and everything is aligned just fine.

[dopplershift]

Here was a quick hack we put together at AMS:

heights = np.array([1, 2, 3, 4, 5]) * units.km
std_pressures = mpcalc.height_to_pressure_std(heights)
for height_tick, p_tick in zip(heights, std_pressures):
    trans, _, _ = skew.ax.get_yaxis_text1_transform(0)
    skew.ax.text(0.02, p_tick, '---{:~d}'.format(height_tick), transform=trans)

[dopplershift]

Matplotlib has a "secondary axis" (provisional) functionality that might show a good path forward on this.

[dopplershift]

Looks like the secondary axis support in matplotlib is pretty useful. I just put this together to answer a SO question:

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

import metpy.calc as mpcalc
from metpy.cbook import get_test_data
from metpy.plots import SkewT
from metpy.units import units

col_names = ['pressure', 'height', 'temperature', 'dewpoint', 'direction', 'speed']
df = pd.read_fwf(get_test_data('jan20_sounding.txt', as_file_obj=False),
                 skiprows=5, usecols=[0, 1, 2, 3, 6, 7], names=col_names)
df = df.dropna(subset=('temperature', 'dewpoint', 'direction', 'speed'
                       ), how='all').reset_index(drop=True)

p = df['pressure'].values * units.hPa
T = df['temperature'].values * units.degC
Td = df['dewpoint'].values * units.degC
wind_speed = df['speed'].values * units.knots
wind_dir = df['direction'].values * units.degrees
u, v = mpcalc.wind_components(wind_speed, wind_dir)

# Standard Skew-T Plot
skew = SkewT()
skew.plot(p, T, 'r')
skew.plot(p, Td, 'g')
skew.plot_dry_adiabats()
skew.plot_moist_adiabats()
skew.plot_mixing_lines()
skew.ax.set_ylim(1015, 100)
skew.ax.set_ylabel('Pressure (hPa)')

# Add a secondary axis that automatically converts between pressure and height
# assuming a standard atmosphere. The value of -0.12 puts the secondary axis
# 0.12 normalized (0 to 1) coordinates left of the original axis.
secax = skew.ax.secondary_yaxis(-0.12,
    functions=(lambda p: mpcalc.pressure_to_height_std(units.Quantity(p, 'mbar')).m_as('km'),
               lambda h: mpcalc.height_to_pressure_std(units.Quantity(h, 'km')).m))
secax.yaxis.set_major_locator(plt.FixedLocator([0, 1, 3, 6, 9, 12, 15]))
secax.yaxis.set_minor_locator(plt.NullLocator())
secax.yaxis.set_major_formatter(plt.ScalarFormatter())
secax.set_ylabel('Height (km)')

I think we should bake that last block into the SkewT class as a method to add the height scale. There are a few API questions there, but I think that should be pretty easy to add.

[guidocioni]

Hey @dopplershift, this may be a stupid question but...what if I already have the height values in the original sounding?

Using standard atmosphere can cause some large discrepancies sometimes, but I cannot find a way to use an array for the values of the secondary y axis instead than a function.

[dopplershift]

If you have the data, then you could just go ahead and manually annotate the plot with that using matplotlib's text()/annotate() methods. You could also use the code above and do some kind of linear/nearest neighbor interpolation as the functions.

One simple approach is, given that you have pairs of height/pressure, you could use the minor ticks with fixed labels for these like this:

pres = [977, 909, 707, 575, 325, 190, 120]
labels = ['sfc', '1km', '3km', '6km', '9km', '12km', '15km']
skew.ax.tick_params(axis='y', which='minor', color='red', pad=20, length=20)
skew.ax.yaxis.set_ticks(pres, labels, minor=True, color='red')

This gives:

Note that with this approach, if you have a minor tick that duplicates a major tick, the minor tick won't show.

[guidocioni]

If you have the data, then you could just go ahead and manually annotate the plot with that using matplotlib's text()/annotate() methods. You could also use the code above and do some kind of linear/nearest neighbor interpolation as the functions.

One simple approach is, given that you have pairs of height/pressure, you could use the minor ticks with fixed labels for these like this:

pres = [977, 909, 707, 575, 325, 190, 120]
labels = ['sfc', '1km', '3km', '6km', '9km', '12km', '15km']
skew.ax.tick_params(axis='y', which='minor', color='red', pad=20, length=20)
skew.ax.yaxis.set_ticks(pres, labels, minor=True, color='red')

This gives:

Note that with this approach, if you have a minor tick that duplicates a major tick, the minor tick won't show.

Ah ok, I was hoping for a more clean way of doing it directly into metpy api without using direct text annotations :)

I will maintain the secondary_y_axis approach for now, thanks!

I tried to create a function that just returns the array of height no matter the input but that doesn't seem to work because I don't really know when and with which arguments those transforming functions are called. Also a function that returns just the nearest value does not work... The doc says the function should expect numpy arrays as input but to me it is not clear what is passed.