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

@@ -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

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

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')

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,17 +795,16 @@ 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)
     dudp = first_derivative(u, x=pressure, axis=vertical_dim)
     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