Adding secondary y axis using python matplotlib with metpy

I know this question seems similar to a lot of other questions here, but I've tried them and unfortunately none of them are addressing the problem I'm currently facing when trying to add a secondary y-axis.

The problem is quite simple, but I can't find anything that would fix it: Adding a secondary y-axis on my SkewT plot changes the y-limits of the plot instead of just adding the axis.

Basically, I wish to add a secondary y-axis since the height is presented using pressure levels inside a SkewT, but it should also be possible to show that height in km as well. I want to tell the second y-axis that:

1. It should scale between 1015 and 100hPa (just like the original y-axis);
2. I want to show only 0, 1, 3, 6, 9, 12, 15 km on the secondary y-axis (simple pressure (hPa) to height (km) conversion);
3. I want the 0km to start on the first pressure level and scale from there;
4. The secondary y-axis should also use log scaling in Y.

Here's an example of what it looks like with the secondary axis, you can see that the scaling if off compared to the first axis:

Here's the bit of code I added to get the change, although just adding the first line changes the graph:

twin = skew.ax.twinx()
twin.set_yscale('log')
twin.spines['right'].set_position(('axes', 0))
twin.set_frame_on(True)
twin.patch.set_visible(False)
twin.set_ylim(skew.ax.get_ylim())

Here is a simpler example so you can test it out yourselves using Metpy's Simple Sounding code example found here:

import matplotlib.pyplot as plt
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

plt.rcParams['figure.figsize'] = (9, 9)

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)

# Drop any rows with all NaN values for T, Td, winds
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)

skew = SkewT()

# Plot the data using normal plotting functions, in this case using
# log scaling in Y, as dictated by the typical meteorological plot
skew.plot(p, T, 'r')
skew.plot(p, Td, 'g')
skew.plot_barbs(p, u, v)

# Add the relevant special lines
skew.plot_dry_adiabats()
skew.plot_moist_adiabats()
skew.plot_mixing_lines()
skew.ax.set_ylim(1000, 100)
# twin = skew.ax.twinx()
# twin.set_yscale('log')
# twin.spines['right'].set_position(('axes', 0))
# twin.set_frame_on(True)
# twin.patch.set_visible(False)
# twin.set_ylim(skew.ax.get_ylim())

plt.savefig("metpy_base.png")

It could just be a simple mistake or there could be something with Metpy itself that makes it so twinx() and the the likes aren't doing what I want them to do. I'm trying to find a solution that just lets me have a second y-axis with the exact same pressure values an scaling as the first one, where I can then show only certain ticks and replace those tick labels with their appropriate km equivalent.

Thanks!

This is most easily accomplished now using Matplotlib's "experimental" (since 3.1) API for adding a secondary axis. It uses MetPy's standard atmosphere calculation functions to transform between height and pressure:

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.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, 'hPa')).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)')

This sample code gives the following image (I stripped it down to keep the code complete but straightforward):

Since this is assuming a standard atmosphere, the 0km tick corresponds exactly to the 1013.25 mb level.