使用 python matplotlib 和 metpy 添加辅助 y 轴

Question

我知道这个问题似乎与这里的许多其他问题类似，但我已经尝试过，不幸的是它们都没有解决我在尝试添加辅助 y 轴时当前面临的问题。

问题非常简单，但我找不到任何可以解决该问题的方法：在 SkewT 图上添加辅助 y 轴会更改图的 y 限制，而不仅仅是添加轴。

基本上，我希望添加辅助 y 轴，因为高度是使用 SkewT 内的压力水平来表示的，但也应该可以以公里为单位显示该高度。我想告诉第二个 y 轴：

1. 它的刻度应在 1015 到 100hPa 之间（就像原始 y 轴一样）；
2. 我只想在辅助 y 轴上显示 0、1、3、6、9、12、15 公里（简单的气压 (hPa) 到高度 (km) 转换）；
3. 我希望 0 公里从第一个压力水平开始并从那里开始；
4. 辅助 y 轴也应在 Y 中使用对数缩放。

这是辅助轴的示例，您可以看到与第一个轴相比缩放比例关闭：

这是我为获得更改而添加的代码，尽管仅添加第一行就会更改图表：

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())

这是一个更简单的示例，因此您可以使用此处找到的 Metpy 的简单发声代码示例自行测试：

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")

这可能只是一个简单的错误，或者 Metpy 本身可能存在一些问题，导致它变得如此

twinx()

，而类似的人并没有做我希望他们做的事情。我正在尝试找到一种解决方案，让我拥有第二个 y 轴，其压力值和缩放比例与第一个轴完全相同，然后我可以仅显示某些刻度并将这些刻度标签替换为相应的公里等效值。

谢谢！

Answer 1

现在使用 Matplotlib 的“实验性”（自 3.1 起）API 添加辅助轴可以最轻松地完成此操作。它使用MetPy的标准大气计算函数来在高度和压力之间进行转换：

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)')

此示例代码给出了以下图像（我将其剥离以保持代码完整但简单）：

由于这是假设标准大气，因此 0 公里刻度正好对应于 1013.25 mb 水平。

Answer 2

我看到 @dopplershift 的响应仅在下 ylim 设置为 1013.25 hPa 或更低时才有效。

否则，辅助轴似乎保留父轴的对数刻度，因此它不能具有负值。

如果我们设置：

skew.ax.set_ylim(1100, 100)

我们得到了不一致的次轴。请注意，1100 hPa 的值设置为 0 米。

图链接

Answer 3

# Get the surface pressure
surface_pressure = p[0].values

# Function to convert height to pressure using the log scale
def height_to_pressure_at_surface(heights, units, surface_pressure):
    pressures = [mpcalc.height_to_pressure_std(h * units).to('hPa').magnitude + surface_pressure - 1013.25 for h in heights]
    return pressures

# Heights in meters and feet
height_ticks_m = [0, 700, 1500, 3000, 5000, 7000, 9000, 12000, 15000]
height_labels_m = ['0m', '700m', '1.5km', '3km', '5km', '7km', '9km', '12km', '15km']
height_ticks_ft = [0, 2500, 5000, 10000, 18000, 30000, 34000, 39000, 50000]
height_labels_ft = ['0ft', '2,500ft', '5,000ft', '10,000ft', '18,000ft', '30,000ft', '34,000ft', '39,000ft', '50,000ft']

# Convert heights to corresponding pressures
pressure_ticks_m = height_to_pressure_at_surface(height_ticks_m, units.m, surface_pressure)
pressure_ticks_ft = height_to_pressure_at_surface(height_ticks_ft, units.ft, surface_pressure)

# Add a secondary axis for heights in meters
secax = skew.ax.secondary_yaxis('left')
secax.spines['left'].set_position(('outward', 60))  # Adjust position to prevent overlap
secax.set_yticks(pressure_ticks_m)
secax.set_yticklabels(height_labels_m)


# Add a third axis for heights in feet
thirdax = skew.ax.secondary_yaxis('left')
thirdax.spines['left'].set_position(('outward', 100))  # Adjust position to prevent overlap
thirdax.set_yticks(pressure_ticks_ft)
thirdax.set_yticklabels(height_labels_ft)

# Add a horizontal line at surface pressure
skew.ax.axhline(surface_pressure, color='black', linestyle='--', linewidth=1)

# Change to adjust data limits and give it a semblance of what we want
skew.ax.set_ylim(1070, 100)
skew.ax.set_xlim(-50, 50)