Make Quantities work with scipy.interpolate

[olebole]

Description

Currently, one can not use astropy.units.Quantity as within scipys interp1d or interp2d. In interp1d, the units are ignored everywhere:

>>> import numpy as np
>>> import astropy.units as u
>>> from scipy.interpolate import interp1d
>>> x = np.arange(1., 10., 2.)
>>> y = 2 * x
>>> interp1d(x, y)(1.3)
array(2.6)
>>> interp1d(x*u.nm, y)(1.3)
array(2.6)
>>> interp1d(x*u.nm, y)(1.3*u.AA)
array(2.6)
>>> interp1d(x, y*u.Jy)(1.3)
array(2.6)

In interp2d, the unit in the z axis is just ignored, and the others axes don't accept a unit at all:

>>> import numpy as np
>>> import astropy.units as u
>>> from scipy.interpolate import interp2d
>>> x = np.array([0,1,2,0,1,2])
>>> y = np.array([0,0,0,3,3,3])
>>> z = np.array([1,2,3,4,5,6])
>>> interp2d(x, y, z)(1.3, 2.5)
array([4.8])
>>> interp2d(x, y, z*u.Jy)(1.3, 2.5)
array([4.8])
>>> interp2d(x*u.nm, y, z)
Traceback (most recent call last):
  File "/usr/lib/python3/dist-packages/astropy/units/quantity.py", line 727, in to_value
    scale = self.unit._to(unit)
  File "/usr/lib/python3/dist-packages/astropy/units/core.py", line 950, in _to
    raise UnitConversionError(
astropy.units.core.UnitConversionError: 'Unit("nm")' is not a scaled version of 'Unit(dimensionless)'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/usr/lib/python3/dist-packages/astropy/units/quantity.py", line 1073, in __float__
    return float(self.to_value(dimensionless_unscaled))
  File "/usr/lib/python3/dist-packages/astropy/units/quantity.py", line 730, in to_value
    value = self._to_value(unit, equivalencies)
  File "/usr/lib/python3/dist-packages/astropy/units/quantity.py", line 660, in _to_value
    return self.unit.to(unit, self.view(np.ndarray),
  File "/usr/lib/python3/dist-packages/astropy/units/core.py", line 987, in to
    return self._get_converter(other, equivalencies=equivalencies)(value)
  File "/usr/lib/python3/dist-packages/astropy/units/core.py", line 918, in _get_converter
    raise exc
  File "/usr/lib/python3/dist-packages/astropy/units/core.py", line 903, in _get_converter
    return self._apply_equivalencies(
  File "/usr/lib/python3/dist-packages/astropy/units/core.py", line 886, in _apply_equivalencies
    raise UnitConversionError(
astropy.units.core.UnitConversionError: 'nm' (length) and '' (dimensionless) are not convertible

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python3/dist-packages/scipy/interpolate/interpolate.py", line 229, in __init__
    self.tck = fitpack.bisplrep(x, y, z, kx=kx, ky=ky, s=0.0)
  File "/usr/lib/python3/dist-packages/scipy/interpolate/_fitpack_impl.py", line 958, in bisplrep
    tx, ty, c, o = _fitpack._surfit(x, y, z, w, xb, xe, yb, ye, kx, ky,
  File "/usr/lib/python3/dist-packages/astropy/units/quantity.py", line 1075, in __float__
    raise TypeError('only dimensionless scalar quantities can be '
TypeError: only dimensionless scalar quantities can be converted to Python scalars

Additional context

The tests were done with astropy 4.0.1.post1, scipy 1.5.2, numpy 1.19.1.

I am not sure whether this is a bug or a missing feature.

[pllim]

From interp2d doc, x and y are supposed to be "data point coordinates." I am not sure if it is appropriate to try to apply unit to coordinates here. What if someone pass in pixel in one dimension but arcsec in another? How does it translate with the missing plate scale information? Perhaps it not crashing on interp1d was pure luck (despite unit being ignored).

As for z, I feel like native Quantity support would not be trivial as it farms out the actual calculations to lower level internal functions. Maybe @mhvk can comment further. For now, strip out the unit, apply your interpolation, and then re-apply the unit. Hope this helps!

[olebole]

Coordinates may have naturally units. Specific use case: I have a table containing a grid

• moon-sun-separation [deg]
• wavelength [nm]
• some brightness factor

and I want to interpolate the factor on places that are not in the table:

sep = np.array(…) * u.deg
wavelength = np.array(…) * u.nm
fac = np.array(…)
moon_brightness = scipy.interpolate.interp2d(sep, wavelength, fac)

b600_grey = moon_brightness(90*u.deg, 600*u.nm)

wavelength and separation are here data point coordinates, and they naturally come with a unit.

[mhvk]

@olebole - I think unit support is in itself entirely logical, and indeed we already support np.interp, but to get quantities to work requires scipy to start supporting __array_function__ or some other way to override the functions. (I haven't yet looked at what it does right now, but often code does a blunt np.asarray, in which case one just loses the units...)

[mhvk]

Now checked in a bit more detail and, as I feared, interp1d just passes its inputs through np.array ( https://github.com/scipy/scipy/blob/02563b2c873f042f377ec19b08ac6c681b826897/scipy/interpolate/interpolate.py#L442-L443), which means the units get stripped immediately. Oddly enough, for interp2d, the units are stripped of z, but x and y are passed through ravel instead of asarray and hence they keep their units. Basically, like for numpy, the code was not written with subclasses in mind. It is in principle changeable, but it will take time and effort (beyond what I have, I fear; EDIT: I did check that with some changes to integration/integrate.py, integration/polyint.py and _lib/_util.py I could get your 1D examples to work).

I note that __array_function__ support is on the roadmap at https://scipy.github.io/devdocs/roadmap.html, but that may also be a while...

[olebole]

@mhvk, whow, thank you very much for your efforts!

[nstarman]

If you're interested in unit-compatible interpolated splines, I made a gist of some code I use quite often.
It's equivalent to the scipy classes, but with unit support.

https://gist.github.com/nstarman/0628c4a1f80cea32437378f04cee957a

@mhvk
Would something like this be worth making a PR for? I also have a test file with ~100% codecov. I would, of course, modify the code style to more match astropy rather than black.