Help needed to calculate height for NCL_h_lat_6.py

[sjsprecious]

I am working on the visualization of the CESM output with the GeoCAT template and one type of plot is the vertical distribution (see example at https://geocat-examples.readthedocs.io/en/latest/gallery/Contours/NCL_h_lat_6.html#sphx-glr-gallery-contours-ncl-h-lat-6-py).

According to the notes in the example above, the height values on the right y-axis are not calculated based on the model output, but just arbitrary data. Is there any existing function I could use to calculate the height values based on my model output first and plot them correspondingly? Thanks.

[jukent]

I can look at this tomorrow!

[jukent]

Hi @sjsprecious I did not make this example so I might have to confer with some more senior GeoCAT team members but I am confused. You're referring to this part of the example (?):

# Create second y-axis to show geo-potential height.
# Currently we're using arbitrary values for height as we haven't figured out
# how to make this work properly yet
axRHS = ax.twinx()

# Use geocat.viz.util convenience function to set axes tick values
gv.set_axes_limits_and_ticks(axRHS, ylim=(0, 32), yticks=np.arange(4, 32, 4))
axRHS.tick_params(labelsize=12)  # manually set tick label size

So since everything is technically plotted against the left-hand y-axis, the right y-axis just shows the conversion between Pressure and height? It has no bearing on the plot, so ultimately we'd want to do the conversion math to create our desired yticks numpy array.

However, I'm wondering is this conversion constant across the different latitude values on our x-axis? I thought changes in temperature or water vapor would change the hight/pressure relationship. If it isn't a constant conversion, we can't accurately display both height and pressure as y-axis anyway - we'd have to pick one.

[jukent]

We will look into using the "International Standard Atmosphere" to do this conversion. A good option is to use MetPy: https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.pressure_to_height_std.html

[sjsprecious]

Thanks @jukent for looking into this issue. I agree that the height may not be a constant across the latitude values but an approximate estimation should be good enough for me. I will try the function you suggested and let you know how it works.

[jukent]

I am going to work on improving the example to use this function as well.

[jukent]

@sjsprecious I have made the change in https://github.com/NCAR/GeoCAT-examples/pull/500

[sjsprecious]

Hi @jukent , thanks for your help. I am able to follow your example to add the height values now.

[jukent]

Great @sjsprecious we're actually debating whether this answer is true. It assumes that heights are evenly spaced, but they aren't exactly. I will work on this more and ping you.

[jukent]

I think the new method is much better! @sjsprecious but I will check with the team next week.

[jukent]

@sjsprecious we've agreed on the solution now. Merging the PR auto-closed this issue, but I want to make sure you are satisfied with the steps and can follow along. Essentially here is the code for making the height axis properly aligned with the pressure axis:

# Create second y-axis to show geo-potential height.
axRHS = ax.twinx()

# Select heights to display as tick labels
heights_nice = np.arange(0,32,4)
# Send "nice" height values back to pressure  as tick locations
pressures_nice = mpcalc.height_to_pressure_std(heights_nice * units('km')).magnitude

axRHS.set_yscale('log')
axRHS.minorticks_off()

gv.set_axes_limits_and_ticks(axRHS,
                            ylim=ax.get_ylim(),
                            yticks=pressures_nice,
                            yticklabels=heights_nice)

So technically both the left and right y-axis are in Pressure, and the limits of the right are inherited from the left. We just picked pretty height values (yticklabels), and found their corresponding pressure values to indicate where they go (yticks).

There are ways to auto-detect what the range for your height values should be (by converting max and min pressure to height using the similar metpy function pressure_to_height_std).

Setting the yscale to logrithmic is vital, but the log has auto minor tick labels that are at natural pressure values so they look confusing with our selected height values (because they are subdividing a different unit).

I hope that makes sense, and please reach out if you have any more questions!

[sjsprecious]

Hi @jukent , thanks for your & teammate's efforts on this. The new method indeed looks better and I could confirm that I am able to follow your example code to make the desired plot. Thank you again for your efficient help!