Repository: Unidata/MetPy

Issue:
get_layer heights unit issue
Using `get_layer` on some sounding data works oddly when asking for data in height. A request for data from 1-5km and 1000-5000m do not return the same data!

```python
import metpy.calc as mpcalc
from datetime import datetime
from metpy.units import units, pandas_dataframe_to_unit_arrays
from siphon.simplewebservice.wyoming import WyomingUpperAir

date = datetime(2021, 6, 27, 0)
station = 'FWD'

df = WyomingUpperAir.request_data(date, station)

data = pandas_dataframe_to_unit_arrays(df)

# Works as expected
p_subset, speed_subset, dir_subset = mpcalc.get_layer(data['pressure'],
                                                      data['speed'],
                                                      data['direction'],
                                                      bottom = 1000 * units.m,
                                                      depth = 5000 * units.m,
                                                      height = data['height'])

# Does not work as expected
p_subset, speed_subset, dir_subset = mpcalc.get_layer(data['pressure'],
                                                      data['speed'],
                                                      data['direction'],
                                                      bottom = 1 * units.km,
                                                      depth = 5 * units.km,
                                                      height = data['height'])

# Also works as expected when not providing heights from the data
p_subset, speed_subset, dir_subset = mpcalc.get_layer(data['pressure'],
                                                      data['speed'],
                                                      data['direction'],
                                                      bottom = 1 * units.km,
                                                      depth = 5 * units.km)
```

Discussion:

Patch:
diff --git a/src/metpy/calc/tools.py b/src/metpy/calc/tools.py
index 73d1de8863..91727ea35b 100644
--- a/src/metpy/calc/tools.py
+++ b/src/metpy/calc/tools.py
@@ -388,9 +388,8 @@ def _get_bound_pressure_height(pressure, bound, height=None, interpolate=True):
                     # Need to cast back to the input type since interp (up to at least numpy
                     # 1.13 always returns float64. This can cause upstream users problems,
                     # resulting in something like np.append() to upcast.
-                    bound_pressure = units.Quantity(np.interp(np.atleast_1d(bound.m), height.m,
-                                                    pressure.m).astype(np.result_type(bound)),
-                                                    pressure.units)
+                    bound_pressure = np.interp(np.atleast_1d(bound),
+                                               height, pressure).astype(np.result_type(bound))
                 else:
                     idx = (np.abs(height - bound)).argmin()
                     bound_pressure = pressure[idx]


Test patch:
diff --git a/tests/calc/test_calc_tools.py b/tests/calc/test_calc_tools.py
index d2876a51dd..b4eab64878 100644
--- a/tests/calc/test_calc_tools.py
+++ b/tests/calc/test_calc_tools.py
@@ -358,6 +358,20 @@ def test_get_layer(pressure, variable, heights, bottom, depth, interp, expected)
     assert_array_almost_equal(y_layer, expected[1], 3)
 
 
+def test_get_layer_units():
+    """Test get_layer when height profile has different units from bottom and depth."""
+    pressure, temperature = layer_test_data()
+    height = units.Quantity(np.linspace(100, 50000, len(pressure)), 'm')
+    pres_subset, temp_subset = get_layer(pressure, temperature,
+                                         bottom=units.Quantity(1, 'km'),
+                                         depth=units.Quantity(5, 'km'),
+                                         height=height)
+    pres_expected = units.Quantity([983, 900, 893], 'hPa')
+    temp_expected = units.Quantity([23.64385006, 16.66666667, 16.11422559], 'degC')
+    assert_array_almost_equal(pres_subset, pres_expected, 2)
+    assert_array_almost_equal(temp_subset, temp_expected, 3)
+
+
 def test_get_layer_masked():
     """Test get_layer with masked arrays as input."""
     p = units.Quantity(np.ma.array([1000, 500, 400]), 'hPa')