adding squared angles?

[dcherian]

(I really don't know what the title of this should be; I also don't know if this should be in cf_xarray or here haha)
I am trying to find the closest point but I can't add degrees_east² and degrees_north² even though they're both dimensionally the same?

import pint_xarray
import xarray as xr
from cf_xarray.units import units

lon = xr.DataArray([-130, -140, -150], attrs={"units": "degrees_east"}).pint.quantify()
lat = xr.DataArray([0, 0, 0], attrs={"units": "degrees_north"}).pint.quantify()

lon_target = -140 * units.degrees_east
lat_target = 0 * units.degrees_north


np.argmin(((lat - lat_target) ** 2 + (lon - lon_target) ** 2))
# This works
# np.argmin(
#    ((lat - lat_target).data.magnitude ** 2 + (lon - lon_target).data.magnitude ** 2)
# )

gives this weird error: ValueError: Integers to negative integer powers are not allowed.

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-35-46b602164f7a> in <module>
     11 
     12 np.argmin(
---> 13     ((lat - lat_target) ** 2 + (lon - lon_target) ** 2)
     14 )
     15 #np.argmin(

~/python/xarray/xarray/core/_typed_ops.py in __add__(self, other)
    204 
    205     def __add__(self, other):
--> 206         return self._binary_op(other, operator.add)
    207 
    208     def __sub__(self, other):

~/python/xarray/xarray/core/dataarray.py in _binary_op(self, other, f, reflexive)
   3024 
   3025         variable = (
-> 3026             f(self.variable, other_variable)
   3027             if not reflexive
   3028             else f(other_variable, self.variable)

~/python/xarray/xarray/core/_typed_ops.py in __add__(self, other)
    394 
    395     def __add__(self, other):
--> 396         return self._binary_op(other, operator.add)
    397 
    398     def __sub__(self, other):

~/python/xarray/xarray/core/variable.py in _binary_op(self, other, f, reflexive)
   2335         with np.errstate(all="ignore"):
   2336             new_data = (
-> 2337                 f(self_data, other_data) if not reflexive else f(other_data, self_data)
   2338             )
   2339         result = Variable(dims, new_data, attrs=attrs)

~/miniconda3/envs/dcpy/lib/python3.8/site-packages/pint/quantity.py in __add__(self, other)
   1076             return self.to_timedelta() + other
   1077         else:
-> 1078             return self._add_sub(other, operator.add)
   1079 
   1080     __radd__ = __add__

~/miniconda3/envs/dcpy/lib/python3.8/site-packages/pint/quantity.py in wrapped(self, *args, **kwargs)
    113         elif isinstance(other, list) and other and isinstance(other[0], type(self)):
    114             return NotImplemented
--> 115         return f(self, *args, **kwargs)
    116 
    117     return wrapped

~/miniconda3/envs/dcpy/lib/python3.8/site-packages/pint/quantity.py in _add_sub(self, other, op)
   1013             else:
   1014                 units = self._units
-> 1015                 magnitude = op(self._magnitude, other.to(self._units).magnitude)
   1016 
   1017         elif (

~/miniconda3/envs/dcpy/lib/python3.8/site-packages/pint/quantity.py in to(self, other, *contexts, **ctx_kwargs)
    658         other = to_units_container(other, self._REGISTRY)
    659 
--> 660         magnitude = self._convert_magnitude_not_inplace(other, *contexts, **ctx_kwargs)
    661 
    662         return self.__class__(magnitude, other)

~/miniconda3/envs/dcpy/lib/python3.8/site-packages/pint/quantity.py in _convert_magnitude_not_inplace(self, other, *contexts, **ctx_kwargs)
    607                 return self._REGISTRY.convert(self._magnitude, self._units, other)
    608 
--> 609         return self._REGISTRY.convert(self._magnitude, self._units, other)
    610 
    611     def _convert_magnitude(self, other, *contexts, **ctx_kwargs):

~/miniconda3/envs/dcpy/lib/python3.8/site-packages/pint/registry.py in convert(self, value, src, dst, inplace)
    945             return value
    946 
--> 947         return self._convert(value, src, dst, inplace)
    948 
    949     def _convert(self, value, src, dst, inplace=False, check_dimensionality=True):

~/miniconda3/envs/dcpy/lib/python3.8/site-packages/pint/registry.py in _convert(self, value, src, dst, inplace)
   1828                 value, src = src._magnitude, src._units
   1829 
-> 1830         return super()._convert(value, src, dst, inplace)
   1831 
   1832     def _get_compatible_units(self, input_units, group_or_system):

~/miniconda3/envs/dcpy/lib/python3.8/site-packages/pint/registry.py in _convert(self, value, src, dst, inplace)
   1431 
   1432         if not (src_offset_unit or dst_offset_unit):
-> 1433             return super()._convert(value, src, dst, inplace)
   1434 
   1435         src_dim = self._get_dimensionality(src)

~/miniconda3/envs/dcpy/lib/python3.8/site-packages/pint/registry.py in _convert(self, value, src, dst, inplace, check_dimensionality)
    982         # Here src and dst have only multiplicative units left. Thus we can
    983         # convert with a factor.
--> 984         factor, _ = self._get_root_units(src / dst)
    985 
    986         # factor is type float and if our magnitude is type Decimal then

~/miniconda3/envs/dcpy/lib/python3.8/site-packages/pint/registry.py in _get_root_units(self, input_units, check_nonmult)
    822 
    823         accumulators = [1, defaultdict(int)]
--> 824         self._get_root_units_recurse(input_units, 1, accumulators)
    825 
    826         factor = accumulators[0]

~/miniconda3/envs/dcpy/lib/python3.8/site-packages/pint/registry.py in _get_root_units_recurse(self, ref, exp, accumulators)
    871                 accumulators[1][key] += exp2
    872             else:
--> 873                 accumulators[0] *= reg._converter.scale ** exp2
    874                 if reg.reference is not None:
    875                     self._get_root_units_recurse(reg.reference, exp2, accumulators)

ValueError: Integers to negative integer powers are not allowed.

[jthielen]

This looks to be a very strange issue in Pint with building the right conversion factor between squared "dimensionless dimensionality" units when called internally from a calculation. Try running ((lon - lon_target) ** 2).to('degrees_north ** 2'), and then re-run your np.argmin line, and see if forcing the conversion into the cache takes care of it. In any case, another good issue for upstream in Pint!