Fix up internal multiplication by units

setup.cfg

@@ -71,7 +71,7 @@ max-line-length = 95
 
 [flake8]
 max-line-length = 95
-application-import-names = metpy
+application-import-names = metpy flake8_metpy
 import-order-style = google
 copyright-check = True
 copyright-author = MetPy Developers
@@ -83,15 +83,20 @@ docstring-convention = numpy
 exclude = docs build src/metpy/io/metar_parser.py
 select = A B C D E F H I M Q RST S T W B902
 ignore = F405 W503 RST902 SIM106
-per-file-ignores = examples/*.py: D T003 T001
-                   tutorials/*.py: D T003 T001
+per-file-ignores = examples/*.py: D MPY001 T003 T001
+                   tutorials/*.py: D MPY001 T003 T001
                    src/metpy/xarray.py: RST304
                    src/metpy/deprecation.py: C801
                    src/metpy/calc/*.py: RST306
                    src/metpy/interpolate/*.py: RST306
                    src/metpy/future.py: RST307
-		   src/metpy/constants.py: RST306
+                   src/metpy/constants.py: RST306
                    docs/doc-server.py: T001
+                   tests/*.py: MPY001
+
+[flake8:local-plugins]
+extension = MPY = flake8_metpy:MetPyChecker
+paths = ./tools/flake8_metpy
 
 [tool:pytest]
 # https://github.com/matplotlib/pytest-mpl/issues/69

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