Healpix Representation

[nstarman]

I don't have time to work on this much in the foreseeable future. So feel free to take this over the finish line or leave as is, no rush either way.
This PR is an early and incomplete implementation of healpix as a Representation. No more need for special methods like xyz_to_healpix. Using the representation machinery a healpix rep can transform to any Astropy representation, e.g. PhysicsSpherical. Also, since Representations can be inside a Frame or SkyCoord, this is super easy to work with.
If implemented, this would be nice to use as the backend in HEALPix. IMO there's a lot of functions in this package that can be organized into methods on data structures. I've used some of them here, but the actual implementation could be moved to the Representation object and the separate function deprecated.

Example usage:

>>> HEALPixNested8Representation = HEALPixBaseRepresentation.get_subclass(order="nested", nside=8)
>>> r = HEALPixNested8Representation([10, 12], 2*u.kpc, 0.5, 0.5)
>>> r
<HEALPixNested8Representation (indices, distance, dx, dy) in (, kpc, , )
    [(10, 2., 0.5, 0.5), (12, 2., 0.5, 0.5)]>

>>> c = coord.ICRS(r, representation_type=HEALPixNested8Representation)
>>> c
<ICRS Coordinate: (indices, distance, dx, dy) in (, kpc, , )
    [(10, 2., 0.5, 0.5), (12, 2., 0.5, 0.5)]>

>>> r.represent_as(coord.PhysicsSphericalRepresentation)
<PhysicsSphericalRepresentation (phi, theta, r) in (rad, rad, kpc)
    [(0.49087385, 1.23095942, 2.), (0.78539816, 1.1410209 , 2.)]>

Signed-off-by: nstarman nstarkman@protonmail.com

[lpsinger]

I love this idea! I tried to do something like this many years ago, but never got it over the finish line. I am excited to review this.

[pllim]

@lpsinger , do you want this in your 1.0 ?

[lpsinger]

• This needs documentation.
• Can we support nside being an array? This would be important for supporting multi-resolution data sets.
• Expose pixel area, resolution, etc. as instance methods.
• Indices should be of dtype np.intp so that they are suitable for indexing. What do we need to do to enforce that? Currently they are a (presumably dimensionless) u.Quantity.
• What should the units for dx and dy be?

[lpsinger]

@lpsinger , do you want this in your 1.0 ?

Hmm, I forgot about this! Do you think that this would eventually replace the existing high-level HEALPix class?

[nstarman]

I regret making this a classmethod. This should be a standalone function.

[nstarman]

Thanks @lpsinger. I opened this PR mostly as a demonstration; I think it's a good idea, but requires a fair bit of work and time I don't have ATM. I'm happy to chat / review PRs.

• This needs documentation.

Agreed.

• Can we support nside being an array? This would be important for supporting multi-resolution data sets.

Yes, with caveats. The Representation machinery requires its subclasses to have the coordinates as the only constructor arguments. e.g. x,y,z -> rho, theta,phi. In the Healpix Representation, these are indices, dx, dy.
nside and norder are not coordinates, but part of the definition of the Representation type, which is why they are set by get_subclass. So a HealpixRepresentation[nside=4, order=6] is a different representation type then HealpixRepresentation[nside=4, order=8], though their constructor arguments are the same.
We can absolutely relax the requirement that nside is a single int to instead nside: int | tuple[int, ...]. The only complication is then that nside would have to be the same length as the coordinates — len(nside) == len(indices. So the class type HealpixRepresentation[nside=(4, 5, 6, 8), order=6] can only accept coordinate arrays of length 4, matching nside.

• Expose pixel area, resolution, etc. as instance methods.

Sounds good. Though see the last point.

• Indices should be of dtype np.intp so that they are suitable for indexing. What do we need to do to enforce that? Currently they are a (presumably dimensionless) u.Quantity.

Aren't they already set as such? In __init__ I have

self._indices = (self._indices << u.one).astype(int)

So indices is a Quantity with dtype int.

• What should the units for dx and dy be?

Currently the Representation machinery requires all coordinates to be Quantity. I have the units set to dimensionless_unscaled ATM, but that can be changed. Should they be radians?

Do you think that this would eventually replace the existing high-level HEALPix class?

Possibly! Or, like CoordinateFrame is a composed structure, containing Representation objects, the high-level HEALPix class might hold a HEALPixRepresentation object and expose useful methods, e.g. pixel_area, resolution, etc.
Or probably better, the high-level HEALPix class is more akin to SkyCoord, which holds a CoordinateFrame, which holds a Representation. A nice thing about the HealpixRepresentation in this PR is that it works with CoordinateFrame!
My inclination is to keep the Representation atomic, like the ones in Astropy core. If the representation doesn't hold all these convenience functions, then they must be stand-alone functions or attributes/methods on another class.

[lpsinger]

Thanks @lpsinger. I opened this PR mostly as a demonstration; I think it's a good idea, but requires a fair bit of work and time I don't have ATM. I'm happy to chat / review PRs.

OK, in that case, I'll proceed with releasing v1.0.0 with what we have now.