Repository: Unidata/MetPy

Issue:
UndefinedUnitError: 'dBZ' is not defined in the unit registry
### What went wrong?

Hello everyone! 
I want to create a radar reflectivity profile, but I encountered an issue when using the cross_section function. It seems that the unit 'dBZ' for reflectivity is not defined, although the dataset used has already defined its unit. I would like to know the reason for this problem and hope that you can provide assistance.
![微信截图_20240102014557](https://github.com/Unidata/MetPy/assets/152379644/72281de2-5a8b-4498-977e-fb9a5b7e193f)



### Operating System

Windows

### Version

1.6.0

### Python Version

3.11.4

### Code to Reproduce

```python
da = xr.open_dataset(r'E:\data\0415\20210415_090551_vdras.nc')
data = da.metpy.assign_crs(
    grid_mapping_name='lambert_conformal_conic',
    latitude_of_projection_origin=39.808887,
    longitude_of_central_meridian=116.47195,
    standard_parallel=39.808887,
    earth_radius=6371.229
).metpy.assign_latitude_longitude()
start = (116.6,39.87)
end = (117.03,39.6)
print(data)
cross = cross_section(data, start, end).set_coords(('lat', 'lon'))
```


### Errors, Traceback, and Logs

```pytb
UndefinedUnitError                        Traceback (most recent call last)
Cell In[75], line 12
     10 end = (117.03,39.6)
     11 print(data)
---> 12 cross = cross_section(data, start, end).set_coords(('lat', 'lon'))

File d:\anaconda3\Lib\site-packages\metpy\interpolate\slices.py:147, in cross_section(data, start, end, steps, interp_type)
    109 r"""Obtain an interpolated cross-sectional slice through gridded data.
    110 
    111 Utilizing the interpolation functionality in xarray, this function takes a vertical
   (...)
    143 
    144 """
    145 if isinstance(data, xr.Dataset):
    146     # Recursively apply to dataset
--> 147     return data.map(cross_section, True, (start, end), steps=steps,
    148                     interp_type=interp_type)
    149 elif data.ndim == 0:
    150     # This has no dimensions, so it is likely a projection variable. In any case, there
    151     # are no data here to take the cross section with. Therefore, do nothing.
    152     return data

File d:\anaconda3\Lib\site-packages\xarray\core\dataset.py:6037, in Dataset.map(self, func, keep_attrs, args, **kwargs)
   6035 if keep_attrs is None:
   6036     keep_attrs = _get_keep_attrs(default=False)
-> 6037 variables = {
   6038     k: maybe_wrap_array(v, func(v, *args, **kwargs))
   6039     for k, v in self.data_vars.items()
   6040 }
   6041 if keep_attrs:
   6042     for k, v in variables.items():

File d:\anaconda3\Lib\site-packages\xarray\core\dataset.py:6038, in <dictcomp>(.0)
   6035 if keep_attrs is None:
   6036     keep_attrs = _get_keep_attrs(default=False)
   6037 variables = {
-> 6038     k: maybe_wrap_array(v, func(v, *args, **kwargs))
   6039     for k, v in self.data_vars.items()
   6040 }
   6041 if keep_attrs:
   6042     for k, v in variables.items():

File d:\anaconda3\Lib\site-packages\metpy\interpolate\slices.py:173, in cross_section(data, start, end, steps, interp_type)
    170     points_cross[points_cross[:, 0] < 0, 0] += 360.
    172 # Return the interpolated data
--> 173 return interpolate_to_slice(data, points_cross, interp_type=interp_type)

File d:\anaconda3\Lib\site-packages\metpy\interpolate\slices.py:59, in interpolate_to_slice(data, points, interp_type)
     53 data_sliced = data.interp({
     54     x.name: xr.DataArray(points[:, 0], dims='index', attrs=x.attrs),
     55     y.name: xr.DataArray(points[:, 1], dims='index', attrs=y.attrs)
     56 }, method=interp_type)
     57 data_sliced.coords['index'] = range(len(points))
---> 59 return data_sliced.metpy.quantify()

File d:\anaconda3\Lib\site-packages\metpy\xarray.py:235, in MetPyDataArrayAccessor.quantify(self)
    222 """Return a new DataArray with the data converted to a `pint.Quantity`.
    223 
    224 Notes
   (...)
    228 subsetting!
    229 """
    230 if (
    231     not is_quantity(self._data_array.data)
    232     and np.issubdtype(self._data_array.data.dtype, np.number)
    233 ):
    234     # Only quantify if not already quantified and is quantifiable
--> 235     quantified_dataarray = self._data_array.copy(data=self.unit_array)
    236     if 'units' in quantified_dataarray.attrs:
    237         del quantified_dataarray.attrs['units']

File d:\anaconda3\Lib\site-packages\metpy\xarray.py:165, in MetPyDataArrayAccessor.unit_array(self)
    163     return self._data_array.data
    164 else:
--> 165     return units.Quantity(self._data_array.data, self.units)

File d:\anaconda3\Lib\site-packages\metpy\xarray.py:142, in MetPyDataArrayAccessor.units(self)
    140 else:
    141     default_unit = 'dimensionless'
--> 142 return units.parse_units(self._data_array.attrs.get('units', default_unit))

File d:\anaconda3\Lib\site-packages\pint\facets\plain\registry.py:1189, in GenericPlainRegistry.parse_units(self, input_string, as_delta, case_sensitive)
   1161 def parse_units(
   1162     self,
   1163     input_string: str,
   1164     as_delta: Optional[bool] = None,
   1165     case_sensitive: Optional[bool] = None,
   1166 ) -> UnitT:
   1167     """Parse a units expression and returns a UnitContainer with
   1168     the canonical names.
   1169 
   (...)
   1185 
   1186     """
   1188     return self.Unit(
-> 1189         self.parse_units_as_container(input_string, as_delta, case_sensitive)
   1190     )

File d:\anaconda3\Lib\site-packages\pint\facets\nonmultiplicative\registry.py:70, in GenericNonMultiplicativeRegistry.parse_units_as_container(self, input_string, as_delta, case_sensitive)
     67 if as_delta is None:
     68     as_delta = self.default_as_delta
---> 70 return super().parse_units_as_container(input_string, as_delta, case_sensitive)

File d:\anaconda3\Lib\site-packages\pint\facets\plain\registry.py:1204, in GenericPlainRegistry.parse_units_as_container(self, input_string, as_delta, case_sensitive)
   1198 as_delta = (
   1199     as_delta if as_delta is not None else True
   1200 )  # TODO This only exists in nonmultiplicative
   1201 case_sensitive = (
   1202     case_sensitive if case_sensitive is not None else self.case_sensitive
   1203 )
-> 1204 return self._parse_units_as_container(input_string, as_delta, case_sensitive)

File d:\anaconda3\Lib\site-packages\pint\facets\plain\registry.py:1239, in GenericPlainRegistry._parse_units_as_container(self, input_string, as_delta, case_sensitive)
   1237 many = len(units) > 1
   1238 for name in units:
-> 1239     cname = self.get_name(name, case_sensitive=case_sensitive)
   1240     value = units[name]
   1241     if not cname:

File d:\anaconda3\Lib\site-packages\pint\facets\plain\registry.py:647, in GenericPlainRegistry.get_name(self, name_or_alias, case_sensitive)
    645 candidates = self.parse_unit_name(name_or_alias, case_sensitive)
    646 if not candidates:
--> 647     raise UndefinedUnitError(name_or_alias)
    649 prefix, unit_name, _ = candidates[0]
    650 if len(candidates) > 1:

UndefinedUnitError: 'dBZ' is not defined in the unit registry
```


Discussion:

Patch:
diff --git a/.gitignore b/.gitignore
index 3b8f0987ca..b4ee237550 100644
--- a/.gitignore
+++ b/.gitignore
@@ -29,3 +29,4 @@ docs/api/generated
 docs/api/areas.rst
 examples/scripts
 test_output/
+docs/sg_execution_times.rst
diff --git a/src/metpy/interpolate/slices.py b/src/metpy/interpolate/slices.py
index b11199565e..2e8799f2fa 100644
--- a/src/metpy/interpolate/slices.py
+++ b/src/metpy/interpolate/slices.py
@@ -7,6 +7,7 @@ import numpy as np
 import xarray as xr
 
 from ..package_tools import Exporter
+from ..units import is_quantity
 from ..xarray import check_axis
 
 exporter = Exporter(globals())
@@ -49,6 +50,7 @@ def interpolate_to_slice(data, points, interp_type='linear'):
                          'your data has been parsed by MetPy with proper x and y '
                          'dimension coordinates.') from None
 
+    need_quantify = is_quantity(data.data)
     data = data.metpy.dequantify()
     data_sliced = data.interp({
         x.name: xr.DataArray(points[:, 0], dims='index', attrs=x.attrs),
@@ -56,7 +58,7 @@ def interpolate_to_slice(data, points, interp_type='linear'):
     }, method=interp_type)
     data_sliced.coords['index'] = range(len(points))
 
-    return data_sliced.metpy.quantify()
+    return data_sliced.metpy.quantify() if need_quantify else data_sliced
 
 
 @exporter.export


Test patch:
diff --git a/tests/interpolate/test_slices.py b/tests/interpolate/test_slices.py
index 41eb6dc01f..5407857c0b 100644
--- a/tests/interpolate/test_slices.py
+++ b/tests/interpolate/test_slices.py
@@ -95,16 +95,24 @@ def test_ds_xy():
     return ds.metpy.parse_cf()
 
 
-def test_interpolate_to_slice_against_selection(test_ds_lonlat):
+@pytest.mark.parametrize('bad_units', [False, True])
+def test_interpolate_to_slice_against_selection(test_ds_lonlat, bad_units):
     """Test interpolate_to_slice on a simple operation."""
     data = test_ds_lonlat['temperature']
+
+    # interpolate_to_slice shouldn't care about units
+    if bad_units:
+        # Needed so we can go back to using attribute metadata
+        data = data.metpy.dequantify()
+        data.attrs['units'] = 'my_bad_units'
+
     path = np.array([[265.0, 30.],
                      [265.0, 36.],
                      [265.0, 42.]])
     test_slice = interpolate_to_slice(data, path)
     true_slice = data.sel({'lat': [30., 36., 42.], 'lon': 265.0})
     # Coordinates differ, so just compare the data
-    assert_array_almost_equal(true_slice.metpy.unit_array, test_slice.metpy.unit_array, 5)
+    assert_array_almost_equal(true_slice.data, test_slice.data, 5)
 
 
 @needs_cartopy