Support arrays for nside

[cdeil]

One of the issues trying to switch Gammapy over to using astropy-healpix is that we call ang2pix with an array of nside values: https://gist.github.com/cdeil/754b76dc7f22511a5504fbbe74dccd62#file-gistfile1-txt-L989

Minimum test case that works with healpy

>>> import healpy as hp
>>> hp.ang2pix(nside=[2, 4], theta=[0, 0], phi=[0, 0])
array([0, 0])

but doesn't work with astropy_healpix.healpy (casting to numpy arrays is missing):

>>> from astropy_healpix import healpy as hp
>>> hp.ang2pix(nside=[2, 4], theta=[0, 0], phi=[0, 0])
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/deil/code/astropy-healpix/astropy_healpix/healpy.py", line 105, in ang2pix
    return lonlat_to_healpix(lon, lat, nside, order='nested' if nest else 'ring')
  File "/Users/deil/code/astropy-healpix/astropy_healpix/core.py", line 341, in lonlat_to_healpix
    nside = int(nside)
TypeError: int() argument must be a string, a bytes-like object or a number, not 'list'

and also support for array-valued nside is missing:

>>> hp.ang2pix(nside=np.array([2, 4]), theta=np.array([0, 0]), phi=np.array([0, 0]))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/deil/code/astropy-healpix/astropy_healpix/healpy.py", line 105, in ang2pix
    return lonlat_to_healpix(lon, lat, nside, order='nested' if nest else 'ring')
  File "/Users/deil/code/astropy-healpix/astropy_healpix/core.py", line 341, in lonlat_to_healpix
    nside = int(nside)
TypeError: only length-1 arrays can be converted to Python scalars

I didn't look yet how much work it would be to make this work for ang2pix or possibly also similar other functions, I'm just going through and see what's missing for Gammapy today.

[astrofrog]

In principle there's no reason we can't make the Cython wrappers take a vector nside and it would be easy. My main concern is that we need to keep it efficient for the scalar case and avoid creating an unnecessary array for nside.

Could you explain why gammapy needs this, to see if there is another way we can do this efficiently?

[astrofrog]

I guess what I'm thinking is that for now it wouldn't be that inefficient to simply call the astropy-healpix routines once for each nside value, since there are only 30 possible nside values.

[astrofrog]

Actually, maybe before worrying we should just try and change the Cython wrappers to use vector nside values and then if a scalar nside is passed in core.py, convert to an array. We might find the performance penalty is not actually large.

[lpsinger]

I need this too because I use healpy for supporting NUNIQ indexing.

[cdeil]

As commented by @lpsinger in https://github.com/astropy/astropy-healpix/pull/71#issuecomment-401027171 - this should still be implemented, but we should change to ufuncs first.