From c7b0b9b9db3b895dbc38fc231412d3b1c6baa2fd Mon Sep 17 00:00:00 2001
From: Ryan May <rmay@ucar.edu>
Date: Thu, 22 Apr 2021 14:12:49 -0600
Subject: [PATCH] BUG: Fix uses of unit multiplication internally

These are problematic with masked arrays. Also, using the Quantity()
constructor is faster across the board. Therefore, just get rid of all
the places where we multiply by units.

This also adds a test for get_layer where it was failing with masked
arrays, which is what originally motivated this.
---
 src/metpy/calc/basic.py          | 69 ++++++++++++------------
 src/metpy/calc/cross_sections.py |  8 +--
 src/metpy/calc/indices.py        | 61 ++++++++++++----------
 src/metpy/calc/kinematics.py     | 21 +++++---
 src/metpy/calc/thermo.py         | 90 +++++++++++++++++++-------------
 src/metpy/calc/tools.py          | 35 ++++++-------
 src/metpy/constants.py           | 44 ++++++++--------
 src/metpy/io/metar.py            | 28 +++++-----
 src/metpy/io/station_data.py     |  3 +-
 src/metpy/plots/declarative.py   | 11 ++--
 src/metpy/plots/skewt.py         | 28 +++++-----
 src/metpy/testing.py             | 16 +++---
 src/metpy/units.py               |  2 +-
 src/metpy/xarray.py              | 16 +++---
 tests/calc/test_calc_tools.py    | 11 ++++
 15 files changed, 241 insertions(+), 202 deletions(-)

diff --git a/src/metpy/calc/basic.py b/src/metpy/calc/basic.py
index 4e8a131025f..71e37d6bd9d 100644
--- a/src/metpy/calc/basic.py
+++ b/src/metpy/calc/basic.py
@@ -24,9 +24,9 @@
 exporter = Exporter(globals())
 
 # The following variables are constants for a standard atmosphere
-t0 = 288. * units.kelvin
-p0 = 1013.25 * units.hPa
-gamma = 6.5 * units('K/km')
+t0 = units.Quantity(288., 'kelvin')
+p0 = units.Quantity(1013.25, 'hPa')
+gamma = units.Quantity(6.5, 'K/km')
 
 
 @exporter.export
@@ -89,24 +89,24 @@ def wind_direction(u, v, convention='from'):
     of 0.
 
     """
-    wdir = 90. * units.deg - np.arctan2(-v, -u)
+    wdir = units.Quantity(90., 'deg') - np.arctan2(-v, -u)
     origshape = wdir.shape
     wdir = np.atleast_1d(wdir)
 
     # Handle oceanographic convection
     if convention == 'to':
-        wdir -= 180 * units.deg
+        wdir -= units.Quantity(180., 'deg')
     elif convention not in ('to', 'from'):
         raise ValueError('Invalid kwarg for "convention". Valid options are "from" or "to".')
 
     mask = wdir <= 0
     if np.any(mask):
-        wdir[mask] += 360. * units.deg
+        wdir[mask] += units.Quantity(360., 'deg')
     # avoid unintended modification of `pint.Quantity` by direct use of magnitude
     calm_mask = (np.asarray(u.magnitude) == 0.) & (np.asarray(v.magnitude) == 0.)
     # np.any check required for legacy numpy which treats 0-d False boolean index as zero
     if np.any(calm_mask):
-        wdir[calm_mask] = 0. * units.deg
+        wdir[calm_mask] = units.Quantity(0., 'deg')
     return wdir.reshape(origshape).to('degrees')
 
 
@@ -199,7 +199,7 @@ def windchill(temperature, speed, face_level_winds=False, mask_undefined=True):
         # noinspection PyAugmentAssignment
         speed = speed * 1.5
 
-    temp_limit, speed_limit = 10. * units.degC, 3 * units.mph
+    temp_limit, speed_limit = units.Quantity(10., 'degC'), units.Quantity(3, 'mph')
     speed_factor = speed.to('km/hr').magnitude ** 0.16
     wcti = units.Quantity((0.6215 + 0.3965 * speed_factor) * temperature.to('degC').magnitude
                           - 11.37 * speed_factor + 13.12, units.degC).to(temperature.units)
@@ -257,38 +257,39 @@ def heat_index(temperature, relative_humidity, mask_undefined=True):
     relative_humidity = np.atleast_1d(relative_humidity)
     # assign units to relative_humidity if they currently are not present
     if not hasattr(relative_humidity, 'units'):
-        relative_humidity = relative_humidity * units.dimensionless
-    delta = temperature.to(units.degF) - 0. * units.degF
+        relative_humidity = units.Quantity(relative_humidity, 'dimensionless')
+    delta = temperature.to(units.degF) - units.Quantity(0., 'degF')
     rh2 = relative_humidity * relative_humidity
     delta2 = delta * delta
 
     # Simplifed Heat Index -- constants converted for relative_humidity in [0, 1]
-    a = -10.3 * units.degF + 1.1 * delta + 4.7 * units.delta_degF * relative_humidity
+    a = (units.Quantity(-10.3, 'degF') + 1.1 * delta
+         + units.Quantity(4.7, 'delta_degF') * relative_humidity)
 
     # More refined Heat Index -- constants converted for relative_humidity in [0, 1]
-    b = (-42.379 * units.degF
+    b = (units.Quantity(-42.379, 'degF')
          + 2.04901523 * delta
-         + 1014.333127 * units.delta_degF * relative_humidity
+         + units.Quantity(1014.333127, 'delta_degF') * relative_humidity
          - 22.475541 * delta * relative_humidity
-         - 6.83783e-3 / units.delta_degF * delta2
-         - 5.481717e2 * units.delta_degF * rh2
-         + 1.22874e-1 / units.delta_degF * delta2 * relative_humidity
+         - units.Quantity(6.83783e-3, '1/delta_degF') * delta2
+         - units.Quantity(5.481717e2, 'delta_degF') * rh2
+         + units.Quantity(1.22874e-1, '1/delta_degF') * delta2 * relative_humidity
          + 8.5282 * delta * rh2
-         - 1.99e-2 / units.delta_degF * delta2 * rh2)
+         - units.Quantity(1.99e-2, '1/delta_degF') * delta2 * rh2)
 
     # Create return heat index
-    hi = np.full(np.shape(temperature), np.nan) * units.degF
+    hi = units.Quantity(np.full(np.shape(temperature), np.nan), 'degF')
     # Retain masked status of temperature with resulting heat index
     if hasattr(temperature, 'mask'):
         hi = masked_array(hi)
 
     # If T <= 40F, Heat Index is T
-    sel = (temperature <= 40. * units.degF)
+    sel = (temperature <= units.Quantity(40., 'degF'))
     if np.any(sel):
         hi[sel] = temperature[sel].to(units.degF)
 
     # If a < 79F and hi is unset, Heat Index is a
-    sel = (a < 79. * units.degF) & np.isnan(hi)
+    sel = (a < units.Quantity(79., 'degF')) & np.isnan(hi)
     if np.any(sel):
         hi[sel] = a[sel]
 
@@ -298,26 +299,28 @@ def heat_index(temperature, relative_humidity, mask_undefined=True):
         hi[sel] = b[sel]
 
     # Adjustment for RH <= 13% and 80F <= T <= 112F
-    sel = ((relative_humidity <= 13. * units.percent) & (temperature >= 80. * units.degF)
-           & (temperature <= 112. * units.degF))
+    sel = ((relative_humidity <= units.Quantity(13., 'percent'))
+           & (temperature >= units.Quantity(80., 'degF'))
+           & (temperature <= units.Quantity(112., 'degF')))
     if np.any(sel):
         rh15adj = ((13. - relative_humidity[sel] * 100.) / 4.
-                   * np.sqrt((17. * units.delta_degF
-                              - np.abs(delta[sel] - 95. * units.delta_degF))
-                             / 17. * units.delta_degF))
-        hi[sel] = hi[sel] - rh15adj
+                   * np.sqrt((units.Quantity(17., 'delta_degF')
+                              - np.abs(delta[sel] - units.Quantity(95., 'delta_degF')))
+                             / units.Quantity(17., 'delta_degF')))
+        hi[sel] = hi[sel] - units.Quantity(rh15adj, 'delta_degF')
 
     # Adjustment for RH > 85% and 80F <= T <= 87F
-    sel = ((relative_humidity > 85. * units.percent) & (temperature >= 80. * units.degF)
-           & (temperature <= 87. * units.degF))
+    sel = ((relative_humidity > units.Quantity(85., 'percent'))
+           & (temperature >= units.Quantity(80., 'degF'))
+           & (temperature <= units.Quantity(87., 'degF')))
     if np.any(sel):
-        rh85adj = 0.02 * (relative_humidity[sel] * 100. - 85.) * (87. * units.delta_degF
-                                                                  - delta[sel])
+        rh85adj = (0.02 * (relative_humidity[sel] * 100. - 85.)
+                   * (units.Quantity(87., 'delta_degF') - delta[sel]))
         hi[sel] = hi[sel] + rh85adj
 
     # See if we need to mask any undefined values
     if mask_undefined:
-        mask = np.array(temperature < 80. * units.degF)
+        mask = np.array(temperature < units.Quantity(80., 'degF'))
         if mask.any():
             hi = masked_array(hi, mask=mask)
     return hi
@@ -398,7 +401,7 @@ def apparent_temperature(temperature, relative_humidity, speed, face_level_winds
     # If no values are masked and provided temperature does not have a mask
     # we should return a non-masked array
     if not np.any(app_temperature.mask) and not hasattr(temperature, 'mask'):
-        app_temperature = np.array(app_temperature.m) * temperature.units
+        app_temperature = units.Quantity(np.array(app_temperature.m), temperature.units)
 
     if is_not_scalar:
         return app_temperature
@@ -1103,7 +1106,7 @@ def altimeter_to_station_pressure(altimeter_value, height):
 
     return ((altimeter_value ** n
              - ((p0.to(altimeter_value.units) ** n * gamma * height) / t0)) ** (1 / n)
-            + 0.3 * units.hPa)
+            + units.Quantity(0.3, 'hPa'))
 
 
 @exporter.export
diff --git a/src/metpy/calc/cross_sections.py b/src/metpy/calc/cross_sections.py
index eb05dd830ad..6d5cbe40b29 100644
--- a/src/metpy/calc/cross_sections.py
+++ b/src/metpy/calc/cross_sections.py
@@ -48,8 +48,8 @@ def distances_from_cross_section(cross):
         y = distance * np.cos(np.deg2rad(forward_az))
 
         # Build into DataArrays
-        x = xr.DataArray(x * units.meter, coords=lon.coords, dims=lon.dims)
-        y = xr.DataArray(y * units.meter, coords=lat.coords, dims=lat.dims)
+        x = xr.DataArray(units.Quantity(x, 'meter'), coords=lon.coords, dims=lon.dims)
+        y = xr.DataArray(units.Quantity(y, 'meter'), coords=lat.coords, dims=lat.dims)
 
     elif check_axis(cross.metpy.x, 'x') and check_axis(cross.metpy.y, 'y'):
 
@@ -88,8 +88,8 @@ def latitude_from_cross_section(cross):
             inverse=True,
             radians=False
         )[1]
-        latitude = xr.DataArray(latitude * units.degrees_north, coords=y.coords, dims=y.dims)
-        return latitude
+        return xr.DataArray(units.Quantity(latitude, 'degrees_north'), coords=y.coords,
+                            dims=y.dims)
 
 
 @exporter.export
diff --git a/src/metpy/calc/indices.py b/src/metpy/calc/indices.py
index 7277485e67e..a1c975b0978 100644
--- a/src/metpy/calc/indices.py
+++ b/src/metpy/calc/indices.py
@@ -78,8 +78,7 @@ def precipitable_water(pressure, dewpoint, *, bottom=None, top=None):
     w = mixing_ratio(saturation_vapor_pressure(dewpoint_layer), pres_layer)
 
     # Since pressure is in decreasing order, pw will be the opposite sign of that expected.
-    pw = -1. * (np.trapz(w.magnitude, pres_layer.magnitude) * (w.units * pres_layer.units)
-                / (mpconsts.g * mpconsts.rho_l))
+    pw = -np.trapz(w, pres_layer) / (mpconsts.g * mpconsts.rho_l)
     return pw.to('millimeters')
 
 
@@ -187,26 +186,29 @@ 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=6000 * units('meter'))
-    wind_mean = [np.mean(u_mean).m, np.mean(v_mean).m] * u_mean.units
+    _, 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)
 
     # mean wind from sfc-500m
-    _, u_500m, v_500m = get_layer(pressure, u, v, height=height, depth=500 * units('meter'))
-    wind_500m = [np.mean(u_500m).m, np.mean(v_500m).m] * u_500m.units
+    _, 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)
 
     # mean wind from 5.5-6km
-    _, u_5500m, v_5500m = get_layer(pressure, u, v, height=height, depth=500 * units('meter'),
-                                    bottom=height[0] + 5500 * units('meter'))
-    wind_5500m = [np.mean(u_5500m).m, np.mean(v_5500m).m] * u_5500m.units
+    _, 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)
 
     # Calculate the shear vector from sfc-500m to 5.5-6km
     shear = wind_5500m - wind_500m
 
     # Take the cross product of the wind shear and k, and divide by the vector magnitude and
-    # multiply by the deviaton empirically calculated in Bunkers (2000) (7.5 m/s)
+    # multiply by the deviation empirically calculated in Bunkers (2000) (7.5 m/s)
     shear_cross = concatenate([shear[1], -shear[0]])
-    shear_mag = np.hypot(*(arg.magnitude for arg in shear)) * shear.units
-    rdev = shear_cross * (7.5 * units('m/s').to(u.units) / shear_mag)
+    shear_mag = units.Quantity(np.hypot(*(arg.magnitude for arg in shear)), shear.units)
+    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
     right_mover = wind_mean + rdev
@@ -309,12 +311,12 @@ def supercell_composite(mucape, effective_storm_helicity, effective_shear):
         Supercell composite
 
     """
-    effective_shear = np.clip(np.atleast_1d(effective_shear), None, 20 * units('m/s'))
-    effective_shear[effective_shear < 10 * units('m/s')] = 0 * units('m/s')
-    effective_shear = effective_shear / (20 * units('m/s'))
+    effective_shear = np.clip(np.atleast_1d(effective_shear), None, units.Quantity(20, 'm/s'))
+    effective_shear[effective_shear < units.Quantity(10, 'm/s')] = units.Quantity(0, 'm/s')
+    effective_shear = effective_shear / units.Quantity(20, 'm/s')
 
-    return ((mucape / (1000 * units('J/kg')))
-            * (effective_storm_helicity / (50 * units('m^2/s^2')))
+    return ((mucape / units.Quantity(1000, 'J/kg'))
+            * (effective_storm_helicity / units.Quantity(50, 'm^2/s^2'))
             * effective_shear).to('dimensionless')
 
 
@@ -361,17 +363,18 @@ def significant_tornado(sbcape, surface_based_lcl_height, storm_helicity_1km, sh
 
     """
     surface_based_lcl_height = np.clip(np.atleast_1d(surface_based_lcl_height),
-                                       1000 * units.m, 2000 * units.m)
-    surface_based_lcl_height[surface_based_lcl_height > 2000 * units.m] = 0 * units.m
-    surface_based_lcl_height = ((2000. * units.m - surface_based_lcl_height)
-                                / (1000. * units.m))
-    shear_6km = np.clip(np.atleast_1d(shear_6km), None, 30 * units('m/s'))
-    shear_6km[shear_6km < 12.5 * units('m/s')] = 0 * units('m/s')
-    shear_6km /= 20 * units('m/s')
-
-    return ((sbcape / (1500. * units('J/kg')))
+                                       units.Quantity(1000., 'm'), units.Quantity(2000., 'm'))
+    surface_based_lcl_height[surface_based_lcl_height > units.Quantity(2000., 'm')] = \
+        units.Quantity(0, 'm')
+    surface_based_lcl_height = ((units.Quantity(2000., 'm') - surface_based_lcl_height)
+                                / units.Quantity(1000., 'm'))
+    shear_6km = np.clip(np.atleast_1d(shear_6km), None, units.Quantity(30, 'm/s'))
+    shear_6km[shear_6km < units.Quantity(12.5, 'm/s')] = units.Quantity(0, 'm/s')
+    shear_6km /= units.Quantity(20, 'm/s')
+
+    return ((sbcape / units.Quantity(1500., 'J/kg'))
             * surface_based_lcl_height
-            * (storm_helicity_1km / (150. * units('m^2/s^2')))
+            * (storm_helicity_1km / units.Quantity(150., 'm^2/s^2'))
             * shear_6km)
 
 
@@ -436,7 +439,7 @@ def critical_angle(pressure, u, v, height, u_storm, v_storm):
     v = v[sort_inds]
 
     # Calculate sfc-500m shear vector
-    shr5 = bulk_shear(pressure, u, v, height=height, depth=500 * units('meter'))
+    shr5 = bulk_shear(pressure, u, v, height=height, depth=units.Quantity(500, 'meter'))
 
     # Make everything relative to the sfc wind orientation
     umn = u_storm - u[0]
@@ -445,6 +448,6 @@ def critical_angle(pressure, u, v, height, u_storm, v_storm):
     vshr = np.asarray([shr5[0].magnitude, shr5[1].magnitude])
     vsm = np.asarray([umn.magnitude, vmn.magnitude])
     angle_c = np.dot(vshr, vsm) / (np.linalg.norm(vshr) * np.linalg.norm(vsm))
-    critical_angle = np.arccos(angle_c) * units('radian')
+    critical_angle = units.Quantity(np.arccos(angle_c), 'radian')
 
     return critical_angle.to('degrees')
diff --git a/src/metpy/calc/kinematics.py b/src/metpy/calc/kinematics.py
index 5e2281cdefc..5de5a8501a7 100644
--- a/src/metpy/calc/kinematics.py
+++ b/src/metpy/calc/kinematics.py
@@ -563,8 +563,7 @@ def montgomery_streamfunction(height, temperature):
 @preprocess_and_wrap()
 @check_units('[length]', '[speed]', '[speed]', '[length]',
              bottom='[length]', storm_u='[speed]', storm_v='[speed]')
-def storm_relative_helicity(height, u, v, depth, *, bottom=0 * units.m,
-                            storm_u=0 * units('m/s'), storm_v=0 * units('m/s')):
+def storm_relative_helicity(height, u, v, depth, *, bottom=None, storm_u=None, storm_v=None):
     # Partially adapted from similar SharpPy code
     r"""Calculate storm relative helicity.
 
@@ -622,6 +621,13 @@ def storm_relative_helicity(height, u, v, depth, *, bottom=0 * units.m,
        ``storm_v`` parameters to keyword-only arguments
 
     """
+    if bottom is None:
+        bottom = units.Quantity(0, 'm')
+    if storm_u is None:
+        storm_u = units.Quantity(0, 'm/s')
+    if storm_v is None:
+        storm_v = units.Quantity(0, 'm/s')
+
     _, u, v = get_layer_heights(height, depth, u, v, with_agl=True, bottom=bottom)
 
     storm_relative_u = u - storm_u
@@ -634,10 +640,10 @@ def storm_relative_helicity(height, u, v, depth, *, bottom=0 * units.m,
     # mask will return a masked value rather than 0. See numpy/numpy#11736
     positive_srh = int_layers[int_layers.magnitude > 0.].sum()
     if np.ma.is_masked(positive_srh):
-        positive_srh = 0.0 * units('meter**2 / second**2')
+        positive_srh = units.Quantity(0.0, 'meter**2 / second**2')
     negative_srh = int_layers[int_layers.magnitude < 0.].sum()
     if np.ma.is_masked(negative_srh):
-        negative_srh = 0.0 * units('meter**2 / second**2')
+        negative_srh = units.Quantity(0.0, 'meter**2 / second**2')
 
     return (positive_srh.to('meter ** 2 / second ** 2'),
             negative_srh.to('meter ** 2 / second ** 2'),
@@ -789,8 +795,8 @@ def potential_vorticity_baroclinic(
         (np.shape(potential_temperature)[y_dim] == 1)
         and (np.shape(potential_temperature)[x_dim] == 1)
     ):
-        dthtady = 0 * units.K / units.m  # axis=y_dim only has one dimension
-        dthtadx = 0 * units.K / units.m  # axis=x_dim only has one dimension
+        dthtady = units.Quantity(0, 'K/m')  # axis=y_dim only has one dimension
+        dthtadx = units.Quantity(0, 'K/m')  # axis=x_dim only has one dimension
     else:
         dthtady = first_derivative(potential_temperature, delta=dy, axis=y_dim)
         dthtadx = first_derivative(potential_temperature, delta=dx, axis=x_dim)
@@ -798,8 +804,7 @@ def potential_vorticity_baroclinic(
     dvdp = first_derivative(v, x=pressure, axis=vertical_dim)
 
     return (-mpconsts.g * (dudp * dthtady - dvdp * dthtadx
-                           + avor * dthtadp)).to(units.kelvin * units.meter**2
-                                                 / (units.second * units.kilogram))
+                           + avor * dthtadp)).to('K * m**2 / (s * kg)')
 
 
 @exporter.export
diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py
index 4357c580e9d..4bb2791eaba 100644
--- a/src/metpy/calc/thermo.py
+++ b/src/metpy/calc/thermo.py
@@ -21,7 +21,7 @@
 
 exporter = Exporter(globals())
 
-sat_pressure_0c = 6.112 * units.millibar
+sat_pressure_0c = units.Quantity(6.112, 'millibar')
 
 
 @exporter.export
@@ -408,7 +408,8 @@ def _lcl_iter(p, p0, w, t):
 
     # np.isclose needed if surface is LCL due to precision error with np.log in dewpoint.
     # Causes issues with parcel_profile_with_lcl if removed. Issue #1187
-    lcl_p = np.where(np.isclose(lcl_p, pressure.m), pressure.m, lcl_p) * pressure.units
+    lcl_p = units.Quantity(np.where(np.isclose(lcl_p, pressure.m), pressure.m, lcl_p),
+                           pressure.units)
 
     return lcl_p, globals()['dewpoint'](vapor_pressure(lcl_p, w)).to(temperature.units)
 
@@ -508,7 +509,8 @@ def lfc(pressure, temperature, dewpoint, parcel_temperature_profile=None, dewpoi
         mask = pressure < this_lcl[0]
         if np.all(_less_or_close(parcel_temperature_profile[mask], temperature[mask])):
             # LFC doesn't exist
-            x, y = np.nan * pressure.units, np.nan * temperature.units
+            x = units.Quantity(np.nan, pressure.units)
+            y = units.Quantity(np.nan, temperature.units)
         else:  # LFC = LCL
             x, y = this_lcl
         return x, y
@@ -522,7 +524,8 @@ def lfc(pressure, temperature, dewpoint, parcel_temperature_profile=None, dewpoi
                                             temperature[1:], direction='decreasing',
                                             log_x=True)
             if np.min(el_pres) > this_lcl[0]:
-                x, y = np.nan * pressure.units, np.nan * temperature.units
+                x = units.Quantity(np.nan, pressure.units)
+                y = units.Quantity(np.nan, temperature.units)
             else:
                 x, y = this_lcl
             return x, y
@@ -666,7 +669,8 @@ def el(pressure, temperature, dewpoint, parcel_temperature_profile=None, which='
 
     # If the top of the sounding parcel is warmer than the environment, there is no EL
     if parcel_temperature_profile[-1] > temperature[-1]:
-        return np.nan * pressure.units, np.nan * temperature.units
+        return (units.Quantity(np.nan, pressure.units),
+                units.Quantity(np.nan, temperature.units))
 
     # Interpolate in log space to find the appropriate pressure - units have to be stripped
     # and reassigned to allow np.log() to function properly.
@@ -679,7 +683,8 @@ def el(pressure, temperature, dewpoint, parcel_temperature_profile=None, which='
                                         parcel_temperature_profile, temperature, dewpoint,
                                         intersect_type='EL')
     else:
-        return np.nan * pressure.units, np.nan * temperature.units
+        return (units.Quantity(np.nan, pressure.units),
+                units.Quantity(np.nan, temperature.units))
 
 
 @exporter.export
@@ -896,7 +901,7 @@ def _insert_lcl_level(pressure, temperature, lcl_pressure):
     # Pressure needs to be increasing for searchsorted, so flip it and then convert
     # the index back to the original array
     loc = pressure.size - pressure[::-1].searchsorted(lcl_pressure)
-    return temperature.units * np.insert(temperature.m, loc, interp_temp.m)
+    return units.Quantity(np.insert(temperature.m, loc, interp_temp.m), temperature.units)
 
 
 @exporter.export
@@ -971,8 +976,8 @@ def saturation_vapor_pressure(temperature):
     """
     # Converted from original in terms of C to use kelvin. Using raw absolute values of C in
     # a formula plays havoc with units support.
-    return sat_pressure_0c * np.exp(17.67 * (temperature - 273.15 * units.kelvin)
-                                    / (temperature - 29.65 * units.kelvin))
+    return sat_pressure_0c * np.exp(17.67 * (temperature - units.Quantity(273.15, 'kelvin'))
+                                    / (temperature - units.Quantity(29.65, 'kelvin')))
 
 
 @exporter.export
@@ -1041,7 +1046,8 @@ def dewpoint(vapor_pressure):
 
     """
     val = np.log(vapor_pressure / sat_pressure_0c)
-    return 0. * units.degC + 243.5 * units.delta_degC * val / (17.67 - val)
+    return (units.Quantity(0., 'degC')
+            + units.Quantity(243.5, 'delta_degC') * val / (17.67 - val))
 
 
 @exporter.export
@@ -1380,7 +1386,7 @@ def density(pressure, temperature, mixing_ratio, molecular_weight_ratio=mpconsts
 
     """
     virttemp = virtual_temperature(temperature, mixing_ratio, molecular_weight_ratio)
-    return (pressure / (mpconsts.Rd * virttemp)).to(units.kilogram / units.meter ** 3)
+    return (pressure / (mpconsts.Rd * virttemp)).to('kg/m**3')
 
 
 @exporter.export
@@ -1441,7 +1447,7 @@ def relative_humidity_wet_psychrometric(pressure, dry_bulb_temperature, wet_bulb
 )
 @check_units('[pressure]', '[temperature]', '[temperature]')
 def psychrometric_vapor_pressure_wet(pressure, dry_bulb_temperature, wet_bulb_temperature,
-                                     psychrometer_coefficient=6.21e-4 / units.kelvin):
+                                     psychrometer_coefficient=None):
     r"""Calculate the vapor pressure with wet bulb and dry bulb temperatures.
 
     This uses a psychrometric relationship as outlined in [WMO8]_, with
@@ -1490,6 +1496,8 @@ def psychrometric_vapor_pressure_wet(pressure, dry_bulb_temperature, wet_bulb_te
     saturation_vapor_pressure
 
     """
+    if psychrometer_coefficient is None:
+        psychrometer_coefficient = units.Quantity(6.21e-4, '1/K')
     return (saturation_vapor_pressure(wet_bulb_temperature) - psychrometer_coefficient
             * pressure * (dry_bulb_temperature - wet_bulb_temperature).to('kelvin'))
 
@@ -1788,7 +1796,7 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom'
 
     # If there is no LFC, no need to proceed.
     if np.isnan(lfc_pressure):
-        return 0 * units('J/kg'), 0 * units('J/kg')
+        return units.Quantity(0, 'J/kg'), units.Quantity(0, 'J/kg')
     else:
         lfc_pressure = lfc_pressure.magnitude
 
@@ -1814,7 +1822,7 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom'
     x_clipped = x[p_mask].magnitude
     y_clipped = y[p_mask].magnitude
     cape = (mpconsts.Rd
-            * (np.trapz(y_clipped, np.log(x_clipped)) * units.degK)).to(units('J/kg'))
+            * units.Quantity(np.trapz(y_clipped, np.log(x_clipped)), 'K')).to(units('J/kg'))
 
     # CIN
     # Only use data between the surface and LFC for calculation
@@ -1822,11 +1830,11 @@ def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom'
     x_clipped = x[p_mask].magnitude
     y_clipped = y[p_mask].magnitude
     cin = (mpconsts.Rd
-           * (np.trapz(y_clipped, np.log(x_clipped)) * units.degK)).to(units('J/kg'))
+           * units.Quantity(np.trapz(y_clipped, np.log(x_clipped)), 'K')).to(units('J/kg'))
 
     # Set CIN to 0 if it's returned as a positive value (#1190)
-    if cin > 0 * units('J/kg'):
-        cin = 0 * units('J/kg')
+    if cin > units.Quantity(0, 'J/kg'):
+        cin = units.Quantity(0, 'J/kg')
     return cape, cin
 
 
@@ -1853,7 +1861,8 @@ def _find_append_zero_crossings(x, y):
         Y values of data
 
     """
-    crossings = find_intersections(x[1:], y[1:], y.units * np.zeros_like(y[1:]), log_x=True)
+    crossings = find_intersections(x[1:], y[1:],
+                                   units.Quantity(np.zeros_like(y[1:]), y.units), log_x=True)
     x = concatenate((x, crossings[0]))
     y = concatenate((y, crossings[1]))
 
@@ -1872,8 +1881,8 @@ def _find_append_zero_crossings(x, y):
 @exporter.export
 @preprocess_and_wrap()
 @check_units('[pressure]', '[temperature]', '[temperature]')
-def most_unstable_parcel(pressure, temperature, dewpoint, height=None,
-                         bottom=None, depth=300 * units.hPa):
+def most_unstable_parcel(pressure, temperature, dewpoint, height=None, bottom=None,
+                         depth=None):
     """
     Determine the most unstable parcel in a layer.
 
@@ -1924,6 +1933,8 @@ def most_unstable_parcel(pressure, temperature, dewpoint, height=None,
        Renamed ``heights`` parameter to ``height``
 
     """
+    if depth is None:
+        depth = units.Quantity(300, 'hPa')
     p_layer, t_layer, td_layer = get_layer(pressure, temperature, dewpoint, bottom=bottom,
                                            depth=depth, height=height, interpolate=False)
     theta_e = equivalent_potential_temperature(p_layer, t_layer, td_layer)
@@ -2016,7 +2027,8 @@ def _isen_iter(iter_log_p, isentlevs_nd, ka, a, b, pok):
     slices = [np.newaxis] * ndim
     slices[vertical_dim] = slice(None)
     slices = tuple(slices)
-    pres = np.broadcast_to(pres[slices].magnitude, temperature.shape) * pres.units
+    pres = units.Quantity(np.broadcast_to(pres[slices].magnitude, temperature.shape),
+                          pres.units)
 
     # Sort input data
     sort_pres = np.argsort(pres.m, axis=vertical_dim)
@@ -2077,11 +2089,11 @@ def _isen_iter(iter_log_p, isentlevs_nd, ka, a, b, pok):
     isentprs[~(good & _less_or_close(isentprs, np.max(pres.m)))] = np.nan
 
     # create list for storing output data
-    ret = [isentprs * units.hPa]
+    ret = [units.Quantity(isentprs, 'hPa')]
 
     # if temperature_out = true, calculate temperature and output as last item in list
     if temperature_out:
-        ret.append((isentlevs_nd / ((mpconsts.P0.m / isentprs) ** ka)) * units.kelvin)
+        ret.append(units.Quantity((isentlevs_nd / ((mpconsts.P0.m / isentprs) ** ka)), 'K'))
 
     # do an interpolation for each additional argument
     if args:
@@ -2350,7 +2362,7 @@ def mixed_layer_cape_cin(pressure, temperature, dewpoint, **kwargs):
     Quantities even when given xarray DataArray profiles.
 
     """
-    depth = kwargs.get('depth', 100 * units.hPa)
+    depth = kwargs.get('depth', units.Quantity(100, 'hPa'))
     parcel_pressure, parcel_temp, parcel_dewpoint = mixed_parcel(pressure, temperature,
                                                                  dewpoint, **kwargs)
 
@@ -2370,7 +2382,7 @@ def mixed_layer_cape_cin(pressure, temperature, dewpoint, **kwargs):
 @preprocess_and_wrap()
 @check_units('[pressure]', '[temperature]', '[temperature]')
 def mixed_parcel(pressure, temperature, dewpoint, parcel_start_pressure=None,
-                 height=None, bottom=None, depth=100 * units.hPa, interpolate=True):
+                 height=None, bottom=None, depth=None, interpolate=True):
     r"""Calculate the properties of a parcel mixed from a layer.
 
     Determines the properties of an air parcel that is the result of complete mixing of a
@@ -2429,6 +2441,9 @@ def mixed_parcel(pressure, temperature, dewpoint, parcel_start_pressure=None,
     if not parcel_start_pressure:
         parcel_start_pressure = pressure[0]
 
+    if depth is None:
+        depth = units.Quantity(100, 'hPa')
+
     # Calculate the potential temperature and mixing ratio over the layer
     theta = potential_temperature(pressure, temperature)
     mixing_ratio = saturation_mixing_ratio(pressure, dewpoint)
@@ -2455,8 +2470,7 @@ def mixed_parcel(pressure, temperature, dewpoint, parcel_start_pressure=None,
 @exporter.export
 @preprocess_and_wrap()
 @check_units('[pressure]')
-def mixed_layer(pressure, *args, height=None, bottom=None, depth=100 * units.hPa,
-                interpolate=True):
+def mixed_layer(pressure, *args, height=None, bottom=None, depth=None, interpolate=True):
     r"""Mix variable(s) over a layer, yielding a mass-weighted average.
 
     This function will integrate a data variable with respect to pressure and determine the
@@ -2499,6 +2513,8 @@ def mixed_layer(pressure, *args, height=None, bottom=None, depth=100 * units.hPa
        Renamed ``p``, ``heights`` parameters to ``pressure``, ``height``
 
     """
+    if depth is None:
+        depth = units.Quantity(100, 'hPa')
     layer = get_layer(pressure, *args, height=height, bottom=bottom,
                       depth=depth, interpolate=interpolate)
     p_layer = layer[0]
@@ -2507,8 +2523,8 @@ def mixed_layer(pressure, *args, height=None, bottom=None, depth=100 * units.hPa
     ret = []
     for datavar_layer in datavars_layer:
         actual_depth = abs(p_layer[0] - p_layer[-1])
-        ret.append((-1. / actual_depth.m) * np.trapz(datavar_layer.m, p_layer.m)
-                   * datavar_layer.units)
+        ret.append(units.Quantity(np.trapz(datavar_layer.m, p_layer.m) / -actual_depth.m,
+                   datavar_layer.units))
     return ret
 
 
@@ -2675,8 +2691,9 @@ def thickness_hydrostatic(pressure, temperature, mixing_ratio=None,
             layer_virttemp = virtual_temperature(layer_temp, layer_w, molecular_weight_ratio)
 
     # Take the integral (with unit handling) and return the result in meters
-    return (- mpconsts.Rd / mpconsts.g * np.trapz(
-        layer_virttemp.m_as('K'), x=np.log(layer_p.m_as('hPa'))) * units.K).to('m')
+    integral = units.Quantity(np.trapz(layer_virttemp.m_as('K'), np.log(layer_p.m_as('hPa'))),
+                              units.K)
+    return (-mpconsts.Rd / mpconsts.g * integral).to('m')
 
 
 @exporter.export
@@ -2949,17 +2966,16 @@ def wet_bulb_temperature(pressure, temperature, dewpoint):
                    flags=['buffered'])
 
     for press, lpress, ltemp, ret in it:
-        press = press * pressure.units
-        lpress = lpress * lcl_press.units
-        ltemp = ltemp * lcl_temp.units
-        moist_adiabat_temperatures = moist_lapse(press, ltemp, lpress)
+        moist_adiabat_temperatures = moist_lapse(units.Quantity(press, pressure.units),
+                                                 units.Quantity(ltemp, lcl_temp.units),
+                                                 units.Quantity(lpress, lcl_press.units))
         ret[...] = moist_adiabat_temperatures.magnitude
 
     # If we started with a scalar, return a scalar
     ret = it.operands[3]
     if ret.size == 1:
         ret = ret[0]
-    return ret * moist_adiabat_temperatures.units
+    return units.Quantity(ret, moist_adiabat_temperatures.units)
 
 
 @exporter.export
@@ -3210,7 +3226,7 @@ def lifted_index(pressure, temperature, parcel_profile):
 
     """
     # find the index for the 500 hPa pressure level.
-    idx = np.where(pressure == 500 * units.hPa)
+    idx = np.where(pressure == units.Quantity(500, 'hPa'))
     # find the measured temperature at 500 hPa.
     T500 = temperature[idx]
     # find the parcel profile temperature at 500 hPa.
diff --git a/src/metpy/calc/tools.py b/src/metpy/calc/tools.py
index 3e10d89d065..5b92f08432b 100644
--- a/src/metpy/calc/tools.py
+++ b/src/metpy/calc/tools.py
@@ -30,8 +30,8 @@
     UND
 )  # note the order matters!
 
-MAX_DEGREE_ANGLE = 360 * units.degree
-BASE_DEGREE_MULTIPLIER = 22.5 * units.degree
+MAX_DEGREE_ANGLE = units.Quantity(360, 'degree')
+BASE_DEGREE_MULTIPLIER = units.Quantity(22.5, 'degree')
 
 DIR_DICT = {dir_str: i * BASE_DEGREE_MULTIPLIER for i, dir_str in enumerate(DIR_STRS)}
 DIR_DICT[UND] = np.nan
@@ -388,9 +388,9 @@ 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 = (np.interp(np.atleast_1d(bound.m), height.m,
-                                                pressure.m).astype(np.result_type(bound))
-                                      * pressure.units)
+                    bound_pressure = units.Quantity(np.interp(np.atleast_1d(bound.m), height.m,
+                                                    pressure.m).astype(np.result_type(bound)),
+                                                    pressure.units)
                 else:
                     idx = (np.abs(height - bound)).argmin()
                     bound_pressure = pressure[idx]
@@ -521,8 +521,7 @@ def get_layer_heights(height, depth, *args, bottom=None, interpolate=True, with_
 @exporter.export
 @preprocess_and_wrap()
 @check_units('[pressure]')
-def get_layer(pressure, *args, height=None, bottom=None, depth=100 * units.hPa,
-              interpolate=True):
+def get_layer(pressure, *args, height=None, bottom=None, depth=None, interpolate=True):
     r"""Return an atmospheric layer from upper air data with the requested bottom and depth.
 
     This function will subset an upper air dataset to contain only the specified layer. The
@@ -566,7 +565,7 @@ def get_layer(pressure, *args, height=None, bottom=None, depth=100 * units.hPa,
     """
     # If we get the depth kwarg, but it's None, set it to the default as well
     if depth is None:
-        depth = 100 * units.hPa
+        depth = units.Quantity(100, 'hPa')
 
     # Make sure pressure and datavars are the same length
     for datavar in args:
@@ -607,9 +606,11 @@ def get_layer(pressure, *args, height=None, bottom=None, depth=100 * units.hPa,
     if interpolate:
         # If we don't have the bottom or top requested, append them
         if not np.any(np.isclose(top_pressure, p_interp)):
-            p_interp = np.sort(np.append(p_interp.m, top_pressure.m)) * pressure.units
+            p_interp = units.Quantity(np.sort(np.append(p_interp.m, top_pressure.m)),
+                                      pressure.units)
         if not np.any(np.isclose(bottom_pressure, p_interp)):
-            p_interp = np.sort(np.append(p_interp.m, bottom_pressure.m)) * pressure.units
+            p_interp = units.Quantity(np.sort(np.append(p_interp.m, bottom_pressure.m)),
+                                      pressure.units)
 
     ret.append(p_interp[::-1])
 
@@ -847,7 +848,7 @@ def lat_lon_grid_deltas(longitude, latitude, x_dim=-1, y_dim=-2, geod=None):
                               latitude[take_x(slice(1, None))])
     dx[(forward_az < 0.) | (forward_az > 180.)] *= -1
 
-    return dx * units.meter, dy * units.meter
+    return units.Quantity(dx, 'meter'), units.Quantity(dy, 'meter')
 
 
 @exporter.export
@@ -1290,7 +1291,7 @@ def _process_deriv_args(f, axis, x, delta):
             delta_units = getattr(delta, 'units', None)
             delta = np.broadcast_to(delta, diff_size, subok=True)
             if not hasattr(delta, 'units') and delta_units is not None:
-                delta = delta * delta_units
+                delta = units.Quantity(delta, delta_units)
         else:
             delta = _broadcast_to_axis(delta, axis, n)
     elif x is not None:
@@ -1393,13 +1394,9 @@ def angle_to_direction(input_angle, full=False, level=3):
     else:
         scalar = False
 
-    # clean any numeric strings, negatives, and None
-    # does not handle strings with alphabet
-    input_angle = np.array(input_angle).astype(float)
-    with np.errstate(invalid='ignore'):  # warns about the np.nan
-        input_angle[np.where(input_angle < 0)] = np.nan
-
-    input_angle = input_angle * origin_units
+    # clean any numeric strings, negatives, and None does not handle strings with alphabet
+    input_angle = units.Quantity(np.array(input_angle).astype(float), origin_units)
+    input_angle[input_angle < 0] = units.Quantity(np.nan, origin_units)
 
     # normalizer used for angles > 360 degree to normalize between 0 - 360
     normalizer = np.array(input_angle.m / MAX_DEGREE_ANGLE.m, dtype=int)
diff --git a/src/metpy/constants.py b/src/metpy/constants.py
index 8d46ae5ddf2..a4bc64c0bf3 100644
--- a/src/metpy/constants.py
+++ b/src/metpy/constants.py
@@ -78,49 +78,49 @@
 # Export all the variables defined in this block
 with exporter:
     # Earth
-    earth_gravity = g = 9.80665 * units('m / s^2')
-    Re = earth_avg_radius = 6371008.7714 * units('m')
-    G = gravitational_constant = 6.67430e-11 * units('m^3 / kg / s^2')
-    GM = geocentric_gravitational_constant = 3986005e8 * units('m^3 / s^2')
-    omega = earth_avg_angular_vel = 7292115e-11 * units('rad / s')
-    d = earth_sfc_avg_dist_sun = 149597870700. * units('m')
-    S = earth_solar_irradiance = 1360.8 * units('W / m^2')
-    delta = earth_max_declination = 23.45 * units('degrees')
-    earth_orbit_eccentricity = 0.0167 * units('dimensionless')
+    earth_gravity = g = units.Quantity(9.80665, 'm / s^2')
+    Re = earth_avg_radius = units.Quantity(6371008.7714, 'm')
+    G = gravitational_constant = units.Quantity(6.67430e-11, 'm^3 / kg / s^2')
+    GM = geocentric_gravitational_constant = units.Quantity(3986005e8, 'm^3 / s^2')
+    omega = earth_avg_angular_vel = units.Quantity(7292115e-11, 'rad / s')
+    d = earth_sfc_avg_dist_sun = units.Quantity(149597870700., 'm')
+    S = earth_solar_irradiance = units.Quantity(1360.8, 'W / m^2')
+    delta = earth_max_declination = units.Quantity(23.45, 'degrees')
+    earth_orbit_eccentricity = units.Quantity(0.0167, 'dimensionless')
     earth_mass = me = geocentric_gravitational_constant / gravitational_constant
 
     # molar gas constant
-    R = 8.314462618 * units('J / mol / K')
+    R = units.Quantity(8.314462618, 'J / mol / K')
 
     # Water
-    Mw = water_molecular_weight = 18.015268 * units('g / mol')
+    Mw = water_molecular_weight = units.Quantity(18.015268, 'g / mol')
     Rv = water_gas_constant = R / Mw
-    rho_l = density_water = 999.97495 * units('kg / m^3')
-    wv_specific_heat_ratio = 1.330 * units('dimensionless')
+    rho_l = density_water = units.Quantity(999.97495, 'kg / m^3')
+    wv_specific_heat_ratio = units.Quantity(1.330, 'dimensionless')
     Cp_v = wv_specific_heat_press = (
         wv_specific_heat_ratio * Rv / (wv_specific_heat_ratio - 1)
     )
     Cv_v = wv_specific_heat_vol = Cp_v / wv_specific_heat_ratio
-    Cp_l = water_specific_heat = 4.2194 * units('kJ / kg / K')
-    Lv = water_heat_vaporization = 2.50084e6 * units('J / kg')
-    Lf = water_heat_fusion = 3.337e5 * units('J / kg')
-    Cp_i = ice_specific_heat = 2090 * units('J / kg / K')
-    rho_i = density_ice = 917 * units('kg / m^3')
+    Cp_l = water_specific_heat = units.Quantity(4.2194, 'kJ / kg / K')
+    Lv = water_heat_vaporization = units.Quantity(2.50084e6, 'J / kg')
+    Lf = water_heat_fusion = units.Quantity(3.337e5, 'J / kg')
+    Cp_i = ice_specific_heat = units.Quantity(2090, 'J / kg / K')
+    rho_i = density_ice = units.Quantity(917, 'kg / m^3')
 
     # Dry air
-    Md = dry_air_molecular_weight = 28.96546e-3 * units('kg / mol')
+    Md = dry_air_molecular_weight = units.Quantity(28.96546e-3, 'kg / mol')
     Rd = dry_air_gas_constant = R / Md
-    dry_air_spec_heat_ratio = 1.4 * units('dimensionless')
+    dry_air_spec_heat_ratio = units.Quantity(1.4, 'dimensionless')
     Cp_d = dry_air_spec_heat_press = (
         dry_air_spec_heat_ratio * Rd / (dry_air_spec_heat_ratio - 1)
     )
     Cv_d = dry_air_spec_heat_vol = Cp_d / dry_air_spec_heat_ratio
     rho_d = dry_air_density_stp = (
-        (1000. * units('mbar')) / (Rd * 273.15 * units('K'))
+        units.Quantity(1000., 'mbar') / (Rd * units.Quantity(273.15, 'K'))
     ).to('kg / m^3')
 
     # General meteorology constants
-    P0 = pot_temp_ref_press = 1000. * units('mbar')
+    P0 = pot_temp_ref_press = units.Quantity(1000., 'mbar')
     kappa = poisson_exponent = (Rd / Cp_d).to('dimensionless')
     gamma_d = dry_adiabatic_lapse_rate = g / Cp_d
     epsilon = molecular_weight_ratio = (Mw / Md).to('dimensionless')
diff --git a/src/metpy/io/metar.py b/src/metpy/io/metar.py
index 67a60069431..da1680b4a42 100644
--- a/src/metpy/io/metar.py
+++ b/src/metpy/io/metar.py
@@ -158,17 +158,16 @@ def parse_metar_to_dataframe(metar_text, *, year=None, month=None):
     try:
         # Create a field for sea-level pressure and make sure it is a float
         df['air_pressure_at_sea_level'] = float(altimeter_to_sea_level_pressure(
-            df.altimeter.values * units('inHg'),
-            df.elevation.values * units('meters'),
-            df.temperature.values * units('degC')).to('hPa').magnitude)
+            units.Quantity(df.altimeter.values, 'inHg'),
+            units.Quantity(df.elevation.values, 'meters'),
+            units.Quantity(df.temperature.values, 'degC')).to('hPa').magnitude)
     except AttributeError:
         df['air_pressure_at_sea_level'] = [np.nan]
 
     # Use get wind components and assign them to u and v variables
-    df['eastward_wind'], df['northward_wind'] = wind_components((df.wind_speed.values
-                                                                 * units.kts),
-                                                                df.wind_direction.values
-                                                                * units.degree)
+    df['eastward_wind'], df['northward_wind'] = wind_components(
+        units.Quantity(df.wind_speed.values, 'kts'),
+        units.Quantity(df.wind_direction.values, 'degree'))
 
     # Round the altimeter and sea-level pressure values
     df['altimeter'] = df.altimeter.round(2)
@@ -424,7 +423,7 @@ def parse_metar_to_named_tuple(metar_text, station_metadata, year, month):
         if (float(tree.altim.text.strip()[1:5])) > 1100:
             altim = float(tree.altim.text.strip()[1:5]) / 100
         else:
-            altim = (int(tree.altim.text.strip()[1:5]) * units.hPa).to('inHg').magnitude
+            altim = units.Quantity(int(tree.altim.text.strip()[1:5]), 'hPa').to('inHg').m
 
     # Returns a named tuple with all the relevant variables
     return Metar(station_id, lat, lon, elev, date_time, wind_dir, wind_spd,
@@ -617,15 +616,14 @@ def merge(x, key='     '):
 
     # Calculate sea-level pressure from function in metpy.calc
     df['air_pressure_at_sea_level'] = altimeter_to_sea_level_pressure(
-        altim * units('inHg'),
-        elev * units('meters'),
-        temp * units('degC')).to('hPa').magnitude
+        units.Quantity(altim, 'inHg'),
+        units.Quantity(elev, 'meters'),
+        units.Quantity(temp, 'degC')).to('hPa').magnitude
 
     # Use get wind components and assign them to eastward and northward winds
-    df['eastward_wind'], df['northward_wind'] = wind_components((df.wind_speed.values
-                                                                 * units.kts),
-                                                                df.wind_direction.values
-                                                                * units.degree)
+    df['eastward_wind'], df['northward_wind'] = wind_components(
+        units.Quantity(df.wind_speed.values, 'kts'),
+        units.Quantity(df.wind_direction.values, 'degree'))
 
     # Drop duplicate values
     df = df.drop_duplicates(subset=['date_time', 'latitude', 'longitude'], keep='last')
diff --git a/src/metpy/io/station_data.py b/src/metpy/io/station_data.py
index 5d31a1c14ab..a22b2dbf498 100644
--- a/src/metpy/io/station_data.py
+++ b/src/metpy/io/station_data.py
@@ -113,7 +113,8 @@ def _read_airports_file(input_file=None):
     station_map = pd.DataFrame({'id': df.ident.values, 'synop_id': 99999,
                                 'latitude': df.latitude_deg.values,
                                 'longitude': df.longitude_deg.values,
-                                'altitude': ((df.elevation_ft.values * units.ft).to('m')).m,
+                                'altitude': units.Quantity(
+                                    df.elevation_ft.values, 'ft').to('m').m,
                                 'country': df.iso_region.str.split('-', n=1,
                                                                    expand=True)[1].values,
                                 'source': input_file
diff --git a/src/metpy/plots/declarative.py b/src/metpy/plots/declarative.py
index a2b49d41555..686da82bcb5 100644
--- a/src/metpy/plots/declarative.py
+++ b/src/metpy/plots/declarative.py
@@ -1570,7 +1570,8 @@ def _build(self):
             else:
                 field_kwargs['plot_units'] = self.plot_units[0]
             if hasattr(self.data, 'units') and (field_kwargs['plot_units'] is not None):
-                parameter = data[ob_type][subset].values * units(self.data.units[ob_type])
+                parameter = units.Quantity(data[ob_type][subset].values,
+                                           self.data.units[ob_type])
             else:
                 parameter = data[ob_type][subset]
             if field_kwargs['formatter'] is not None:
@@ -1592,10 +1593,10 @@ def _build(self):
             vector_kwargs['color'] = self.vector_field_color
             vector_kwargs['plot_units'] = self.vector_plot_units
             if hasattr(self.data, 'units') and (vector_kwargs['plot_units'] is not None):
-                u = (data[self.vector_field[0]][subset].values
-                     * units(self.data.units[self.vector_field[0]]))
-                v = (data[self.vector_field[1]][subset].values
-                     * units(self.data.units[self.vector_field[1]]))
+                u = units.Quantity(data[self.vector_field[0]][subset].values,
+                                   self.data.units[self.vector_field[0]])
+                v = units.Quantity(data[self.vector_field[1]][subset].values,
+                                   self.data.units[self.vector_field[1]])
             else:
                 vector_kwargs.pop('plot_units')
                 u = data[self.vector_field[0]][subset]
diff --git a/src/metpy/plots/skewt.py b/src/metpy/plots/skewt.py
index d08fbacb45a..304d6dd0e0d 100644
--- a/src/metpy/plots/skewt.py
+++ b/src/metpy/plots/skewt.py
@@ -485,14 +485,15 @@ def plot_dry_adiabats(self, t0=None, pressure=None, **kwargs):
         # Determine set of starting temps if necessary
         if t0 is None:
             xmin, xmax = self.ax.get_xlim()
-            t0 = np.arange(xmin, xmax + 1, 10) * units.degC
+            t0 = units.Quantity(np.arange(xmin, xmax + 1, 10), 'degC')
 
         # Get pressure levels based on ylims if necessary
         if pressure is None:
-            pressure = np.linspace(*self.ax.get_ylim()) * units.mbar
+            pressure = units.Quantity(np.linspace(*self.ax.get_ylim()), 'mbar')
 
         # Assemble into data for plotting
-        t = dry_lapse(pressure, t0[:, np.newaxis], 1000. * units.mbar).to(units.degC)
+        t = dry_lapse(pressure, t0[:, np.newaxis],
+                      units.Quantity(1000., 'mbar')).to(units.degC)
         linedata = [np.vstack((ti.m, pressure.m)).T for ti in t]
 
         # Add to plot
@@ -542,15 +543,16 @@ def plot_moist_adiabats(self, t0=None, pressure=None, **kwargs):
         # Determine set of starting temps if necessary
         if t0 is None:
             xmin, xmax = self.ax.get_xlim()
-            t0 = np.concatenate((np.arange(xmin, 0, 10),
-                                 np.arange(0, xmax + 1, 5))) * units.degC
+            t0 = units.Quantity(np.concatenate((np.arange(xmin, 0, 10),
+                                                np.arange(0, xmax + 1, 5))), 'degC')
 
         # Get pressure levels based on ylims if necessary
         if pressure is None:
-            pressure = np.linspace(*self.ax.get_ylim()) * units.mbar
+            pressure = units.Quantity(np.linspace(*self.ax.get_ylim()), 'mbar')
 
         # Assemble into data for plotting
-        t = moist_lapse(pressure, t0[:, np.newaxis], 1000. * units.mbar).to(units.degC)
+        t = moist_lapse(pressure, t0[:, np.newaxis],
+                        units.Quantity(1000., 'mbar')).to(units.degC)
         linedata = [np.vstack((ti.m, pressure.m)).T for ti in t]
 
         # Add to plot
@@ -600,7 +602,7 @@ def plot_mixing_lines(self, mixing_ratio=None, pressure=None, **kwargs):
 
         # Set pressure range if necessary
         if pressure is None:
-            pressure = np.linspace(600, max(self.ax.get_ylim())) * units.mbar
+            pressure = units.Quantity(np.linspace(600, max(self.ax.get_ylim())), 'mbar')
 
         # Assemble data for plotting
         td = dewpoint(vapor_pressure(pressure, mixing_ratio))
@@ -930,10 +932,10 @@ def plot_colormapped(self, u, v, c, intervals=None, colors=None, **kwargs):
             if intervals.dimensionality == {'[length]': 1.0}:
 
                 # Find any intervals not in the data and interpolate them
-                interpolation_heights = [bound.m for bound in intervals if bound not in c]
-                interpolation_heights = np.array(interpolation_heights) * intervals.units
-                interpolation_heights = (np.sort(interpolation_heights.magnitude)
-                                         * interpolation_heights.units)
+                interpolation_heights = np.array([bound.m for bound in intervals
+                                                  if bound not in c])
+                interpolation_heights = units.Quantity(np.sort(interpolation_heights),
+                                                       intervals.units)
                 (interpolated_heights, interpolated_u,
                  interpolated_v) = interpolate_1d(interpolation_heights, c, c, u, v)
 
@@ -952,7 +954,7 @@ def plot_colormapped(self, u, v, c, intervals=None, colors=None, **kwargs):
                 intervals = intervals.to_base_units()
             # If segmenting by anything else, do not interpolate, just use the data
             else:
-                intervals = np.asarray(intervals) * intervals.units
+                intervals = units.Quantity(np.asarray(intervals), intervals.units)
 
             norm = mcolors.BoundaryNorm(intervals.magnitude, cmap.N)
             cmap.set_over('none')
diff --git a/src/metpy/testing.py b/src/metpy/testing.py
index d1556e35411..ec56a1f2663 100644
--- a/src/metpy/testing.py
+++ b/src/metpy/testing.py
@@ -84,12 +84,12 @@ def to_float(s):
 
     p, z, t, td, direc, spd = np.array(arr_data).T
 
-    p = p * units.hPa
-    z = z * units.meters
-    t = t * units.degC
-    td = td * units.degC
-    direc = direc * units.degrees
-    spd = spd * units.knots
+    p = units.Quantity(p, 'hPa')
+    z = units.Quantity(z, 'meters')
+    t = units.Quantity(t, 'degC')
+    td = units.Quantity(td, 'degC')
+    direc = units.Quantity(direc, 'degrees')
+    spd = units.Quantity(spd, 'knots')
 
     u, v = wind_components(spd, direc)
 
@@ -158,9 +158,9 @@ def check_mask(actual, desired):
     np.testing.assert_array_equal(actual_mask, desired_mask)
 
 
-def assert_nan(value, units):
+def assert_nan(value, value_units):
     """Check for nan with proper units."""
-    value, _ = check_and_drop_units(value, np.nan * units)
+    value, _ = check_and_drop_units(value, units.Quantity(np.nan, value_units))
     assert np.isnan(value)
 
 
diff --git a/src/metpy/units.py b/src/metpy/units.py
index 8e172d18da2..5915bbd6259 100644
--- a/src/metpy/units.py
+++ b/src/metpy/units.py
@@ -92,7 +92,7 @@ def pandas_dataframe_to_unit_arrays(df, column_units=None):
     res = {}
     for column in df:
         if column in column_units and column_units[column]:
-            res[column] = df[column].values * units(column_units[column])
+            res[column] = units.Quantity(df[column].values, column_units[column])
         else:
             res[column] = df[column].values
     return res
diff --git a/src/metpy/xarray.py b/src/metpy/xarray.py
index 41e9dc01a0b..29b6409c917 100644
--- a/src/metpy/xarray.py
+++ b/src/metpy/xarray.py
@@ -444,8 +444,8 @@ def coordinates_identical(self, other):
     @property
     def time_deltas(self):
         """Return the time difference of the data in seconds (to microsecond precision)."""
-        return (np.diff(self._data_array.values).astype('timedelta64[us]').astype('int64')
-                / 1e6 * units.s)
+        us_diffs = np.diff(self._data_array.values).astype('timedelta64[us]').astype('int64')
+        return units.Quantity(us_diffs / 1e6, 's')
 
     def find_axis_name(self, axis):
         """Return the name of the axis corresponding to the given identifier.
@@ -790,7 +790,7 @@ def _rebuild_coords(self, var, crs):
                         new_coord_var = coord_var.copy(
                             data=(
                                 coord_var.metpy.unit_array
-                                * (height * units.meter)
+                                * units.Quantity(height, 'meter')
                             ).m_as('meter')
                         )
                         new_coord_var.attrs['units'] = 'meter'
@@ -1411,10 +1411,12 @@ def grid_deltas_from_dataarray(f, kind='default'):
         dy_units = units(y.attrs.get('units'))
 
     # Broadcast to input and attach units
-    dx = dx_var.set_dims(f.dims, shape=[dx_var.sizes[dim] if dim in dx_var.dims else 1
-                                        for dim in f.dims]).data * dx_units
-    dy = dy_var.set_dims(f.dims, shape=[dy_var.sizes[dim] if dim in dy_var.dims else 1
-                                        for dim in f.dims]).data * dy_units
+    dx_var = dx_var.set_dims(f.dims, shape=[dx_var.sizes[dim] if dim in dx_var.dims else 1
+                                            for dim in f.dims])
+    dx = units.Quantity(dx_var.data, dx_units)
+    dy_var = dy_var.set_dims(f.dims, shape=[dy_var.sizes[dim] if dim in dy_var.dims else 1
+                                            for dim in f.dims])
+    dy = units.Quantity(dy_var.data, dy_units)
 
     return dx, dy
 
diff --git a/tests/calc/test_calc_tools.py b/tests/calc/test_calc_tools.py
index d18346556ab..d2876a51ddc 100644
--- a/tests/calc/test_calc_tools.py
+++ b/tests/calc/test_calc_tools.py
@@ -358,6 +358,17 @@ def test_get_layer(pressure, variable, heights, bottom, depth, interp, expected)
     assert_array_almost_equal(y_layer, expected[1], 3)
 
 
+def test_get_layer_masked():
+    """Test get_layer with masked arrays as input."""
+    p = units.Quantity(np.ma.array([1000, 500, 400]), 'hPa')
+    u = units.Quantity(np.arange(3), 'm/s')
+    p_layer, u_layer = get_layer(p, u, depth=units.Quantity(6000, 'm'))
+    true_p_layer = units.Quantity([1000., 500., 464.4742], 'hPa')
+    true_u_layer = units.Quantity([0., 1., 1.3303], 'm/s')
+    assert_array_almost_equal(p_layer, true_p_layer, 4)
+    assert_array_almost_equal(u_layer, true_u_layer, 4)
+
+
 def test_greater_or_close():
     """Test floating point greater or close to."""
     x = np.array([0.0, 1.0, 1.49999, 1.5, 1.5000, 1.7])