Repository: Unidata/MetPy

Issue:
Add non-weighted pressure-based mean function
While we have `mean_pressure_weighted`, we should add a corresponding integral-based mean calculation to go along with it. For instance, something like `bunkers_storm_motion` should probably be using that instead of `numpy.mean`, since the latter would be biased by uneven sampling of levels in a vertical profile. It should be as simple as taking `mean_pressure_weighted` and ripping out the weighting (and adjusting the denominator to `p[0] - p[-1]`).

Not sure what to call it...`mean_vertical`? `vertical_mean`? `mean_non_weighted`?

Discussion:

Patch:
diff --git a/src/metpy/calc/indices.py b/src/metpy/calc/indices.py
index 4678740a41..ef258242e8 100644
--- a/src/metpy/calc/indices.py
+++ b/src/metpy/calc/indices.py
@@ -284,9 +284,9 @@ def bunkers_storm_motion(pressure, u, v, height):
     >>> sped = [5, 15, 20, 30, 50, 60] * units.knots
     >>> u, v = wind_components(sped, wdir)
     >>> bunkers_storm_motion(p, u, v, h)
-    (<Quantity([22.73539654 10.27331352], 'knot')>,
-    <Quantity([ 7.27954821 34.99751866], 'knot')>,
-    <Quantity([15.00747237 22.63541609], 'knot')>)
+    (<Quantity([22.09618172 12.43406736], 'knot')>,
+    <Quantity([ 6.02861839 36.76517865], 'knot')>,
+    <Quantity([14.06240005 24.599623  ], 'knot')>)
 
     Notes
     -----
@@ -299,20 +299,24 @@ def bunkers_storm_motion(pressure, u, v, height):
 
     """
     # mean wind from sfc-6km
-    _, u_mean, v_mean = get_layer(pressure, u, v, height=height,
-                                  depth=units.Quantity(6000, 'meter'))
-    wind_mean = units.Quantity([np.mean(u_mean).m, np.mean(v_mean).m], u_mean.units)
+    wind_mean = weighted_continuous_average(pressure, u, v, height=height,
+                                            depth=units.Quantity(6000, 'meter'))
+
+    wind_mean = units.Quantity.from_list(wind_mean)
 
     # mean wind from sfc-500m
-    _, u_500m, v_500m = get_layer(pressure, u, v, height=height,
-                                  depth=units.Quantity(500, 'meter'))
-    wind_500m = units.Quantity([np.mean(u_500m).m, np.mean(v_500m).m], u_500m.units)
+    wind_500m = weighted_continuous_average(pressure, u, v, height=height,
+                                            depth=units.Quantity(500, 'meter'))
+
+    wind_500m = units.Quantity.from_list(wind_500m)
 
     # mean wind from 5.5-6km
-    _, u_5500m, v_5500m = get_layer(pressure, u, v, height=height,
-                                    depth=units.Quantity(500, 'meter'),
-                                    bottom=height[0] + units.Quantity(5500, 'meter'))
-    wind_5500m = units.Quantity([np.mean(u_5500m).m, np.mean(v_5500m).m], u_5500m.units)
+    wind_5500m = weighted_continuous_average(
+        pressure, u, v, height=height,
+        depth=units.Quantity(500, 'meter'),
+        bottom=height[0] + units.Quantity(5500, 'meter'))
+
+    wind_5500m = units.Quantity.from_list(wind_5500m)
 
     # Calculate the shear vector from sfc-500m to 5.5-6km
     shear = wind_5500m - wind_500m
@@ -320,7 +324,7 @@ def bunkers_storm_motion(pressure, u, v, height):
     # Take the cross product of the wind shear and k, and divide by the vector magnitude and
     # multiply by the deviation empirically calculated in Bunkers (2000) (7.5 m/s)
     shear_cross = concatenate([shear[1], -shear[0]])
-    shear_mag = units.Quantity(np.hypot(*(arg.magnitude for arg in shear)), shear.units)
+    shear_mag = np.hypot(*shear)
     rdev = shear_cross * (units.Quantity(7.5, 'm/s').to(u.units) / shear_mag)
 
     # Add the deviations to the layer average wind to get the RM motion


Test patch:
diff --git a/tests/calc/test_indices.py b/tests/calc/test_indices.py
index 744edeb63e..b2283c3e7d 100644
--- a/tests/calc/test_indices.py
+++ b/tests/calc/test_indices.py
@@ -171,8 +171,8 @@ def test_bunkers_motion():
     motion = concatenate(bunkers_storm_motion(data['pressure'],
                          data['u_wind'], data['v_wind'],
                          data['height']))
-    truth = [2.18346161, 0.86094706, 11.6006767, 12.53639395, 6.89206916,
-             6.69867051] * units('m/s')
+    truth = [2.062733, 0.96246913, 11.22554254, 12.83861839, 6.64413777,
+             6.90054376] * units('m/s')
     assert_almost_equal(motion.flatten(), truth, 8)