# Licensed under a 3-clause BSD style license - see LICENSE.rst
import math
import numpy as np
from astropy import units as u
from astropy.coordinates import Longitude, Latitude
from . import _core
__all__ = [
'nside_to_pixel_area',
'nside_to_pixel_resolution',
'pixel_resolution_to_nside',
'nside_to_npix',
'npix_to_nside',
'level_to_nside',
'nside_to_level',
'level_ipix_to_uniq',
'uniq_to_level_ipix',
'lonlat_to_healpix',
'healpix_to_lonlat',
'xyz_to_healpix',
'healpix_to_xyz',
'bilinear_interpolation_weights',
'interpolate_bilinear_lonlat',
'boundaries_lonlat',
'neighbours',
]
_NUMPY_COPY_IF_NEEDED = False if np.__version__.startswith("1.") else None
def _validate_order(order):
# We also support upper-case, to support directly the values
# ORDERING = {'RING', 'NESTED'} in FITS headers
# This is currently undocumented in the docstrings.
if order == 'nested' or order == 'NESTED':
return 'nested'
elif order == 'ring' or order == 'RING':
return 'ring'
else:
raise ValueError("order must be 'nested' or 'ring'")
def _validate_offset(label, offset):
offset = np.asarray(offset)
if np.any((offset < 0) | (offset > 1)):
raise ValueError(f'd{label} must be in the range [0:1]')
def _validate_level(level):
if np.any(level < 0):
raise ValueError('level must be positive')
def _validate_nside(nside):
log_2_nside = np.round(np.log2(nside))
if not np.all(2 ** log_2_nside == nside):
raise ValueError('nside must be a power of two')
def _validate_npix(level, ipix):
if not np.all(ipix < (3 << 2*(level + 1))):
raise ValueError('ipix for a specific level must be inferior to npix')
[docs]
def level_to_nside(level):
"""
Find the pixel dimensions of the top-level HEALPix tiles.
This is given by ``nside = 2**level``.
Parameters
----------
level : int
The resolution level
Returns
-------
nside : int
The number of pixels on the side of one of the 12 'top-level' HEALPix tiles.
"""
level = np.asarray(level, dtype=np.int64)
_validate_level(level)
return 2 ** level
[docs]
def nside_to_level(nside):
"""
Find the HEALPix level for a given nside.
This is given by ``level = log2(nside)``.
This function is the inverse of `level_to_nside`.
Parameters
----------
nside : int
The number of pixels on the side of one of the 12 'top-level' HEALPix tiles.
Must be a power of two.
Returns
-------
level : int
The level of the HEALPix cells
"""
nside = np.asarray(nside, dtype=np.int64)
_validate_nside(nside)
return np.log2(nside).astype(np.int64)
[docs]
def uniq_to_level_ipix(uniq):
"""
Convert a HEALPix cell uniq number to its (level, ipix) equivalent.
A uniq number is a 64 bits integer equaling to : ipix + 4*(4**level). Please read
this `paper <http://ivoa.net/documents/MOC/20140602/REC-MOC-1.0-20140602.pdf>`_
for more details about uniq numbers.
Parameters
----------
uniq : int
The uniq number of a HEALPix cell.
Returns
-------
level, ipix: int, int
The level and index of the HEALPix cell computed from ``uniq``.
"""
uniq = np.asarray(uniq, dtype=np.int64)
level = (np.log2(uniq//4)) // 2
level = level.astype(np.int64)
_validate_level(level)
ipix = uniq - (1 << 2*(level + 1))
_validate_npix(level, ipix)
return level, ipix
[docs]
def level_ipix_to_uniq(level, ipix):
"""
Convert a level and HEALPix index into a uniq number representing the cell.
This function is the inverse of `uniq_to_level_ipix`.
Parameters
----------
level : int
The level of the HEALPix cell
ipix : int
The index of the HEALPix cell
Returns
-------
uniq : int
The uniq number representing the HEALPix cell.
"""
level = np.asarray(level, dtype=np.int64)
ipix = np.asarray(ipix, dtype=np.int64)
_validate_level(level)
_validate_npix(level, ipix)
return ipix + (1 << 2*(level + 1))
[docs]
def nside_to_pixel_area(nside):
"""
Find the area of HEALPix pixels given the pixel dimensions of one of
the 12 'top-level' HEALPix tiles.
Parameters
----------
nside : int
The number of pixels on the side of one of the 12 'top-level' HEALPix tiles.
Returns
-------
pixel_area : :class:`~astropy.units.Quantity`
The area of the HEALPix pixels
"""
nside = np.asanyarray(nside, dtype=np.int64)
_validate_nside(nside)
npix = 12 * nside * nside
pixel_area = 4 * math.pi / npix * u.sr
return pixel_area
[docs]
def nside_to_pixel_resolution(nside):
"""
Find the resolution of HEALPix pixels given the pixel dimensions of one of
the 12 'top-level' HEALPix tiles.
Parameters
----------
nside : int
The number of pixels on the side of one of the 12 'top-level' HEALPix tiles.
Returns
-------
resolution : :class:`~astropy.units.Quantity`
The resolution of the HEALPix pixels
See also
--------
pixel_resolution_to_nside
"""
nside = np.asanyarray(nside, dtype=np.int64)
_validate_nside(nside)
return (nside_to_pixel_area(nside) ** 0.5).to(u.arcmin)
[docs]
def pixel_resolution_to_nside(resolution, round='nearest'):
"""Find closest HEALPix nside for a given angular resolution.
This function is the inverse of `nside_to_pixel_resolution`,
for the default rounding scheme of ``round='nearest'``.
If you choose ``round='up'``, you'll get HEALPix pixels that
have at least the requested resolution (usually a bit better
due to rounding).
Pixel resolution is defined as square root of pixel area.
Parameters
----------
resolution : `~astropy.units.Quantity`
Angular resolution
round : {'up', 'nearest', 'down'}
Which way to round
Returns
-------
nside : int
The number of pixels on the side of one of the 12 'top-level' HEALPix tiles.
Always a power of 2.
Examples
--------
>>> from astropy import units as u
>>> from astropy_healpix import pixel_resolution_to_nside
>>> pixel_resolution_to_nside(13 * u.arcmin)
256
>>> pixel_resolution_to_nside(13 * u.arcmin, round='up')
512
"""
resolution = resolution.to(u.rad).value
pixel_area = resolution * resolution
npix = 4 * math.pi / pixel_area
nside = np.sqrt(npix / 12)
# Now we have to round to the closest ``nside``
# Since ``nside`` must be a power of two,
# we first compute the corresponding ``level = log2(nside)`
# round the level and then go back to nside
level = np.log2(nside)
if round == 'up':
level = np.ceil(level)
elif round == 'nearest':
level = np.round(level)
elif round == 'down':
level = np.floor(level)
else:
raise ValueError(f'Invalid value for round: {round!r}')
# For very low requested resolution (i.e. large angle values), we
# return ``level=0``, i.e. ``nside=1``, i.e. the lowest resolution
# that exists with HEALPix
level = np.clip(level.astype(int), 0, None)
return level_to_nside(level)
[docs]
def nside_to_npix(nside):
"""
Find the number of pixels corresponding to a HEALPix resolution.
Parameters
----------
nside : int
The number of pixels on the side of one of the 12 'top-level' HEALPix tiles.
Returns
-------
npix : int
The number of pixels in the HEALPix map.
"""
nside = np.asanyarray(nside, dtype=np.int64)
_validate_nside(nside)
return 12 * nside ** 2
[docs]
def npix_to_nside(npix):
"""
Find the number of pixels on the side of one of the 12 'top-level' HEALPix
tiles given a total number of pixels.
Parameters
----------
npix : int
The number of pixels in the HEALPix map.
Returns
-------
nside : int
The number of pixels on the side of one of the 12 'top-level' HEALPix tiles.
"""
npix = np.asanyarray(npix, dtype=np.int64)
if not np.all(npix % 12 == 0):
raise ValueError('Number of pixels must be divisible by 12')
square_root = np.sqrt(npix / 12)
if not np.all(square_root ** 2 == npix / 12):
raise ValueError('Number of pixels is not of the form 12 * nside ** 2')
return np.round(square_root).astype(int)
[docs]
def healpix_to_lonlat(healpix_index, nside, dx=None, dy=None, order='ring'):
"""
Convert HEALPix indices (optionally with offsets) to longitudes/latitudes.
If no offsets (``dx`` and ``dy``) are provided, the coordinates will default
to those at the center of the HEALPix pixels.
Parameters
----------
healpix_index : int or `~numpy.ndarray`
HEALPix indices (as a scalar or array)
nside : int or `~numpy.ndarray`
Number of pixels along the side of each of the 12 top-level HEALPix tiles
dx, dy : float or `~numpy.ndarray`, optional
Offsets inside the HEALPix pixel, which must be in the range [0:1],
where 0.5 is the center of the HEALPix pixels (as scalars or arrays)
order : { 'nested' | 'ring' }, optional
Order of HEALPix pixels
Returns
-------
lon : :class:`~astropy.coordinates.Longitude`
The longitude values
lat : :class:`~astropy.coordinates.Latitude`
The latitude values
"""
_validate_nside(nside)
if _validate_order(order) == 'ring':
func = _core.healpix_ring_to_lonlat
else: # _validate_order(order) == 'nested'
func = _core.healpix_nested_to_lonlat
if dx is None:
dx = 0.5
else:
_validate_offset('x', dx)
if dy is None:
dy = 0.5
else:
_validate_offset('y', dy)
nside = np.asarray(nside, dtype=np.intc)
lon, lat = func(healpix_index, nside, dx, dy)
lon = Longitude(lon, unit=u.rad, copy=_NUMPY_COPY_IF_NEEDED)
lat = Latitude(lat, unit=u.rad, copy=_NUMPY_COPY_IF_NEEDED)
return lon, lat
[docs]
def lonlat_to_healpix(lon, lat, nside, return_offsets=False, order='ring'):
"""
Convert longitudes/latitudes to HEALPix indices
Parameters
----------
lon, lat : :class:`~astropy.units.Quantity`
The longitude and latitude values as :class:`~astropy.units.Quantity`
instances with angle units.
nside : int or `~numpy.ndarray`
Number of pixels along the side of each of the 12 top-level HEALPix tiles
order : { 'nested' | 'ring' }
Order of HEALPix pixels
return_offsets : bool, optional
If `True`, the returned values are the HEALPix pixel indices as well as
``dx`` and ``dy``, the fractional positions inside the pixels. If
`False` (the default), only the HEALPix pixel indices is returned.
Returns
-------
healpix_index : int or `~numpy.ndarray`
The HEALPix indices
dx, dy : `~numpy.ndarray`
Offsets inside the HEALPix pixel in the range [0:1], where 0.5 is the
center of the HEALPix pixels
"""
if _validate_order(order) == 'ring':
func = _core.lonlat_to_healpix_ring
else: # _validate_order(order) == 'nested'
func = _core.lonlat_to_healpix_nested
nside = np.asarray(nside, dtype=np.intc)
lon = lon.to_value(u.rad)
lat = lat.to_value(u.rad)
healpix_index, dx, dy = func(lon, lat, nside)
if return_offsets:
return healpix_index, dx, dy
else:
return healpix_index
[docs]
def healpix_to_xyz(healpix_index, nside, dx=None, dy=None, order='ring'):
"""
Convert HEALPix indices (optionally with offsets) to Cartesian coordinates.
If no offsets (``dx`` and ``dy``) are provided, the coordinates will default
to those at the center of the HEALPix pixels.
Parameters
----------
healpix_index : int or `~numpy.ndarray`
HEALPix indices (as a scalar or array)
nside : int or `~numpy.ndarray`
Number of pixels along the side of each of the 12 top-level HEALPix tiles
dx, dy : float or `~numpy.ndarray`, optional
Offsets inside the HEALPix pixel, which must be in the range [0:1],
where 0.5 is the center of the HEALPix pixels (as scalars or arrays)
order : { 'nested' | 'ring' }, optional
Order of HEALPix pixels
Returns
-------
x, y, z : float or `~numpy.ndarray`
The Cartesian coordinate components
"""
_validate_nside(nside)
if _validate_order(order) == 'ring':
func = _core.healpix_ring_to_xyz
else: # _validate_order(order) == 'nested'
func = _core.healpix_nested_to_xyz
if dx is None:
dx = 0.5
else:
_validate_offset('x', dx)
if dy is None:
dy = 0.5
else:
_validate_offset('y', dy)
nside = np.asarray(nside, dtype=np.intc)
return func(healpix_index, nside, dx, dy)
[docs]
def xyz_to_healpix(x, y, z, nside, return_offsets=False, order='ring'):
"""
Convert longitudes/latitudes to HEALPix indices
Parameters
----------
x, y, z : float or `~numpy.ndarray`
The Cartesian coordinate components
nside : int or `~numpy.ndarray`
Number of pixels along the side of each of the 12 top-level HEALPix tiles
order : { 'nested' | 'ring' }
Order of HEALPix pixels
return_offsets : bool, optional
If `True`, the returned values are the HEALPix pixel indices as well as
``dx`` and ``dy``, the fractional positions inside the pixels. If
`False` (the default), only the HEALPix pixel indices is returned.
Returns
-------
healpix_index : int or `~numpy.ndarray`
The HEALPix indices
dx, dy : `~numpy.ndarray`
Offsets inside the HEALPix pixel in the range [0:1], where 0.5 is the
center of the HEALPix pixels
"""
if _validate_order(order) == 'ring':
func = _core.xyz_to_healpix_ring
else: # _validate_order(order) == 'nested'
func = _core.xyz_to_healpix_nested
nside = np.asarray(nside, dtype=np.intc)
healpix_index, dx, dy = func(x, y, z, nside)
if return_offsets:
return healpix_index, dx, dy
else:
return healpix_index
def nested_to_ring(nested_index, nside):
"""
Convert a HEALPix 'nested' index to a HEALPix 'ring' index
Parameters
----------
nested_index : int or `~numpy.ndarray`
Healpix index using the 'nested' ordering
nside : int or `~numpy.ndarray`
Number of pixels along the side of each of the 12 top-level HEALPix tiles
Returns
-------
ring_index : int or `~numpy.ndarray`
Healpix index using the 'ring' ordering
"""
nside = np.asarray(nside, dtype=np.intc)
return _core.nested_to_ring(nested_index, nside)
def ring_to_nested(ring_index, nside):
"""
Convert a HEALPix 'ring' index to a HEALPix 'nested' index
Parameters
----------
ring_index : int or `~numpy.ndarray`
Healpix index using the 'ring' ordering
nside : int or `~numpy.ndarray`
Number of pixels along the side of each of the 12 top-level HEALPix tiles
Returns
-------
nested_index : int or `~numpy.ndarray`
Healpix index using the 'nested' ordering
"""
nside = np.asarray(nside, dtype=np.intc)
return _core.ring_to_nested(ring_index, nside)
[docs]
def bilinear_interpolation_weights(lon, lat, nside, order='ring'):
"""
Get the four neighbours for each (lon, lat) position and the weight
associated with each one for bilinear interpolation.
Parameters
----------
lon, lat : :class:`~astropy.units.Quantity`
The longitude and latitude values as :class:`~astropy.units.Quantity`
instances with angle units.
nside : int
Number of pixels along the side of each of the 12 top-level HEALPix tiles
order : { 'nested' | 'ring' }
Order of HEALPix pixels
Returns
-------
indices : `~numpy.ndarray`
2-D array with shape (4, N) giving the four indices to use for the
interpolation
weights : `~numpy.ndarray`
2-D array with shape (4, N) giving the four weights to use for the
interpolation
"""
lon = lon.to_value(u.rad)
lat = lat.to_value(u.rad)
_validate_nside(nside)
nside = np.asarray(nside, dtype=np.intc)
result = _core.bilinear_interpolation_weights(lon, lat, nside)
indices = np.stack(result[:4])
weights = np.stack(result[4:])
if _validate_order(order) == 'nested':
indices = ring_to_nested(indices, nside)
return indices, weights
[docs]
def interpolate_bilinear_lonlat(lon, lat, values, order='ring'):
"""
Interpolate values at specific longitudes/latitudes using bilinear interpolation
Parameters
----------
lon, lat : :class:`~astropy.units.Quantity`
The longitude and latitude values as :class:`~astropy.units.Quantity` instances
with angle units.
values : `~numpy.ndarray`
Array with the values in each HEALPix pixel. The first dimension should
have length 12 * nside ** 2 (and nside is determined automatically from
this).
order : { 'nested' | 'ring' }
Order of HEALPix pixels
Returns
-------
result : float `~numpy.ndarray`
The interpolated values
"""
nside = npix_to_nside(values.shape[0])
indices, weights = bilinear_interpolation_weights(lon, lat, nside, order=order)
values = values[indices]
# At this point values has shape (N, M) where both N and M might be several
# dimensions, and weights has shape (N,), so we need to transpose in order
# to benefit from broadcasting, then transpose back so that the dimension
# with length 4 is at the start again, ready for summing.
result = (values.T * weights.T).T
return result.sum(axis=0)
[docs]
def neighbours(healpix_index, nside, order='ring'):
"""
Find all the HEALPix pixels that are the neighbours of a HEALPix pixel
Parameters
----------
healpix_index : `~numpy.ndarray`
Array of HEALPix pixels
nside : int
Number of pixels along the side of each of the 12 top-level HEALPix tiles
order : { 'nested' | 'ring' }
Order of HEALPix pixels
Returns
-------
neigh : `~numpy.ndarray`
Array giving the neighbours starting SW and rotating clockwise. This has
one extra dimension compared to ``healpix_index`` - the first dimension -
which is set to 8. For example if healpix_index has shape (2, 3),
``neigh`` has shape (8, 2, 3).
Notes
-----
Some HEALPix pixels do not have all 8 neighbours. In these cases, the
corresponding entry in the returned array has the value of -1 and Numpy
may print an invalid value warning. To suppress the warning, use
:class:`numpy.errstate`.
"""
_validate_nside(nside)
nside = np.asarray(nside, dtype=np.intc)
if _validate_order(order) == 'ring':
func = _core.neighbours_ring
else: # _validate_order(order) == 'nested'
func = _core.neighbours_nested
return np.stack(func(healpix_index, nside))
def healpix_cone_search(lon, lat, radius, nside, order='ring'):
"""
Find all the HEALPix pixels within a given radius of a longitude/latitude.
Note that this returns all pixels that overlap, including partially, with
the search cone. This function can only be used for a single lon/lat pair at
a time, since different calls to the function may result in a different
number of matches.
Parameters
----------
lon, lat : :class:`~astropy.units.Quantity`
The longitude and latitude to search around
radius : :class:`~astropy.units.Quantity`
The search radius
nside : int
Number of pixels along the side of each of the 12 top-level HEALPix tiles
order : { 'nested' | 'ring' }
Order of HEALPix pixels
Returns
-------
healpix_index : `~numpy.ndarray`
1-D array with all the matching HEALPix pixel indices.
"""
lon = lon.to_value(u.deg)
lat = lat.to_value(u.deg)
radius = radius.to_value(u.deg)
_validate_nside(nside)
order = _validate_order(order)
return _core.healpix_cone_search(lon, lat, radius, nside, order)
[docs]
def boundaries_lonlat(healpix_index, step, nside, order='ring'):
"""
Return the longitude and latitude of the edges of HEALPix pixels
This returns the longitude and latitude of points along the edge of each
HEALPIX pixel. The number of points returned for each pixel is ``4 * step``,
so setting ``step`` to 1 returns just the corners.
Parameters
----------
healpix_index : `~numpy.ndarray`
1-D array of HEALPix pixels
step : int
The number of steps to take along each edge.
nside : int
Number of pixels along the side of each of the 12 top-level HEALPix tiles
order : { 'nested' | 'ring' }
Order of HEALPix pixels
Returns
-------
lon, lat : :class:`~astropy.units.Quantity`
The longitude and latitude, as 2-D arrays where the first dimension is
the same as the ``healpix_index`` input, and the second dimension has
size ``4 * step``.
"""
healpix_index = np.asarray(healpix_index, dtype=np.int64)
step = int(step)
if step < 1:
raise ValueError('step must be at least 1')
# PERF: this could be optimized by writing a Cython routine to do this to
# avoid allocating temporary arrays
frac = np.linspace(0., 1., step + 1)[:-1]
dx = np.hstack([1 - frac, np.repeat(0, step), frac, np.repeat(1, step)])
dy = np.hstack([np.repeat(1, step), 1 - frac, np.repeat(0, step), frac])
healpix_index, dx, dy = np.broadcast_arrays(healpix_index.reshape(-1, 1), dx, dy)
lon, lat = healpix_to_lonlat(healpix_index.ravel(), nside, dx.ravel(), dy.ravel(), order=order)
lon = lon.reshape(-1, 4 * step)
lat = lat.reshape(-1, 4 * step)
return lon, lat