Repository: Unidata/MetPy

Issue:
Improve setup of pint UnitRegistry
There are now functions for `set_application_registry()` and `get_application_registry()`, so we should really look at using those for a registry. This would allow us to play nicely with other libraries/applications using pint, such as xarray-pint. Questions:
* We currently add our own units and loosen the temperature unit rules. What do we do about this?
* Should we modify existing registries?
* Should we ever set the application registry?

Discussion:
dopplershift: @jthielen What do you think of this?
jthielen: Overall, this looks good to me! After testing out various permutations of MetPy (before and after this change), pint-xarray, and cf-xarray, this indeed looks like the best approach we can practically take for compatibility.

One corner case I discovered: if you import cf-xarray's registry (which defines a new registry and sets it with `pint.set_application_registry`), previously defined `Quantity`s start failing some of MetPy's internal `isinstance(thing, units.Quantity)` checks. For example, while contrived, this demonstrates the issue 

```python
import xarray as xr
import pint_xarray

from metpy.units import units

data_pre = xr.DataArray([3, 5, 7] * units('m/s'))

print(data_pre.pint.units, data_pre.metpy.units)

import cf_xarray.units

data_post = xr.DataArray([3, 5, 7] * units('m/s'))

print(data_pre.pint.units, data_pre.metpy.units)
print(data_post.pint.units, data_post.metpy.units)
```

```
meter / second meter / second
meter / second dimensionless
meter / second meter / second
```

This should be fixable by https://github.com/jthielen/MetPy/commit/6f4c090fe32faa5e217605ec0b8e1d9be4cb165f (comparing against pint's top-level `Quantity` rather than our imported registry's `Quantity`). While we definitely can't expect to fully support mixed registries (given that Pint itself usually doesn't), by changing these checks, we can fix errors like "`u` requires "[speed]" but given "dimensionless" that would show up in situations like this and https://github.com/Unidata/MetPy/issues/2672.

Patch:
diff --git a/src/metpy/interpolate/slices.py b/src/metpy/interpolate/slices.py
index 9ac27ed1ff..50bd70f54e 100644
--- a/src/metpy/interpolate/slices.py
+++ b/src/metpy/interpolate/slices.py
@@ -7,7 +7,7 @@ import numpy as np
 import xarray as xr
 
 from ..package_tools import Exporter
-from ..units import units
+from ..units import is_quantity, units
 from ..xarray import check_axis
 
 exporter = Exporter(globals())
@@ -57,10 +57,7 @@ def interpolate_to_slice(data, points, interp_type='linear'):
     data_sliced.coords['index'] = range(len(points))
 
     # Bug in xarray: interp strips units
-    if (
-        isinstance(data.data, units.Quantity)
-        and not isinstance(data_sliced.data, units.Quantity)
-    ):
+    if is_quantity(data.data) and not is_quantity(data_sliced.data):
         data_sliced.data = units.Quantity(data_sliced.data, data.data.units)
 
     return data_sliced
diff --git a/src/metpy/units.py b/src/metpy/units.py
index cddf9161f6..253f2dbbcc 100644
--- a/src/metpy/units.py
+++ b/src/metpy/units.py
@@ -15,6 +15,7 @@ units : :class:`pint.UnitRegistry`
     import this registry and use it to grab units.
 
 """
+import contextlib
 import functools
 from inspect import Parameter, signature
 import logging
@@ -37,38 +38,55 @@ _base_unit_of_dimensionality = {
     '[speed]': 'm s**-1'
 }
 
-# Create registry, with preprocessors for UDUNITS-style powers (m2 s-2) and percent signs
-units = pint.UnitRegistry(
-    autoconvert_offset_to_baseunit=True,
-    preprocessors=[
-        functools.partial(
-            re.sub,
-            r'(?<=[A-Za-z\)])(?![A-Za-z\)])(?<![0-9\-][eE])(?<![0-9\-])(?=[0-9\-])',
-            '**'
-        ),
-        lambda string: string.replace('%', 'percent')
-    ]
-)
-
-# Capture v0.10 NEP 18 warning on first creation
-with warnings.catch_warnings():
-    warnings.simplefilter('ignore')
-    units.Quantity([])
-
-# Add a percent unit
-units.define('percent = 0.01 = %')
-
-# Define commonly encountered units not defined by pint
-units.define('degrees_north = degree = degrees_N = degreesN = degree_north = degree_N '
-             '= degreeN')
-units.define('degrees_east = degree = degrees_E = degreesE = degree_east = degree_E = degreeE')
-
-# Alias geopotential meters (gpm) to just meters
-units.define('@alias meter = gpm')
+
+def _fix_udunits_powers(string):
+    """Replace UDUNITS-style powers (m2 s-2) with exponent symbols (m**2 s**-2)."""
+    return _UDUNIT_POWER.sub('**', string)
+
+
+# Fix UDUNITS-style powers and percent signs
+_UDUNIT_POWER = re.compile(r'(?<=[A-Za-z\)])(?![A-Za-z\)])'
+                           r'(?<![0-9\-][eE])(?<![0-9\-])(?=[0-9\-])')
+_unit_preprocessors = [_fix_udunits_powers, lambda string: string.replace('%', 'percent')]
+
+
+def setup_registry(reg):
+    """Set up a given registry with MetPy's default tweaks and settings."""
+    reg.autoconvert_offset_to_baseunit = True
+
+    # For Pint 0.18.0, need to deal with the fact that the wrapper isn't forwarding on setting
+    # the attribute.
+    with contextlib.suppress(AttributeError):
+        reg.get().autoconvert_offset_to_baseunit = True
+
+    for pre in _unit_preprocessors:
+        if pre not in reg.preprocessors:
+            reg.preprocessors.append(pre)
+
+    # Add a percent unit
+    reg.define('percent = 0.01 = %')
+
+    # Define commonly encountered units not defined by pint
+    reg.define('degrees_north = degree = degrees_N = degreesN = degree_north = degree_N '
+               '= degreeN')
+    reg.define('degrees_east = degree = degrees_E = degreesE = degree_east = degree_E '
+               '= degreeE')
+
+    # Alias geopotential meters (gpm) to just meters
+    reg.define('@alias meter = gpm')
+
+    # Enable pint's built-in matplotlib support
+    reg.setup_matplotlib()
+
+    return reg
+
+
+# Make our modifications using pint's application registry--which allows us to better
+# interoperate with other libraries using Pint.
+units = setup_registry(pint.get_application_registry())
 
 # Silence UnitStrippedWarning
-if hasattr(pint, 'UnitStrippedWarning'):
-    warnings.simplefilter('ignore', category=pint.UnitStrippedWarning)
+warnings.simplefilter('ignore', category=pint.UnitStrippedWarning)
 
 
 def pandas_dataframe_to_unit_arrays(df, column_units=None):
@@ -106,6 +124,11 @@ def pandas_dataframe_to_unit_arrays(df, column_units=None):
     return res
 
 
+def is_quantity(*args):
+    """Check whether an instance is a quantity."""
+    return all(isinstance(a, pint.Quantity) for a in args)
+
+
 def concatenate(arrs, axis=0):
     r"""Concatenate multiple values into a new quantity.
 
@@ -371,9 +394,3 @@ def process_units(
 
         return wrapper
     return dec
-
-
-# Enable pint's built-in matplotlib support
-units.setup_matplotlib()
-
-del pint
diff --git a/src/metpy/xarray.py b/src/metpy/xarray.py
index d9d0cb7edf..598dad5aa8 100644
--- a/src/metpy/xarray.py
+++ b/src/metpy/xarray.py
@@ -28,7 +28,8 @@ from pyproj import CRS, Proj
 import xarray as xr
 
 from ._vendor.xarray import either_dict_or_kwargs, expanded_indexer, is_dict_like
-from .units import _mutate_arguments, DimensionalityError, UndefinedUnitError, units
+from .units import (_mutate_arguments, DimensionalityError, is_quantity, UndefinedUnitError,
+                    units)
 
 __all__ = ('MetPyDataArrayAccessor', 'MetPyDatasetAccessor', 'grid_deltas_from_dataarray')
 metpy_axes = ['time', 'vertical', 'y', 'latitude', 'x', 'longitude']
@@ -130,7 +131,7 @@ class MetPyDataArrayAccessor:
     @property
     def units(self):
         """Return the units of this DataArray as a `pint.Unit`."""
-        if isinstance(self._data_array.variable._data, units.Quantity):
+        if is_quantity(self._data_array.variable._data):
             return self._data_array.variable._data.units
         else:
             return units.parse_units(self._data_array.attrs.get('units', 'dimensionless'))
@@ -138,7 +139,7 @@ class MetPyDataArrayAccessor:
     @property
     def magnitude(self):
         """Return the magnitude of the data values of this DataArray (i.e., without units)."""
-        if isinstance(self._data_array.data, units.Quantity):
+        if is_quantity(self._data_array.data):
             return self._data_array.data.magnitude
         else:
             return self._data_array.data
@@ -153,7 +154,7 @@ class MetPyDataArrayAccessor:
         will be loaded into memory by this operation. Do not utilize on moderate- to
         large-sized remote datasets before subsetting!
         """
-        if isinstance(self._data_array.data, units.Quantity):
+        if is_quantity(self._data_array.data):
             return self._data_array.data
         else:
             return units.Quantity(self._data_array.data, self.units)
@@ -222,7 +223,7 @@ class MetPyDataArrayAccessor:
         subsetting!
         """
         if (
-            not isinstance(self._data_array.data, units.Quantity)
+            not is_quantity(self._data_array.data)
             and np.issubdtype(self._data_array.data.dtype, np.number)
         ):
             # Only quantify if not already quantified and is quantifiable
@@ -235,7 +236,7 @@ class MetPyDataArrayAccessor:
 
     def dequantify(self):
         """Return a new DataArray with the data as magnitude and the units as an attribute."""
-        if isinstance(self._data_array.data, units.Quantity):
+        if is_quantity(self._data_array.data):
             # Only dequantify if quantified
             dequantified_dataarray = self._data_array.copy(
                 data=self._data_array.data.magnitude
@@ -1277,7 +1278,7 @@ def _wrap_output_like_matching_units(result, match):
         return result if output_xarray else result.metpy.unit_array
     else:
         result = (
-            result.to(match_units) if isinstance(result, units.Quantity)
+            result.to(match_units) if is_quantity(result)
             else units.Quantity(result, match_units)
         )
         return (
@@ -1293,10 +1294,8 @@ def _wrap_output_like_not_matching_units(result, match):
         return result if output_xarray else result.metpy.unit_array
     # Determine if need to upcast to Quantity
     if (
-        not isinstance(result, units.Quantity)
-        and (
-            isinstance(match, units.Quantity)
-            or (output_xarray and isinstance(match.data, units.Quantity))
+        not is_quantity(result) and (
+            is_quantity(match) or (output_xarray and is_quantity(match.data))
         )
     ):
         result = units.Quantity(result)


Test patch:
diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py
index 1766b91bb2..37d868b851 100644
--- a/tests/calc/test_thermo.py
+++ b/tests/calc/test_thermo.py
@@ -35,7 +35,7 @@ from metpy.calc import (brunt_vaisala_frequency, brunt_vaisala_frequency_squared
                         virtual_temperature, wet_bulb_temperature)
 from metpy.calc.thermo import _find_append_zero_crossings
 from metpy.testing import assert_almost_equal, assert_array_almost_equal, assert_nan
-from metpy.units import masked_array, units
+from metpy.units import is_quantity, masked_array, units
 
 
 def test_relative_humidity_from_dewpoint():
@@ -632,7 +632,7 @@ def test_equivalent_potential_temperature_masked():
         np.ma.array([311.18586, 313.51781, 315.93971], mask=[False, True, False]),
         units.kelvin
     )
-    assert isinstance(ept, units.Quantity)
+    assert is_quantity(ept)
     assert isinstance(ept.m, np.ma.MaskedArray)
     assert_array_almost_equal(ept, expected, 3)
 
@@ -656,7 +656,7 @@ def test_saturation_equivalent_potential_temperature_masked():
         np.ma.array([335.02750, 338.95813, 343.08740]),
         units.kelvin
     )
-    assert isinstance(s_ept, units.Quantity)
+    assert is_quantity(s_ept)
     assert isinstance(s_ept.m, np.ma.MaskedArray)
     assert_array_almost_equal(s_ept, expected, 3)
 
diff --git a/tests/io/test_metar.py b/tests/io/test_metar.py
index 026a5385ad..341fb27e76 100644
--- a/tests/io/test_metar.py
+++ b/tests/io/test_metar.py
@@ -12,7 +12,7 @@ from metpy.cbook import get_test_data
 from metpy.io import parse_metar_file, parse_metar_to_dataframe
 from metpy.io._metar_parser.metar_parser import parse
 from metpy.io.metar import Metar, parse_metar
-from metpy.units import units
+from metpy.units import is_quantity, units
 
 
 @pytest.mark.parametrize(['metar', 'truth'], [
@@ -365,7 +365,7 @@ def test_parse_no_pint_objects_in_df():
 
     for df in (parse_metar_file(input_file), parse_metar_to_dataframe(metar_str)):
         for column in df:
-            assert not isinstance(df[column][0], units.Quantity)
+            assert not is_quantity(df[column][0])
 
 
 def test_repr():
diff --git a/tests/test_xarray.py b/tests/test_xarray.py
index 0805501739..8f336eab3e 100644
--- a/tests/test_xarray.py
+++ b/tests/test_xarray.py
@@ -13,7 +13,7 @@ import xarray as xr
 from metpy.plots.mapping import CFProjection
 from metpy.testing import (assert_almost_equal, assert_array_almost_equal, assert_array_equal,
                            get_test_data)
-from metpy.units import DimensionalityError, units
+from metpy.units import DimensionalityError, is_quantity, units
 from metpy.xarray import (add_grid_arguments_from_xarray, add_vertical_dim_from_xarray,
                           check_axis, check_matching_coordinates, grid_deltas_from_dataarray,
                           preprocess_and_wrap)
@@ -102,7 +102,7 @@ def test_geodetic(test_var, ccrs):
 def test_unit_array(test_var):
     """Test unit handling through the accessor."""
     arr = test_var.metpy.unit_array
-    assert isinstance(arr, units.Quantity)
+    assert is_quantity(arr)
     assert arr.units == units.kelvin
 
 
@@ -177,7 +177,7 @@ def test_quantify(test_ds_generic):
     """Test quantify method for converting data to Quantity."""
     original = test_ds_generic['test'].values
     result = test_ds_generic['test'].metpy.quantify()
-    assert isinstance(result.data, units.Quantity)
+    assert is_quantity(result.data)
     assert result.data.units == units.kelvin
     assert 'units' not in result.attrs
     np.testing.assert_array_almost_equal(result.data, units.Quantity(original))
@@ -197,7 +197,7 @@ def test_dequantify():
 def test_dataset_quantify(test_ds_generic):
     """Test quantify method for converting data to Quantity on Datasets."""
     result = test_ds_generic.metpy.quantify()
-    assert isinstance(result['test'].data, units.Quantity)
+    assert is_quantity(result['test'].data)
     assert result['test'].data.units == units.kelvin
     assert 'units' not in result['test'].attrs
     np.testing.assert_array_almost_equal(
@@ -1365,7 +1365,7 @@ def test_preprocess_and_wrap_with_variable():
 
     assert isinstance(result_12, xr.DataArray)
     xr.testing.assert_identical(func(data1, data2), expected_12)
-    assert isinstance(result_21, units.Quantity)
+    assert is_quantity(result_21)
     assert_array_equal(func(data2, data1), expected_21)
 
 
diff --git a/tests/units/test_units.py b/tests/units/test_units.py
index 68cac4fa33..ac99827c86 100644
--- a/tests/units/test_units.py
+++ b/tests/units/test_units.py
@@ -9,7 +9,8 @@ import pandas as pd
 import pytest
 
 from metpy.testing import assert_array_almost_equal, assert_array_equal, assert_nan
-from metpy.units import check_units, concatenate, pandas_dataframe_to_unit_arrays, units
+from metpy.units import (check_units, concatenate, is_quantity,
+                         pandas_dataframe_to_unit_arrays, units)
 
 
 def test_concatenate():
@@ -177,6 +178,18 @@ def test_added_degrees_units():
     assert units('degrees_east').to_base_units().units == units.radian
 
 
+def test_is_quantity():
+    """Test is_quantity properly works."""
+    assert is_quantity(1 * units.m)
+    assert not is_quantity(np.array([1]))
+
+
+def test_is_quantity_multiple():
+    """Test is_quantity with multiple inputs."""
+    assert is_quantity(1 * units.m, np.array([4.]) * units.degree)
+    assert not is_quantity(1 * units.second, np.array([5., 2.]))
+
+
 def test_gpm_unit():
     """Test that the gpm unit does alias to meters."""
     x = 1 * units('gpm')