Merge pull request #1325 from jthielen/quantity-in-xarray

Switching to Quantity-in-DataArray/NEP 18 Support over DataArray-with-units-attribute

.gitignore

@@ -2,6 +2,7 @@
 .DS_Store
 *.swp
 .idea/
+.vscode/
 
 *.pyc
 

.travis.yml

@@ -49,7 +49,7 @@ matrix:
   include:
     - python: 3.6
       env:
-        - VERSIONS="numpy==1.13.0 matplotlib==2.1.0 scipy==1.0.0 pint==0.8 xarray==0.12.3 pandas==0.22.0"
+        - VERSIONS="numpy==1.16.0 matplotlib==2.1.0 scipy==1.0.0 pint==0.10.1 xarray==0.13.0 pandas==0.22.0"
         - TASK="coverage"
         - TEST_OUTPUT_CONTROL=""
     - python: "3.8-dev"
@@ -126,6 +126,7 @@ before_script:
 
 script:
   - export TEST_DATA_DIR=${TRAVIS_BUILD_DIR}/staticdata;
+  - export NUMPY_EXPERIMENTAL_ARRAY_FUNCTION=1;
   - if [[ $TASK == "docs" ]]; then
       pushd docs;
       make clean overridecheck html linkcheck O=-W;

docs/installguide.rst

@@ -9,11 +9,11 @@ In general, MetPy tries to support minor versions of dependencies released withi
 years. For Python itself, that means supporting the last two minor releases.
 
 * matplotlib >= 2.1.0
-* numpy >= 1.13.0
+* numpy >= 1.16.0
 * scipy >= 1.0.0
-* pint >= 0.8
+* pint >= 0.10.1
 * pandas >= 0.22.0
-* xarray >= 0.12.3
+* xarray >= 0.13.0
 * traitlets >= 4.3.0
 * pooch >= 0.1
 

examples/Four_Panel_Map.py

@@ -63,9 +63,9 @@ def plot_background(ax):
 
 # Do unit conversions to what we wish to plot
 vort_500 = vort_500 * 1e5
-surface_temp.metpy.convert_units('degF')
-precip_water.metpy.convert_units('inches')
-winds_300.metpy.convert_units('knots')
+surface_temp = surface_temp.metpy.convert_units('degF')
+precip_water = precip_water.metpy.convert_units('inches')
+winds_300 = winds_300.metpy.convert_units('knots')
 
 ###########################################
 

examples/cross_section.py

@@ -66,8 +66,8 @@
                                           dims=specific_humidity.dims,
                                           attrs={'units': rh.units})
 
-cross['u_wind'].metpy.convert_units('knots')
-cross['v_wind'].metpy.convert_units('knots')
+cross['u_wind'] = cross['u_wind'].metpy.convert_units('knots')
+cross['v_wind'] = cross['v_wind'].metpy.convert_units('knots')
 cross['t_wind'], cross['n_wind'] = mpcalc.cross_section_components(cross['u_wind'],
                                                                    cross['v_wind'])
 

setup.cfg

@@ -39,10 +39,10 @@ setup_requires = setuptools_scm
 python_requires = >=3.6
 install_requires =
     matplotlib>=2.1.0
-    numpy>=1.13.0
+    numpy>=1.16.0
     scipy>=1.0
-    pint>=0.8
-    xarray>=0.12.3
+    pint>=0.10.1
+    xarray>=0.13.0
     pooch>=0.1
     traitlets>=4.3.0
     pandas>=0.22.0

src/metpy/calc/basic.py

@@ -18,7 +18,7 @@
 from .tools import wrap_output_like
 from .. import constants as mpconsts
 from ..package_tools import Exporter
-from ..units import atleast_1d, check_units, masked_array, units
+from ..units import check_units, masked_array, units
 from ..xarray import preprocess_xarray
 
 exporter = Exporter(globals())
@@ -90,7 +90,7 @@ def wind_direction(u, v, convention='from'):
     """
     wdir = 90. * units.deg - np.arctan2(-v, -u)
     origshape = wdir.shape
-    wdir = atleast_1d(wdir)
+    wdir = np.atleast_1d(wdir)
 
     # Handle oceanographic convection
     if convention == 'to':
@@ -244,8 +244,8 @@ def heat_index(temperature, relative_humidity, mask_undefined=True):
     windchill
 
     """
-    temperature = atleast_1d(temperature)
-    relative_humidity = atleast_1d(relative_humidity)
+    temperature = np.atleast_1d(temperature)
+    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
@@ -357,9 +357,9 @@ def apparent_temperature(temperature, relative_humidity, speed, face_level_winds
     """
     is_not_scalar = isinstance(temperature.m, (list, tuple, np.ndarray))
 
-    temperature = atleast_1d(temperature)
-    relative_humidity = atleast_1d(relative_humidity)
-    speed = atleast_1d(speed)
+    temperature = np.atleast_1d(temperature)
+    relative_humidity = np.atleast_1d(relative_humidity)
+    speed = np.atleast_1d(speed)
 
     # NB: mask_defined=True is needed to know where computed values exist
     wind_chill_temperature = windchill(temperature, speed, face_level_winds=face_level_winds,
@@ -388,7 +388,7 @@ def apparent_temperature(temperature, relative_humidity, speed, face_level_winds
     if is_not_scalar:
         return app_temperature
     else:
-        return atleast_1d(app_temperature)[0]
+        return np.atleast_1d(app_temperature)[0]
 
 
 @exporter.export

src/metpy/calc/cross_sections.py

@@ -13,6 +13,7 @@
 from .basic import coriolis_parameter
 from .tools import first_derivative
 from ..package_tools import Exporter
+from ..units import units
 from ..xarray import check_axis, check_matching_coordinates
 
 exporter = Exporter(globals())
@@ -49,8 +50,8 @@ def distances_from_cross_section(cross):
         y = distance * np.cos(np.deg2rad(forward_az))
 
         # Build into DataArrays
-        x = xr.DataArray(x, coords=lon.coords, dims=lon.dims, attrs={'units': 'meters'})
-        y = xr.DataArray(y, coords=lat.coords, dims=lat.dims, attrs={'units': 'meters'})
+        x = xr.DataArray(x * units.meter, coords=lon.coords, dims=lon.dims)
+        y = xr.DataArray(y * units.meter, coords=lat.coords, dims=lat.dims)
 
     elif check_axis(cross.metpy.x, 'x') and check_axis(cross.metpy.y, 'y'):
 
@@ -86,8 +87,7 @@ def latitude_from_cross_section(cross):
         latitude = ccrs.Geodetic().transform_points(cross.metpy.cartopy_crs,
                                                     cross.metpy.x.values,
                                                     y.values)[..., 1]
-        latitude = xr.DataArray(latitude, coords=y.coords, dims=y.dims,
-                                attrs={'units': 'degrees_north'})
+        latitude = xr.DataArray(latitude * units.degrees_north, coords=y.coords, dims=y.dims)
         return latitude
 
 
@@ -165,10 +165,6 @@ def cross_section_components(data_x, data_y, index='index'):
     component_tang = data_x * unit_tang[0] + data_y * unit_tang[1]
     component_norm = data_x * unit_norm[0] + data_y * unit_norm[1]
 
-    # Reattach units (only reliable attribute after operation)
-    component_tang.attrs = {'units': data_x.attrs['units']}
-    component_norm.attrs = {'units': data_x.attrs['units']}
-
     return component_tang, component_norm
 
 
@@ -207,9 +203,8 @@ def normal_component(data_x, data_y, index='index'):
     component_norm = data_x * unit_norm[0] + data_y * unit_norm[1]
 
     # Reattach only reliable attributes after operation
-    for attr in ('units', 'grid_mapping'):
-        if attr in data_x.attrs:
-            component_norm.attrs[attr] = data_x.attrs[attr]
+    if 'grid_mapping' in data_x.attrs:
+        component_norm.attrs['grid_mapping'] = data_x.attrs['grid_mapping']
 
     return component_norm
 
@@ -249,9 +244,8 @@ def tangential_component(data_x, data_y, index='index'):
     component_tang = data_x * unit_tang[0] + data_y * unit_tang[1]
 
     # Reattach only reliable attributes after operation
-    for attr in ('units', 'grid_mapping'):
-        if attr in data_x.attrs:
-            component_tang.attrs[attr] = data_x.attrs[attr]
+    if 'grid_mapping' in data_x.attrs:
+        component_tang.attrs['grid_mapping'] = data_x.attrs['grid_mapping']
 
     return component_tang
 
@@ -292,20 +286,16 @@ def absolute_momentum(u, v, index='index'):
 
     """
     # Get the normal component of the wind
-    norm_wind = normal_component(u, v, index=index)
-    norm_wind.metpy.convert_units('m/s')
+    norm_wind = normal_component(u, v, index=index).metpy.convert_units('m/s')
 
     # Get other pieces of calculation (all as ndarrays matching shape of norm_wind)
-    latitude = latitude_from_cross_section(norm_wind)  # in degrees_north
+    latitude = latitude_from_cross_section(norm_wind)
     _, latitude = xr.broadcast(norm_wind, latitude)
-    f = coriolis_parameter(np.deg2rad(latitude.values)).magnitude  # in 1/s
+    f = coriolis_parameter(np.deg2rad(latitude.values))
     x, y = distances_from_cross_section(norm_wind)
-    x.metpy.convert_units('meters')
-    y.metpy.convert_units('meters')
+    x = x.metpy.convert_units('meters')
+    y = y.metpy.convert_units('meters')
     _, x, y = xr.broadcast(norm_wind, x, y)
-    distance = np.hypot(x, y).values  # in meters
-
-    m = norm_wind + f * distance
-    m.attrs = {'units': norm_wind.attrs['units']}
+    distance = np.hypot(x.metpy.quantify(), y.metpy.quantify())
 
-    return m
+    return norm_wind + f * distance

src/metpy/calc/indices.py

@@ -8,7 +8,7 @@
 from .tools import _remove_nans, get_layer
 from .. import constants as mpconsts
 from ..package_tools import Exporter
-from ..units import atleast_1d, check_units, concatenate, units
+from ..units import check_units, concatenate, units
 from ..xarray import preprocess_xarray
 
 exporter = Exporter(globals())
@@ -258,7 +258,7 @@ def supercell_composite(mucape, effective_storm_helicity, effective_shear):
         supercell composite
 
     """
-    effective_shear = np.clip(atleast_1d(effective_shear), None, 20 * units('m/s'))
+    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'))
 
@@ -306,12 +306,12 @@ def significant_tornado(sbcape, surface_based_lcl_height, storm_helicity_1km, sh
         significant tornado parameter
 
     """
-    surface_based_lcl_height = np.clip(atleast_1d(surface_based_lcl_height),
+    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(atleast_1d(shear_6km), None, 30 * units('m/s'))
+    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')
 

src/metpy/calc/kinematics.py

@@ -9,7 +9,7 @@
 from .. import constants as mpconsts
 from ..cbook import iterable
 from ..package_tools import Exporter
-from ..units import atleast_2d, check_units, concatenate, units
+from ..units import check_units, concatenate, units
 from ..xarray import preprocess_xarray
 
 exporter = Exporter(globals())
@@ -260,7 +260,7 @@ def advection(scalar, wind, deltas):
 
     # Make them be at least 2D (handling the 1D case) so that we can do the
     # multiply and sum below
-    grad, wind = atleast_2d(grad, wind)
+    grad, wind = np.atleast_2d(grad, wind)
 
     return (-grad * wind).sum(axis=0)
 

src/metpy/calc/thermo.py

@@ -14,7 +14,7 @@
 from ..cbook import broadcast_indices
 from ..interpolate.one_dimension import interpolate_1d
 from ..package_tools import Exporter
-from ..units import atleast_1d, check_units, concatenate, units
+from ..units import check_units, concatenate, units
 from ..xarray import preprocess_xarray
 
 exporter = Exporter(globals())
@@ -267,7 +267,7 @@ def dt(t, p):
 
     pressure = pressure.to('mbar')
     reference_pressure = reference_pressure.to('mbar')
-    temperature = atleast_1d(temperature)
+    temperature = np.atleast_1d(temperature)
 
     side = 'left'
 
@@ -2374,9 +2374,9 @@ def wet_bulb_temperature(pressure, temperature, dewpoint):
 
     """
     if not hasattr(pressure, 'shape'):
-        pressure = atleast_1d(pressure)
-        temperature = atleast_1d(temperature)
-        dewpoint = atleast_1d(dewpoint)
+        pressure = np.atleast_1d(pressure)
+        temperature = np.atleast_1d(temperature)
+        dewpoint = np.atleast_1d(dewpoint)
 
     it = np.nditer([pressure, temperature, dewpoint, None],
                    op_dtypes=['float', 'float', 'float', 'float'],

src/metpy/calc/tools.py

@@ -16,7 +16,7 @@
 from ..cbook import broadcast_indices, result_type
 from ..interpolate import interpolate_1d, log_interpolate_1d
 from ..package_tools import Exporter
-from ..units import atleast_1d, check_units, concatenate, diff, units
+from ..units import check_units, concatenate, units
 from ..xarray import check_axis, preprocess_xarray
 
 exporter = Exporter(globals())
@@ -1276,7 +1276,7 @@ def _process_deriv_args(f, kwargs):
         if 'x' in kwargs:
             raise ValueError('Cannot specify both "x" and "delta".')
 
-        delta = atleast_1d(kwargs['delta'])
+        delta = np.atleast_1d(kwargs['delta'])
         if delta.size == 1:
             diff_size = list(f.shape)
             diff_size[axis] -= 1
@@ -1288,7 +1288,7 @@ def _process_deriv_args(f, kwargs):
             delta = _broadcast_to_axis(delta, axis, n)
     elif 'x' in kwargs:
         x = _broadcast_to_axis(kwargs['x'], axis, n)
-        delta = diff(x, axis=axis)
+        delta = np.diff(x, axis=axis)
     else:
         raise ValueError('Must specify either "x" or "delta" for value positions.')
 
@@ -1492,7 +1492,7 @@ def wrap_output_like(**wrap_kwargs):
     - As matched output (final returned value):
 
         * ``pint.Quantity``
-        * ``xarray.DataArray``
+        * ``xarray.DataArray`` wrapping a ``pint.Quantity``
 
     (if matched output is not one of these types, we instead treat the match as if it was a
     dimenionless Quantity.)
@@ -1572,16 +1572,27 @@ def wrapper(*args, **kwargs):
 
 def _wrap_output_like_matching_units(result, match):
     """Convert result to be like match with matching units for output wrapper."""
-    match_units = str(getattr(match, 'units', ''))
-    output_xarray = isinstance(match, xr.DataArray)
+    if isinstance(match, xr.DataArray):
+        output_xarray = True
+        match_units = match.metpy.units
+    elif isinstance(match, units.Quantity):
+        output_xarray = False
+        match_units = match.units
+    else:
+        output_xarray = False
+        match_units = ''
+
     if isinstance(result, xr.DataArray):
-        result.metpy.convert_units(match_units)
-        return result if output_xarray else result.metpy.unit_array
+        result = result.metpy.convert_units(match_units)
+        return result.metpy.quantify() if output_xarray else result.metpy.unit_array
     else:
         result = result.m_as(match_units) if isinstance(result, units.Quantity) else result
         if output_xarray:
-            return xr.DataArray(result, dims=match.dims, coords=match.coords,
-                                attrs={'units': match_units})
+            return xr.DataArray(
+                units.Quantity(result, match_units),
+                dims=match.dims,
+                coords=match.coords
+            )
         else:
             return units.Quantity(result, match_units)
 
@@ -1590,7 +1601,7 @@ def _wrap_output_like_not_matching_units(result, match):
     """Convert result to be like match without matching units for output wrapper."""
     output_xarray = isinstance(match, xr.DataArray)
     if isinstance(result, xr.DataArray):
-        return result if output_xarray else result.metpy.unit_array
+        return result.metpy.quantify() if output_xarray else result.metpy.unit_array
     else:
         if isinstance(result, units.Quantity):
             result_magnitude = result.magnitude
@@ -1600,7 +1611,10 @@ def _wrap_output_like_not_matching_units(result, match):
             result_units = ''
 
         if output_xarray:
-            return xr.DataArray(result_magnitude, dims=match.dims, coords=match.coords,
-                                attrs={'units': result_units})
+            return xr.DataArray(
+                units.Quantity(result_magnitude, result_units),
+                dims=match.dims,
+                coords=match.coords
+            )
         else:
             return units.Quantity(result_magnitude, result_units)