Issue adding height axis from pressure with Pa data

[rigoudyg]

Description The function add_height_from_pressure_axis works well if pressure data is in hPa. I work with pressure data in Pa (as required for CMIP6 data). In this case, the height axis is shifted.

To Reproduce

1. Input data can be downloaded on ESGF (whatever the node, for example https://esgf-node.llnl.gov/search/cmip6/), for instance searching for file: CMIP6.CMIP.CNRM-CERFACS.CNRM-CM6-1.historical.r1i1p1f2.Amon.ta.gr
2. Here is a code example using the previous file

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from matplotlib.ticker import ScalarFormatter
import geocat.viz as gv

# Open a netCDF data file using xarray default engine and load the data into xarrays
ds = xr.open_dataset("CMIP/CNRM-CERFACS/CNRM-CM6-1/historical/r1i1p1f2/Amon/ta/gr/latest/ta_Amon_CNRM-CM6-1_historical_r1i1p1f2_gr_185001-201412.nc", decode_cf=True, decode_times=True, decode_coords=True, use_cftime=True)

# Extract variables
T = ds.ta.isel(time=0).sel(lon=0)

# Generate figure (set its size (width, height) in inches) and axes
fig = plt.figure(figsize=(20, 20))
ax = plt.gca()

# Set y-axis to have log-scale
plt.yscale('log')

# Contour-plot T-data
p = T.plot.contour(ax=ax, levels=27, colors='red', extend='neither')
ax.clabel(p, fmt='%d', inline=1, fontsize=14, colors='k')

# Use geocat-viz utility function to add minor ticks to x-axis
gv.add_major_minor_ticks(ax, x_minor_per_major=3, y_minor_per_major=0)

# Use geocat.viz.util convenience function to set axes tick values
# Set y-lim inorder for y-axis to have descending values
gv.set_axes_limits_and_ticks(ax,
                             xticks=np.linspace(-60, 60, 5),
                             xticklabels=['60S', '30S', '0', '30N', '60N'],
                             yticks=T["plev"])
ax.invert_yaxis()

# Change formatter or else we tick values formatted in exponential form
ax.yaxis.set_major_formatter(ScalarFormatter())

# Create second y-axis to show geo-potential height.
axRHS = gv.add_height_from_pressure_axis(ax, pressure_units=T["plev"].units)

gv.set_titles_and_labels(ax=ax, maintitle="Monthly temperature - historical - lon 0 - time 185001")

# Show the plot
plt.show()

4. Look at the plot, the height axis is shifted.

Expected behavior The height axis should begin at 100000Pa (i.e. 1000 hPa) and not at 1000Pa.

Screenshots

Fix I think that the pressure list must be converted back to input pressure units before being set to ticks.

OS: Linux - ubuntu

Environment Install via pip

Package            Version
------------------ ---------
Cartopy            0.22.0
certifi            2023.7.22
cftime             1.6.2
charset-normalizer 3.2.0
click              8.1.7
cloudpickle        2.2.1
cmaps              2.0.1
contourpy          1.1.0
cycler             0.11.0
dask               2023.8.1
fonttools          4.42.1
fsspec             2023.6.0
geocat.viz         2023.7.0
h5netcdf           1.2.0
h5py               3.9.0
idna               3.4
importlib-metadata 6.8.0
joblib             1.3.2
kiwisolver         1.4.5
locket             1.0.0
matplotlib         3.7.2
MetPy              1.5.1
netCDF4            1.6.4
numpy              1.25.2
packaging          23.1
pandas             2.1.0
partd              1.4.0
Pillow             10.0.0
Pint               0.22
pip                22.3.1
platformdirs       3.10.0
pooch              1.7.0
pyparsing          3.0.9
pyproj             3.6.0
pyshp              2.3.1
python-dateutil    2.8.2
pytz               2023.3
PyYAML             6.0.1
requests           2.31.0
scikit-learn       1.3.0
scipy              1.11.2
setuptools         65.5.1
shapely            2.0.1
six                1.16.0
threadpoolctl      3.2.0
toolz              0.12.0
traitlets          5.9.0
typing_extensions  4.7.1
tzdata             2023.3
urllib3            2.0.4
wheel              0.38.4
xarray             2023.8.0
zipp               3.16.2

[kafitzgerald]

Thanks for reporting this! I'll take a look.

[kafitzgerald]

Thanks again for pointing this out. You were correct about the unit conversion. I'm guessing the sample datasets folks were working with all had pressure in hPa.

I put in a quick bug fix over in PR #152. Hopefully we can get this in quickly.

[rigoudyg]

Thank you for your answer, once available, I'll install the next version for working on plot tools. Best regards, Gaëlle

[kafitzgerald]

The fix should now be in the new release on PyPI and conda-forge.

Please go ahead and reopen this issue if the update doesn't fix your problem.

Thanks again for reporting this!