Add a Thickness Calculation #621

jthielen · 2017-11-24T15:41:46Z

A pull request working on issue #295.

This adds a thickness calculation using an integral form of the hypsometric equation with virtual temperature adjustment. Also included are two tests of this calculation, one using matching data from other tests and one using a simple dry isothermal layer.

I don't yet have any text citations for the formula (I was just pulling from my notes), but once I get back to campus I should be able add that if needed.

Any feedback is greatly appreciated!

jrleeman

Thanks for another excellent PR! A couple little comments and the test error are about all that needs addressed.

jrleeman · 2017-11-27T16:10:53Z

metpy/calc/thermo.py

+    mixing = kwargs.pop('mixing', np.zeros_like(temperature))
+    molecular_weight_ratio = kwargs.pop('molecular_weight_ratio', epsilon)
+    bottom = kwargs.pop('bottom', None)
+    depth = kwargs.pop('depth', pressure.max() - pressure.min())


I think depth = kwargs.pop('depth', None) should work here as well.

Unfortunately, get_layer defaults to a depth of 100 hPa when None is specified, which didn't seem to make sense for this function (compared to the data of the whole layer). But, I can definitely still change it to that if you would like.

Good point. Do you really need to call get_layer if you don’t have a depth argument then?

That's a good point, I don't think I really need to in that case. So, I've made the get_layer call conditional now in my latest commit.

jrleeman · 2017-11-27T16:11:51Z

metpy/calc/thermo.py

+    layer_virttemp = virtual_temperature(layer_temp, layer_w, molecular_weight_ratio)
+
+    # Take the integral and return the result in meters
+    return (- Rd / g * np.trapz(layer_virttemp / layer_p, x=layer_p) *


dlnp right?

Thank you for pointing that out...I'm used to the continuous case where d(p)/p and d(ln(p)) are equivalent. But, after playing around with it some more, d(ln(p)) does indeed give the better results compared to what the layer-average form of the hyposometric equation would lead me to expect. I've corrected this in the commit I just made.

jrleeman · 2017-11-27T16:13:22Z

metpy/calc/tests/test_thermo.py

+    assert_almost_equal(thickness, 10725.92 * units.m, 2)
+
+
+def test_thickness_hydrostatic_isothermal():


Travis error is coming from rounding on legacy here:

E AssertionError: E Arrays are not almost equal to 2 decimals E ACTUAL: 5542.1180018995983 E DESIRED: 5542.11

Hopefully this should be fixed with the formula and test updates I was just able to make.

dopplershift

This is an excellent PR! I only have some concerns about whether the handling for optional mixing ratio is a good idea for potentially large datasets.

dopplershift · 2017-11-28T16:11:56Z

metpy/calc/thermo.py

+    pressure_to_height_std, virtual_temperature
+
+    """
+    mixing = kwargs.pop('mixing', np.zeros_like(temperature))


Instead of calculating virtual temp with 0 mixing, could you instead do:

mixing = kwargs.pop(‘mixing’, None) if mixing is None: layer_virttemp = temperature else: layer_virttemp = virtual_temperature(...)

I appreciate the elegance of your original approach, but I feel like for a large enough grid (and deep enough profile), the array of zeros could have significant overhead for both memory and runtime.

Of course now that makes the get_layer call conditional as well, but I think that’s worth saving the extra calculation. You’d also need another test case for the case where you’re not passing in w. Thoughts, @jthielen ?

That sounds good to reduce the overhead of the zeros approach. I have made the calls conditional now in my most recent commit, and added tests to cover both situations. Let me know if that is what you were thinking.

dopplershift · 2017-11-28T16:19:57Z

metpy/calc/thermo.py

+    mixing = kwargs.pop('mixing', np.zeros_like(temperature))
+    molecular_weight_ratio = kwargs.pop('molecular_weight_ratio', epsilon)
+    bottom = kwargs.pop('bottom', None)
+    depth = kwargs.pop('depth', pressure.max() - pressure.min())


Good point. Do you really need to call get_layer if you don’t have a depth argument then?

dopplershift · 2017-11-28T16:21:50Z

metpy/calc/thermo.py

+    temperature : `pint.Quantity`
+        Atmospheric temperature profile
+    mixing : `pint.Quantity`, optional
+        Profile of dimensionless mass mixing ratio. Defaults to zero, so that virtual


Would need to update if the handling of no mixing ratio is changed.

dopplershift · 2017-11-28T16:24:25Z

@jrleeman I don’t recall: does get_layer support > 1D arrays?

jrleeman · 2017-11-28T17:51:34Z

Hmmm.. I don't remember - I'll have to try it, but I believe it does. (We do use it in some of the indices.)

jthielen · 2017-11-28T21:05:30Z

In case it's useful, I wasn't able to get get_layer to work with 2D arrays. Here is a quick example:

import numpy as np
from metpy.calc.tools import get_layer
from metpy.units import units

# Works
pressure = [850, 700, 500] * units.hPa
temperature = [5, 0, -2] * units.degC
print(get_layer(pressure, temperature, bottom=750*units.hPa))

# Fails
pressure = [[850, 700, 500],
            [850, 700, 500],
            [850, 700, 500]] * units.hPa
temperature = [[5, 0, -2],
               [4, 1, 0],
               [4, 0, -1]] * units.degC
print(get_layer(pressure, temperature, bottom=750*units.hPa))

Output:

[<Quantity([750 700 650], 'hectopascal')>, <Quantity([ 1.77673794  0.         -0.44049977], 'degC')>]

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-55-7723eafe4a3d> in <module>()
     15                [4, 1, 0],
     16                [4, 0, -1]] * units.degC
---> 17 print(get_layer(pressure, temperature, bottom=750*units.hPa))

~/miniconda3/envs/devel/lib/python3.6/site-packages/MetPy-0.6.1+12.geed0fca-py3.6.egg/metpy/units.py in wrapper(*args, **kwargs)
    206                                 '    x = x * units.meter / units.second')
    207                     raise ValueError(msg)
--> 208                 return func(*args, **kwargs)
    209 
    210             return wrapper

~/miniconda3/envs/devel/lib/python3.6/site-packages/MetPy-0.6.1+12.geed0fca-py3.6.egg/metpy/calc/tools.py in get_layer(pressure, *args, **kwargs)
    551     bottom_pressure, bottom_height = _get_bound_pressure_height(pressure, bottom,
    552                                                                 heights=heights,
--> 553                                                                 interpolate=interpolate)
    554 
    555     # Calculate the top if whatever units depth is in

~/miniconda3/envs/devel/lib/python3.6/site-packages/MetPy-0.6.1+12.geed0fca-py3.6.egg/metpy/calc/tools.py in _get_bound_pressure_height(pressure, bound, heights, interpolate)
    333     if bound.dimensionality == {'[length]': -1.0, '[mass]': 1.0, '[time]': -2.0}:
    334         # If the bound is in the pressure data, we know the pressure bound exactly
--> 335         if bound in pressure:
    336             bound_pressure = bound
    337             # If we have heights, we know the exact height value, otherwise return standard

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

This commit adds a thickness calculation using an integral form of the hypsometric equation with virtual temperature adjustment. This assumes hydrostatic balance. Also included are two tests of this calculation, one using matching data from other tests and one using a simple dry isothermal layer.

Based on feedback, corrects the thickness formula to use d(ln(pres)) instead of d(pres)/pres. While these are equivalent in the continuous case, they do not give the same results in the discrete case here. Tests have likewise been updated to account for this change in formula (which should hopefully also resolve the legacy test failure).

Based on review feedback, updates the optional mixing ratio to be handled conditionally rather than with an array of zeros (if none is given), as well as conditionally calling get_layer only when a bottom/depth is specified. Adds two additional tests to cover all four scenarios (full layer-moist, subset-moist, full layer-dry, subset-dry). Also adds a citation for the thickness formula.

The recent addition of the static energy calculations unfortunately caused a messy rebase conflict where the new functions got interleaved. I thought I got it all cleaned up, but unfortunately missed this exporter statement.

jrleeman · 2017-11-28T22:30:08Z

Yep - My remembering was incorrect. It isn't grid compatible ☹️

dopplershift

I think this is in good shape now. Thanks for making those changes. Good @jrleeman ?

jrleeman

LGTM

dopplershift · 2017-12-01T22:34:36Z

Thanks for the help! 🎉

jthielen requested review from dopplershift and jrleeman as code owners November 24, 2017 15:41

jrleeman suggested changes Nov 27, 2017

View reviewed changes

dopplershift added Area: Calc Pertains to calculations Type: Feature New functionality labels Nov 28, 2017

dopplershift added this to the 0.7 milestone Nov 28, 2017

dopplershift reviewed Nov 28, 2017

View reviewed changes

jthielen added 4 commits November 28, 2017 15:27

Correction After Rebase with Upstream

c5fdfd8

The recent addition of the static energy calculations unfortunately caused a messy rebase conflict where the new functions got interleaved. I thought I got it all cleaned up, but unfortunately missed this exporter statement.

dopplershift approved these changes Nov 29, 2017

View reviewed changes

jrleeman approved these changes Nov 29, 2017

View reviewed changes

Merge branch 'master' into thickness-hydrostatic

6b70dc7

dopplershift merged commit ab550b8 into Unidata:master Dec 1, 2017

jthielen deleted the thickness-hydrostatic branch December 2, 2017 05:43

jrleeman mentioned this pull request Dec 6, 2017

Height from pressure, temp, relh #295

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add a Thickness Calculation #621

Add a Thickness Calculation #621

jthielen commented Nov 24, 2017

jrleeman left a comment

jrleeman Nov 27, 2017

jthielen Nov 27, 2017

dopplershift Nov 28, 2017

jthielen Nov 28, 2017

jrleeman Nov 27, 2017

jthielen Nov 27, 2017

jrleeman Nov 27, 2017

jthielen Nov 27, 2017

dopplershift left a comment

dopplershift Nov 28, 2017

dopplershift Nov 28, 2017

jthielen Nov 28, 2017

dopplershift Nov 28, 2017

dopplershift Nov 28, 2017

dopplershift commented Nov 28, 2017

jrleeman commented Nov 28, 2017

jthielen commented Nov 28, 2017

jrleeman commented Nov 28, 2017

dopplershift left a comment

jrleeman left a comment

dopplershift commented Dec 1, 2017

		assert_almost_equal(thickness, 10725.92 * units.m, 2)


		def test_thickness_hydrostatic_isothermal():

Add a Thickness Calculation #621

Add a Thickness Calculation #621

Conversation

jthielen commented Nov 24, 2017

jrleeman left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dopplershift left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dopplershift commented Nov 28, 2017

jrleeman commented Nov 28, 2017

jthielen commented Nov 28, 2017

jrleeman commented Nov 28, 2017

dopplershift left a comment

Choose a reason for hiding this comment

jrleeman left a comment

Choose a reason for hiding this comment

dopplershift commented Dec 1, 2017