Repository: Unidata/MetPy

Issue:
Support units in `inverse_distance_to_points`
Currently `inverse_distance_to_points` chokes when given `values` having units because it does:

https://github.com/Unidata/MetPy/blob/749cf07ec2d15e644f822e45796bd5b20dc3dd7c/src/metpy/interpolate/points.py#L257-L258

This means when calling `interpolate_to_grid` using either Cressman or Barnes, you get a nasty traceback:
```pytb
---------------------------------------------------------------------------
DimensionalityError                       Traceback (most recent call last)
~/miniconda3/envs/py39/lib/python3.9/site-packages/pint/quantity.py in __float__(self)
    896             return float(self._convert_magnitude_not_inplace(UnitsContainer()))
--> 897         raise DimensionalityError(self._units, "dimensionless")
    898 

DimensionalityError: Cannot convert from 'kelvin' to 'dimensionless'

The above exception was the direct cause of the following exception:

ValueError                                Traceback (most recent call last)
/var/folders/r8/bhrhn4m56lgdckj491rm6_w40000gp/T/ipykernel_52617/3693995647.py in <module>
      1 grdx = 25000
----> 2 thgridx,thgridy,theta = interpolate_to_grid(x_masked,y_masked, potemp,interp_type='barnes', kappa_star=10, gamma=0.5,hres=grdx)
      3 
      4 print(np.shape(thgridx))
      5 print(np.unique(theta))

~/repos/metpy/src/metpy/pandas.py in wrapper(*args, **kwargs)
     19         kwargs = {name: (v.values if isinstance(v, pd.Series) else v)
     20                   for name, v in kwargs.items()}
---> 21         return func(*args, **kwargs)
     22     return wrapper

~/repos/metpy/src/metpy/interpolate/grid.py in interpolate_to_grid(x, y, z, interp_type, hres, minimum_neighbors, gamma, kappa_star, search_radius, rbf_func, rbf_smooth, boundary_coords)
    298     points_obs = np.array(list(zip(x, y)))
    299     points_grid = generate_grid_coords(grid_x, grid_y)
--> 300     img = interpolate_to_points(points_obs, z, points_grid, interp_type=interp_type,
    301                                 minimum_neighbors=minimum_neighbors, gamma=gamma,
    302                                 kappa_star=kappa_star, search_radius=search_radius,

~/repos/metpy/src/metpy/interpolate/points.py in interpolate_to_points(points, values, xi, interp_type, minimum_neighbors, gamma, kappa_star, search_radius, rbf_func, rbf_smooth)
    363         else:
    364             kappa = tools.calc_kappa(ave_spacing, kappa_star)
--> 365             return inverse_distance_to_points(points, values, xi, search_radius, gamma, kappa,
    366                                               min_neighbors=minimum_neighbors,
    367                                               kind=interp_type)

~/repos/metpy/src/metpy/interpolate/points.py in inverse_distance_to_points(points, values, xi, r, gamma, kappa, min_neighbors, kind)
    268                 img[idx] = cressman_point(dists, values_subset, r)
    269             elif kind == 'barnes':
--> 270                 img[idx] = barnes_point(dists, values_subset, kappa, gamma)
    271 
    272             else:

ValueError: setting an array element with a sequence.
```
as is seen in [this Stack Overflow question](https://stackoverflow.com/questions/70248893/interpolate-to-grid-returns-all-nans).

This is a pretty bad incompatibility in our API.

Discussion:

Patch:
diff --git a/src/metpy/interpolate/points.py b/src/metpy/interpolate/points.py
index aab9f9c5d7..a68c7d4b7b 100644
--- a/src/metpy/interpolate/points.py
+++ b/src/metpy/interpolate/points.py
@@ -3,6 +3,7 @@
 # SPDX-License-Identifier: BSD-3-Clause
 """Interpolate data valid at one set of points to another in multiple dimensions."""
 
+import functools
 import logging
 
 import numpy as np
@@ -11,6 +12,7 @@ from scipy.spatial import cKDTree, ConvexHull, Delaunay, qhull
 
 from . import geometry, tools
 from ..package_tools import Exporter
+from ..units import units
 
 exporter = Exporter(globals())
 
@@ -193,19 +195,21 @@ def natural_neighbor_to_points(points, values, xi):
 
     """
     tri = Delaunay(points)
-
     members, circumcenters = geometry.find_natural_neighbors(tri, xi)
 
-    img = np.empty(shape=(xi.shape[0]), dtype=values.dtype)
-    img.fill(np.nan)
-
-    for ind, (grid, neighbors) in enumerate(members.items()):
+    if hasattr(values, 'units'):
+        org_units = values.units
+        values = values.magnitude
+    else:
+        org_units = None
 
-        if len(neighbors) > 0:
+    img = np.asarray([natural_neighbor_point(*np.array(points).transpose(),
+                                             values, xi[grid], tri, neighbors, circumcenters)
+                      if len(neighbors) > 0 else np.nan
+                      for grid, neighbors in members.items()])
 
-            points_transposed = np.array(points).transpose()
-            img[ind] = natural_neighbor_point(points_transposed[0], points_transposed[1],
-                                              values, xi[grid], tri, neighbors, circumcenters)
+    if org_units:
+        img = units.Quantity(img, org_units)
 
     return img
 
@@ -216,7 +220,7 @@ def inverse_distance_to_points(points, values, xi, r, gamma=None, kappa=None, mi
     r"""Generate an inverse distance weighting interpolation to the given points.
 
     Values are assigned to the given interpolation points based on either [Cressman1959]_ or
-    [Barnes1964]_. The Barnes implementation used here based on [Koch1983]_.
+    [Barnes1964]_. The Barnes implementation used here is based on [Koch1983]_.
 
     Parameters
     ----------
@@ -227,8 +231,7 @@ def inverse_distance_to_points(points, values, xi, r, gamma=None, kappa=None, mi
     xi: array_like, shape (M, 2)
         Points to interpolate the data onto.
     r: float
-        Radius from grid center, within which observations
-        are considered and weighted.
+        Radius from grid center, within which observations are considered and weighted.
     gamma: float
         Adjustable smoothing parameter for the barnes interpolation. Default None.
     kappa: float
@@ -250,27 +253,28 @@ def inverse_distance_to_points(points, values, xi, r, gamma=None, kappa=None, mi
     inverse_distance_to_grid
 
     """
-    obs_tree = cKDTree(points)
+    if kind == 'cressman':
+        interp_func = functools.partial(cressman_point, radius=r)
+    elif kind == 'barnes':
+        interp_func = functools.partial(barnes_point, kappa=kappa, gamma=gamma)
+    else:
+        raise ValueError(f'{kind} interpolation not supported.')
 
+    obs_tree = cKDTree(points)
     indices = obs_tree.query_ball_point(xi, r=r)
 
-    img = np.empty(shape=(xi.shape[0]), dtype=values.dtype)
-    img.fill(np.nan)
-
-    for idx, (matches, grid) in enumerate(zip(indices, xi)):
-        if len(matches) >= min_neighbors:
-
-            x1, y1 = obs_tree.data[matches].T
-            values_subset = values[matches]
-            dists = geometry.dist_2(grid[0], grid[1], x1, y1)
+    if hasattr(values, 'units'):
+        org_units = values.units
+        values = values.magnitude
+    else:
+        org_units = None
 
-            if kind == 'cressman':
-                img[idx] = cressman_point(dists, values_subset, r)
-            elif kind == 'barnes':
-                img[idx] = barnes_point(dists, values_subset, kappa, gamma)
+    img = np.asarray([interp_func(geometry.dist_2(*grid, *obs_tree.data[matches].T),
+                                  values[matches]) if len(matches) >= min_neighbors else np.nan
+                      for matches, grid in zip(indices, xi)])
 
-            else:
-                raise ValueError(f'{kind} interpolation not supported.')
+    if org_units:
+        img = units.Quantity(img, org_units)
 
     return img
 
@@ -312,7 +316,7 @@ def interpolate_to_points(points, values, xi, interp_type='linear', minimum_neig
     gamma: float
         Adjustable smoothing parameter for the barnes interpolation. Default 0.25.
     kappa_star: float
-        Response parameter for barnes interpolation, specified nondimensionally
+        Response parameter for barnes interpolation, specified non-dimensionally
         in terms of the Nyquist. Default 5.052
     search_radius: float
         A search radius to use for the Barnes and Cressman interpolation schemes.
@@ -343,13 +347,24 @@ def interpolate_to_points(points, values, xi, interp_type='linear', minimum_neig
     """
     # If this is a type that `griddata` handles, hand it along to `griddata`
     if interp_type in ['linear', 'nearest', 'cubic']:
-        return griddata(points, values, xi, method=interp_type)
+        if hasattr(values, 'units'):
+            org_units = values.units
+            values = values.magnitude
+        else:
+            org_units = None
+
+        ret = griddata(points, values, xi, method=interp_type)
+
+        if org_units:
+            ret = units.Quantity(ret, org_units)
+
+        return ret
 
     # If this is natural neighbor, hand it along to `natural_neighbor`
     elif interp_type == 'natural_neighbor':
         return natural_neighbor_to_points(points, values, xi)
 
-    # If this is Barnes/Cressman, determine search_radios and hand it along to
+    # If this is Barnes/Cressman, determine search_radius and hand it along to
     # `inverse_distance`
     elif interp_type in ['cressman', 'barnes']:
         ave_spacing = tools.average_spacing(points)
@@ -368,14 +383,25 @@ def interpolate_to_points(points, values, xi, interp_type='linear', minimum_neig
 
     # If this is radial basis function, make the interpolator and apply it
     elif interp_type == 'rbf':
-
         points_transposed = np.array(points).transpose()
         xi_transposed = np.array(xi).transpose()
+
+        if hasattr(values, 'units'):
+            org_units = values.units
+            values = values.magnitude
+        else:
+            org_units = None
+
         rbfi = Rbf(points_transposed[0], points_transposed[1], values, function=rbf_func,
                    smooth=rbf_smooth)
-        return rbfi(xi_transposed[0], xi_transposed[1])
+        ret = rbfi(xi_transposed[0], xi_transposed[1])
+
+        if org_units:
+            ret = units.Quantity(ret, org_units)
+
+        return ret
 
     else:
-        raise ValueError('Interpolation option not available. '
+        raise ValueError(f'Interpolation option {interp_type} not available. '
                          'Try: linear, nearest, cubic, natural_neighbor, '
                          'barnes, cressman, rbf')
diff --git a/src/metpy/io/metar.py b/src/metpy/io/metar.py
index 6c19248e06..53fea540cb 100644
--- a/src/metpy/io/metar.py
+++ b/src/metpy/io/metar.py
@@ -11,7 +11,7 @@ import warnings
 import numpy as np
 import pandas as pd
 
-from ._metar_parser.metar_parser import parse, ParseError, TreeNode
+from ._metar_parser.metar_parser import parse, ParseError
 from ._tools import open_as_needed
 from .station_data import station_info
 from ..package_tools import Exporter
@@ -61,23 +61,6 @@ col_units = {'station_id': None,
              'current_wx3_symbol': None}
 
 
-def _tree_repr_(self):
-    """Produce string representation of a TreeNodex object."""
-    rep = self.__class__.__name__ + '('
-    args = []
-    for key, value in self.__dict__.items():
-        if key == 'elements':
-            continue
-        args.append(key + '=' + repr(value))
-        dict_str = ', '.join(args) + ')'
-    rep += dict_str
-    return rep
-
-
-# Monkey patch to improve debugging
-TreeNode.__repr__ = _tree_repr_
-
-
 @exporter.export
 def parse_metar_to_dataframe(metar_text, *, year=None, month=None):
     """Parse a single METAR report into a Pandas DataFrame.


Test patch:
diff --git a/tests/interpolate/test_grid.py b/tests/interpolate/test_grid.py
index f59fa91542..a2a621a7fe 100644
--- a/tests/interpolate/test_grid.py
+++ b/tests/interpolate/test_grid.py
@@ -7,7 +7,6 @@
 import logging
 
 import numpy as np
-from numpy.testing import assert_array_almost_equal
 import pytest
 
 from metpy.cbook import get_test_data
@@ -15,6 +14,8 @@ from metpy.interpolate.grid import (generate_grid, generate_grid_coords, get_bou
                                     get_xy_range, get_xy_steps, interpolate_to_grid,
                                     interpolate_to_isosurface, inverse_distance_to_grid,
                                     natural_neighbor_to_grid)
+from metpy.testing import assert_array_almost_equal
+from metpy.units import units
 
 logging.getLogger('metpy.interpolate.grid').setLevel(logging.ERROR)
 
@@ -233,9 +234,10 @@ def test_interpolate_to_isosurface():
     assert_array_almost_equal(truth, dt_theta)
 
 
+@pytest.mark.parametrize('assume_units', [None, 'mbar'])
 @pytest.mark.parametrize('method', interp_methods)
 @pytest.mark.parametrize('boundary_coords', boundary_types)
-def test_interpolate_to_grid(method, test_coords, boundary_coords):
+def test_interpolate_to_grid(method, assume_units, test_coords, boundary_coords):
     r"""Test main grid interpolation function."""
     xp, yp = test_coords
 
@@ -257,11 +259,14 @@ def test_interpolate_to_grid(method, test_coords, boundary_coords):
     if boundary_coords is not None:
         extra_kw['boundary_coords'] = boundary_coords
 
-    _, _, img = interpolate_to_grid(xp, yp, z, hres=10, interp_type=method, **extra_kw)
-
     with get_test_data(f'{method}_test.npz') as fobj:
         truth = np.load(fobj)['img']
 
+    if assume_units:
+        z = units.Quantity(z, assume_units)
+        truth = units.Quantity(truth, assume_units)
+
+    _, _, img = interpolate_to_grid(xp, yp, z, hres=10, interp_type=method, **extra_kw)
     assert_array_almost_equal(truth, img)
 
 
diff --git a/tests/interpolate/test_points.py b/tests/interpolate/test_points.py
index 73bc566c04..beb0ac4906 100644
--- a/tests/interpolate/test_points.py
+++ b/tests/interpolate/test_points.py
@@ -7,7 +7,6 @@
 import logging
 
 import numpy as np
-from numpy.testing import assert_almost_equal, assert_array_almost_equal
 import pytest
 from scipy.spatial import cKDTree, Delaunay
 
@@ -16,6 +15,8 @@ from metpy.interpolate import (interpolate_to_points, inverse_distance_to_points
                                natural_neighbor_to_points)
 from metpy.interpolate.geometry import dist_2, find_natural_neighbors
 from metpy.interpolate.points import barnes_point, cressman_point, natural_neighbor_point
+from metpy.testing import assert_almost_equal, assert_array_almost_equal
+from metpy.units import units
 
 logging.getLogger('metpy.interpolate.points').setLevel(logging.ERROR)
 
@@ -51,7 +52,7 @@ def test_nn_point(test_data):
     members, tri_info = find_natural_neighbors(tri,
                                                list(zip(sim_gridx, sim_gridy)))
 
-    val = natural_neighbor_point(xp, yp, z, [sim_gridx[0], sim_gridy[0]],
+    val = natural_neighbor_point(xp, yp, z, (sim_gridx[0], sim_gridy[0]),
                                  tri, members[0], tri_info)
 
     truth = 1.009
@@ -108,43 +109,51 @@ def test_natural_neighbor_to_points(test_data, test_points):
     assert_array_almost_equal(truth, img)
 
 
-interp_methods = ['cressman', 'barnes', 'shouldraise']
+def test_inverse_distance_to_points_invalid(test_data, test_points):
+    """Test that inverse_distance_to_points raises when given an invalid method."""
+    xp, yp, z = test_data
+    obs_points = np.vstack([xp, yp]).transpose()
+    with pytest.raises(ValueError):
+        inverse_distance_to_points(obs_points, z, test_points, kind='shouldraise', r=40)
 
 
-@pytest.mark.parametrize('method', interp_methods)
-def test_inverse_distance_to_points(method, test_data, test_points):
-    r"""Test inverse distance interpolation to grid function."""
+@pytest.mark.parametrize('assume_units', [None, 'mbar'])
+@pytest.mark.parametrize('method', ['cressman', 'barnes'])
+def test_inverse_distance_to_points(method, assume_units, test_data, test_points):
+    r"""Test inverse distance interpolation to points function."""
     xp, yp, z = test_data
     obs_points = np.vstack([xp, yp]).transpose()
 
-    extra_kw = {}
-    if method == 'shouldraise':
-        extra_kw['r'] = 40
-        with pytest.raises(ValueError):
-            inverse_distance_to_points(
-                obs_points, z, test_points, kind=method, **extra_kw)
-    else:
-        test_file = ''
-        if method == 'cressman':
-            extra_kw['r'] = 20
-            extra_kw['min_neighbors'] = 1
-            test_file = 'cressman_r20_mn1.npz'
-        elif method == 'barnes':
-            extra_kw['r'] = 40
-            extra_kw['kappa'] = 100
-            test_file = 'barnes_r40_k100.npz'
+    extra_kw, test_file = {'cressman': ({'r': 20, 'min_neighbors': 1}, 'cressman_r20_mn1.npz'),
+                           'barnes': ({'r': 40, 'kappa': 100}, 'barnes_r40_k100.npz')}[method]
+
+    with get_test_data(test_file) as fobj:
+        truth = np.load(fobj)['img'].reshape(-1)
+
+    if assume_units:
+        z = units.Quantity(z, assume_units)
+        truth = units.Quantity(truth, assume_units)
+
+    img = inverse_distance_to_points(obs_points, z, test_points, kind=method, **extra_kw)
+    assert_array_almost_equal(truth, img)
+
 
-        img = inverse_distance_to_points(obs_points, z, test_points, kind=method, **extra_kw)
+def test_interpolate_to_points_invalid(test_data):
+    """Test that interpolate_to_points raises when given an invalid method."""
+    xp, yp, z = test_data
+    obs_points = np.vstack([xp, yp]).transpose() * 10
 
-        with get_test_data(test_file) as fobj:
-            truth = np.load(fobj)['img'].reshape(-1)
+    with get_test_data('interpolation_test_points.npz') as fobj:
+        test_points = np.load(fobj)['points']
 
-        assert_array_almost_equal(truth, img)
+    with pytest.raises(ValueError):
+        interpolate_to_points(obs_points, z, test_points, interp_type='shouldraise')
 
 
+@pytest.mark.parametrize('assume_units', [None, 'mbar'])
 @pytest.mark.parametrize('method', ['natural_neighbor', 'cressman', 'barnes', 'linear',
-                                    'nearest', 'rbf', 'shouldraise', 'cubic'])
-def test_interpolate_to_points(method, test_data):
+                                    'nearest', 'rbf', 'cubic'])
+def test_interpolate_to_points(method, assume_units, test_data):
     r"""Test main grid interpolation function."""
     xp, yp, z = test_data
     obs_points = np.vstack([xp, yp]).transpose() * 10
@@ -152,23 +161,19 @@ def test_interpolate_to_points(method, test_data):
     with get_test_data('interpolation_test_points.npz') as fobj:
         test_points = np.load(fobj)['points']
 
-    extra_kw = {}
     if method == 'cressman':
-        extra_kw['search_radius'] = 200
-        extra_kw['minimum_neighbors'] = 1
+        extra_kw = {'search_radius': 200, 'minimum_neighbors': 1}
     elif method == 'barnes':
-        extra_kw['search_radius'] = 400
-        extra_kw['minimum_neighbors'] = 1
-        extra_kw['gamma'] = 1
-    elif method == 'shouldraise':
-        with pytest.raises(ValueError):
-            interpolate_to_points(
-                obs_points, z, test_points, interp_type=method, **extra_kw)
-        return
-
-    img = interpolate_to_points(obs_points, z, test_points, interp_type=method, **extra_kw)
+        extra_kw = {'search_radius': 400, 'minimum_neighbors': 1, 'gamma': 1}
+    else:
+        extra_kw = {}
 
     with get_test_data(f'{method}_test.npz') as fobj:
         truth = np.load(fobj)['img'].reshape(-1)
 
+    if assume_units:
+        z = units.Quantity(z, assume_units)
+        truth = units.Quantity(truth, assume_units)
+
+    img = interpolate_to_points(obs_points, z, test_points, interp_type=method, **extra_kw)
     assert_array_almost_equal(truth, img)
diff --git a/tests/io/test_metar.py b/tests/io/test_metar.py
index b2896ab82e..ec342e5631 100644
--- a/tests/io/test_metar.py
+++ b/tests/io/test_metar.py
@@ -10,7 +10,6 @@ import pytest
 
 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
 
@@ -353,29 +352,3 @@ 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)
-
-
-def test_repr():
-    """Test that the TreeNode string representation works."""
-    str1 = 'KSMQ 201953Z AUTO VRB05KT 10SM CLR 03/M10 A3026'
-    tree = parse(str1)
-    rep = repr(tree)
-    assert rep == ("TreeNode1(text='KSMQ 201953Z AUTO VRB05KT 10SM CLR 03/M10 A3026', "
-                   "offset=0, metar=TreeNode(text='', offset=0), "
-                   "siteid=TreeNode(text='KSMQ', offset=0), "
-                   "datetime=TreeNode2(text=' 201953Z', offset=4, "
-                   "sep=TreeNode(text=' ', offset=4)), "
-                   "auto=TreeNode(text=' AUTO', offset=12), "
-                   "wind=TreeNode4(text=' VRB05KT', offset=17, "
-                   "wind_dir=TreeNode(text='VRB', offset=18), "
-                   "wind_spd=TreeNode(text='05', offset=21), "
-                   "gust=TreeNode(text='', offset=23)), "
-                   "vis=TreeNode6(text=' 10SM', offset=25, "
-                   "sep=TreeNode(text=' ', offset=25)), run=TreeNode(text='', offset=30), "
-                   "curwx=TreeNode(text='', offset=30), "
-                   "skyc=TreeNode(text=' CLR', offset=30), "
-                   "temp_dewp=TreeNode13(text=' 03/M10', offset=34, "
-                   "sep=TreeNode(text=' ', offset=34), temp=TreeNode(text='03', offset=35), "
-                   "dewp=TreeNode(text='M10', offset=38)), "
-                   "altim=TreeNode(text=' A3026', offset=41), "
-                   "remarks=TreeNode(text='', offset=47), end=TreeNode(text='', offset=47))")