_hpxUtils.py

Go to the documentation of this file.

# This file is part of afw.
#
# Developed for the LSST Data Management System.
# This product includes software developed by the LSST Project
# (https://www.lsst.org).
# See the COPYRIGHT file at the top-level directory of this distribution
# for details of code ownership.
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <https://www.gnu.org/licenses/>.
#
# This code is partly based on AladinSrc.jar version 11.024.
# AladinSrc.jar is licensed with GPLv3, see http://aladin.u-strasbg.fr/COPYING
import numpy as np
 
from lsst.daf.base import PropertySet
 
from ._geom import makeSkyWcs
 
 
def makeHpxWcs(hips_order, pixel, shift_order=9):
    """
    Make a SkyWcs object with HEALPix grid projection (HPX).
 
    The SkyWcs generated by this function is suitable to be used with a
    Hierarchical Progressive Survey (HiPS) FITS image as described in
    https://www.ivoa.net/documents/HiPS/20170519/REC-HIPS-1.0-20170519.pdf
 
    A HiPS image covers one HEALPix cell, with the HEALPix nside equal to
    2**hips_order. Each cell is 'shift_order' orders deeper than the HEALPix
    cell, with 2**shift_order x 2**shift_order sub-pixels on a side, which
    defines the target resolution of the HiPS image. The IVOA recommends
    shift_order=9, for 2**9=512 pixels on a side.  See Notes below to
    convert from hips_order to image resolution.
 
    Parameters
    ----------
    hips_order : `int`
        HiPS order, such that HEALPix nside=2**hips_order.
        Must be a positive integer.
    pixel : `int`
        Pixel number in the nest ordering scheme.
    shift_order : `int`, optional
        Shift order for subpixels, such that there are 2**shift_order
        sub-pixels on a side of the HiPS cell.
        Must be a positive integer.
 
    Returns
    -------
    wcs : `lsst.geom.SkyWcs`
 
    Raises
    ------
    `ValueError`: Raise if hips_order is <=0, or if shift_order is <=0, or
        if pixel number is out of range for the given hips_order
        (0 <= pixel < 12*nside*nside).
 
    Notes
    -----
    Table 5 from
    https://www.ivoa.net/documents/HiPS/20170519/REC-HIPS-1.0-20170519.pdf
    shows the relationship between hips_order, number of tiles (full
    sky coverage), cell size, and sub-pixel size/image resolution (with
    the default shift_order=9):
 
    +------------+-----------------+--------------+------------------+
    | hips_order | Number of Tiles | Cell Size    | Image Resolution |
    +============+=================+==============+==================+
    | 0          | 12              | 58.63 deg    | 6.871 arcmin     |
    | 1          | 48              | 29.32 deg    | 3.435 arcmin     |
    | 2          | 192             | 14.66 deg    | 1.718 arcmin     |
    | 3          | 768             | 7.329 deg    | 51.53 arcsec     |
    | 4          | 3072            | 3.665 deg    | 25.77 arcsec     |
    | 5          | 12288           | 1.832 deg    | 12.88 arcsec     |
    | 6          | 49152           | 54.97 arcmin | 6.442 arcsec     |
    | 7          | 196608          | 27.48 arcmin | 3.221 arcsec     |
    | 8          | 786432          | 13.74 arcmin | 1.61 arcsec      |
    | 9          | 3145728         | 6.871 arcmin | 805.2mas         |
    | 10         | 12582912        | 3.435 arcmin | 402.6mas         |
    | 11         | 50331648        | 1.718 arcmin | 201.3mas         |
    | 12         | 201326592       | 51.53 arcsec | 100.6mas         |
    | 13         | 805306368       | 25.77 arcsec | 50.32mas         |
    +------------+-----------------+--------------+------------------+
    """
    if shift_order <= 0:
        raise ValueError(f"shift_order {shift_order} must be positive.")
    hips_tilepix = 2**shift_order
 
    if hips_order <= 0:
        raise ValueError(f"order {hips_order} must be positive.")
    nside_cell = 2**hips_order
 
    if pixel < 0 or pixel >= 12*nside_cell*nside_cell:
        raise ValueError(f"pixel value {pixel} out of range.")
 
    # The HEALPix grid projection (HPX) is defined in the FITS standard
    # https://fits.gsfc.nasa.gov/standard40/fits_standard40aa-le.pdf
    # from Calabretta & Roukema (2007)
    # https://ui.adsabs.harvard.edu/abs/2007MNRAS.381..865C/abstract
    # which defines the standard H = 4, K = 3 pixelization parameters
    # encoded in PV2_1 = H, PV2_2 = K.
    # The CRVAL1, CRVAL2 values should always be 0, 0 according to
    # the FITS standard.
    # The CD matrix is defined in wcslib HPXcvt.c.
    # The Calabretta & Roukema (2007) paper and wcslib HPXcvt.c only
    # define full-sky HPX projections.  For single pixels we
    # use the code from AladinSrc.jar Tile2HPX.java to compute
    # CRPIX1, CRPIX2.
 
    # The nside of the sub-pixels is the product of the tile nside
    # and the number of sub-pixels on a side.
    nside_pix = nside_cell*hips_tilepix
    # All tiles are rotated 45 degrees.
    cos45 = np.sqrt(2.0)/2.0
    # This defines the pixel scale.
    scale = 90.0/nside_pix/np.sqrt(2.0)
    cos45_scale = cos45*scale
    # The projected center of the pixel used for the HPX header is
    # a non-trivial computation, and typically is outside of the
    # tile pixel itself.  Therefore, these values are not the same
    # as the values computed from healpy.pix2ang().
    cent_ra_proj, cent_dec_proj = _hpx_projected_center(hips_order, pixel)
 
    md = PropertySet()
    md['CD1_1'] = -cos45_scale
    md['CD1_2'] = -cos45_scale
    md['CD2_1'] = cos45_scale
    md['CD2_2'] = -cos45_scale
    md['CTYPE1'] = 'RA---HPX'
    md['CTYPE2'] = 'DEC--HPX'
    md['CRVAL1'] = 0.0
    md['CRVAL2'] = 0.0
    md['PV2_1'] = 4
    md['PV2_2'] = 3
    md['CRPIX1'] = ((hips_tilepix + 1)/2.0) - 0.5*(-cent_ra_proj/cos45_scale + cent_dec_proj/cos45_scale)
    md['CRPIX2'] = ((hips_tilepix + 1)/2.0) - 0.5*(-cent_ra_proj/cos45_scale - cent_dec_proj/cos45_scale)
 
    return makeSkyWcs(md)
 
 
def _hpx_projected_center(hips_order, pixel):
    """
    Compute the projected center for use in HPX WCS headers.
 
    The values of cent_ra_proj, cent_dec_proj computed by this function are
    typically outside of the cell pixel itself, and are not the same as
    the values computed from 'healpy.pix2ang()'.
 
    Code is adapted from AladinSrc.jar version 11.024, Tile2HPX.java.
    AladinSrc.jar is licensed with GPLv3, see
    http://aladin.u-strasbg.fr/COPYING
 
    Parameters
    ----------
    hips_order : `int`
        HiPS order, such that HEALPix nside=2**hips_order.
    pixel : `int`
        Pixel number in the nest ordering scheme.
 
    Returns
    -------
    cent_ra_proj, cent_dec_proj : `float`
        Projected center ra/dec in degrees.
 
    Raises
    ------
    `ValueError`: Raised if hips_order is <=0, or if pixel number is out of
        range for the given order (0 < 12*nside*nside).
    """
    if hips_order <= 0:
        raise ValueError(f"hips_order {hips_order} must be positive.")
    nside_cell = 2**hips_order
 
    if pixel < 0 or pixel >= 12*nside_cell*nside_cell:
        raise ValueError(f"pixel value {pixel} out of range.")
 
    twice_depth = np.left_shift(hips_order, 1)
    xy_mask = np.left_shift(1, twice_depth) - 1
    fc = _ZOrderCurve2DInt()
 
    d0h = np.int32(np.right_shift(pixel, twice_depth))
    _hash = fc.hash2ij(pixel & xy_mask)
    i_in_d0h = fc.ij2i(_hash)
    j_in_d0h = fc.ij2j(_hash)
    # Compute coordinates from the center of the base pixel
    # with x-axis = W-->E, y-axis = S-->N
    l_in_d0h = i_in_d0h - j_in_d0h
    h_in_d0h = i_in_d0h + j_in_d0h - (nside_cell - 1)
    # Compute coordinates of the base pixel in the projection plane
    d0h_by_4_quotient = np.right_shift(d0h, 2)
    d0h_mod_4 = d0h - np.left_shift(d0h_by_4_quotient, 2)
    h_d0h = 1 - d0h_by_4_quotient
    l_d0h = np.left_shift(d0h_mod_4, 1)
    if ((h_d0h == 0) and ((l_d0h == 6) or ((l_d0h == 4) and (l_in_d0h > 0)))):
        # Equatorial region
        l_d0h -= 8
    elif (h_d0h != 0):
        # Polar caps regions
        l_d0h += 1
        if (l_d0h > 3):
            l_d0h -= 8
    # Finalize
    return (np.rad2deg((np.pi/4.)*(l_d0h + l_in_d0h/float(nside_cell))),
            np.rad2deg((np.pi/4.)*(h_d0h + h_in_d0h/float(nside_cell))))
 
 
class _ZOrderCurve2DInt(object):
    """
    Z-Order 2D curve for 32-bit integer values.
 
    Code is ported from AladinSrc.jar version 11.024,
    ZOrderCurve2DImpls.java.
    AladinSrc.jar is licensed with GPLv3, see
    http://aladin.u-strasbg.fr/COPYING
 
    From the original documentation:
    "Z-Order Curve (ZOC) implementation in which the vertical coordinate
    carry the most significant bit (VMSB). This implementation is based
    on a lookup table (LOOKUP). We assume that each discritized
    coordinates is coded on maximum 32 bits (INT)."
    """
    LUPT_TO_HASH = np.array([
        0x0000, 0x0001, 0x0004, 0x0005, 0x0010, 0x0011, 0x0014, 0x0015, 0x0040, 0x0041, 0x0044,
        0x0045, 0x0050, 0x0051, 0x0054, 0x0055, 0x0100, 0x0101, 0x0104, 0x0105, 0x0110, 0x0111,
        0x0114, 0x0115, 0x0140, 0x0141, 0x0144, 0x0145, 0x0150, 0x0151, 0x0154, 0x0155, 0x0400,
        0x0401, 0x0404, 0x0405, 0x0410, 0x0411, 0x0414, 0x0415, 0x0440, 0x0441, 0x0444, 0x0445,
        0x0450, 0x0451, 0x0454, 0x0455, 0x0500, 0x0501, 0x0504, 0x0505, 0x0510, 0x0511, 0x0514,
        0x0515, 0x0540, 0x0541, 0x0544, 0x0545, 0x0550, 0x0551, 0x0554, 0x0555, 0x1000, 0x1001,
        0x1004, 0x1005, 0x1010, 0x1011, 0x1014, 0x1015, 0x1040, 0x1041, 0x1044, 0x1045, 0x1050,
        0x1051, 0x1054, 0x1055, 0x1100, 0x1101, 0x1104, 0x1105, 0x1110, 0x1111, 0x1114, 0x1115,
        0x1140, 0x1141, 0x1144, 0x1145, 0x1150, 0x1151, 0x1154, 0x1155, 0x1400, 0x1401, 0x1404,
        0x1405, 0x1410, 0x1411, 0x1414, 0x1415, 0x1440, 0x1441, 0x1444, 0x1445, 0x1450, 0x1451,
        0x1454, 0x1455, 0x1500, 0x1501, 0x1504, 0x1505, 0x1510, 0x1511, 0x1514, 0x1515, 0x1540,
        0x1541, 0x1544, 0x1545, 0x1550, 0x1551, 0x1554, 0x1555, 0x4000, 0x4001, 0x4004, 0x4005,
        0x4010, 0x4011, 0x4014, 0x4015, 0x4040, 0x4041, 0x4044, 0x4045, 0x4050, 0x4051, 0x4054,
        0x4055, 0x4100, 0x4101, 0x4104, 0x4105, 0x4110, 0x4111, 0x4114, 0x4115, 0x4140, 0x4141,
        0x4144, 0x4145, 0x4150, 0x4151, 0x4154, 0x4155, 0x4400, 0x4401, 0x4404, 0x4405, 0x4410,
        0x4411, 0x4414, 0x4415, 0x4440, 0x4441, 0x4444, 0x4445, 0x4450, 0x4451, 0x4454, 0x4455,
        0x4500, 0x4501, 0x4504, 0x4505, 0x4510, 0x4511, 0x4514, 0x4515, 0x4540, 0x4541, 0x4544,
        0x4545, 0x4550, 0x4551, 0x4554, 0x4555, 0x5000, 0x5001, 0x5004, 0x5005, 0x5010, 0x5011,
        0x5014, 0x5015, 0x5040, 0x5041, 0x5044, 0x5045, 0x5050, 0x5051, 0x5054, 0x5055, 0x5100,
        0x5101, 0x5104, 0x5105, 0x5110, 0x5111, 0x5114, 0x5115, 0x5140, 0x5141, 0x5144, 0x5145,
        0x5150, 0x5151, 0x5154, 0x5155, 0x5400, 0x5401, 0x5404, 0x5405, 0x5410, 0x5411, 0x5414,
        0x5415, 0x5440, 0x5441, 0x5444, 0x5445, 0x5450, 0x5451, 0x5454, 0x5455, 0x5500, 0x5501,
        0x5504, 0x5505, 0x5510, 0x5511, 0x5514, 0x5515, 0x5540, 0x5541, 0x5544, 0x5545, 0x5550,
        0x5551, 0x5554, 0x5555], dtype=np.int16)
 
    LUPT_TO_IJ_INT = np.array([
        0x000000000, 0x000000001, 0x100000000, 0x100000001, 0x000000002, 0x000000003,
        0x100000002, 0x100000003, 0x200000000, 0x200000001, 0x300000000, 0x300000001,
        0x200000002, 0x200000003, 0x300000002, 0x300000003, 0x000000004, 0x000000005,
        0x100000004, 0x100000005, 0x000000006, 0x000000007, 0x100000006, 0x100000007,
        0x200000004, 0x200000005, 0x300000004, 0x300000005, 0x200000006, 0x200000007,
        0x300000006, 0x300000007, 0x400000000, 0x400000001, 0x500000000, 0x500000001,
        0x400000002, 0x400000003, 0x500000002, 0x500000003, 0x600000000, 0x600000001,
        0x700000000, 0x700000001, 0x600000002, 0x600000003, 0x700000002, 0x700000003,
        0x400000004, 0x400000005, 0x500000004, 0x500000005, 0x400000006, 0x400000007,
        0x500000006, 0x500000007, 0x600000004, 0x600000005, 0x700000004, 0x700000005,
        0x600000006, 0x600000007, 0x700000006, 0x700000007, 0x000000008, 0x000000009,
        0x100000008, 0x100000009, 0x00000000A, 0x00000000B, 0x10000000A, 0x10000000B,
        0x200000008, 0x200000009, 0x300000008, 0x300000009, 0x20000000A, 0x20000000B,
        0x30000000A, 0x30000000B, 0x00000000C, 0x00000000D, 0x10000000C, 0x10000000D,
        0x00000000E, 0x00000000F, 0x10000000E, 0x10000000F, 0x20000000C, 0x20000000D,
        0x30000000C, 0x30000000D, 0x20000000E, 0x20000000F, 0x30000000E, 0x30000000F,
        0x400000008, 0x400000009, 0x500000008, 0x500000009, 0x40000000A, 0x40000000B,
        0x50000000A, 0x50000000B, 0x600000008, 0x600000009, 0x700000008, 0x700000009,
        0x60000000A, 0x60000000B, 0x70000000A, 0x70000000B, 0x40000000C, 0x40000000D,
        0x50000000C, 0x50000000D, 0x40000000E, 0x40000000F, 0x50000000E, 0x50000000F,
        0x60000000C, 0x60000000D, 0x70000000C, 0x70000000D, 0x60000000E, 0x60000000F,
        0x70000000E, 0x70000000F, 0x800000000, 0x800000001, 0x900000000, 0x900000001,
        0x800000002, 0x800000003, 0x900000002, 0x900000003, 0xA00000000, 0xA00000001,
        0xB00000000, 0xB00000001, 0xA00000002, 0xA00000003, 0xB00000002, 0xB00000003,
        0x800000004, 0x800000005, 0x900000004, 0x900000005, 0x800000006, 0x800000007,
        0x900000006, 0x900000007, 0xA00000004, 0xA00000005, 0xB00000004, 0xB00000005,
        0xA00000006, 0xA00000007, 0xB00000006, 0xB00000007, 0xC00000000, 0xC00000001,
        0xD00000000, 0xD00000001, 0xC00000002, 0xC00000003, 0xD00000002, 0xD00000003,
        0xE00000000, 0xE00000001, 0xF00000000, 0xF00000001, 0xE00000002, 0xE00000003,
        0xF00000002, 0xF00000003, 0xC00000004, 0xC00000005, 0xD00000004, 0xD00000005,
        0xC00000006, 0xC00000007, 0xD00000006, 0xD00000007, 0xE00000004, 0xE00000005,
        0xF00000004, 0xF00000005, 0xE00000006, 0xE00000007, 0xF00000006, 0xF00000007,
        0x800000008, 0x800000009, 0x900000008, 0x900000009, 0x80000000A, 0x80000000B,
        0x90000000A, 0x90000000B, 0xA00000008, 0xA00000009, 0xB00000008, 0xB00000009,
        0xA0000000A, 0xA0000000B, 0xB0000000A, 0xB0000000B, 0x80000000C, 0x80000000D,
        0x90000000C, 0x90000000D, 0x80000000E, 0x80000000F, 0x90000000E, 0x90000000F,
        0xA0000000C, 0xA0000000D, 0xB0000000C, 0xB0000000D, 0xA0000000E, 0xA0000000F,
        0xB0000000E, 0xB0000000F, 0xC00000008, 0xC00000009, 0xD00000008, 0xD00000009,
        0xC0000000A, 0xC0000000B, 0xD0000000A, 0xD0000000B, 0xE00000008, 0xE00000009,
        0xF00000008, 0xF00000009, 0xE0000000A, 0xE0000000B, 0xF0000000A, 0xF0000000B,
        0xC0000000C, 0xC0000000D, 0xD0000000C, 0xD0000000D, 0xC0000000E, 0xC0000000F,
        0xD0000000E, 0xD0000000F, 0xE0000000C, 0xE0000000D, 0xF0000000C, 0xF0000000D,
        0xE0000000E, 0xE0000000F, 0xF0000000E, 0xF0000000F], dtype=np.int64)
 
    def __init__(self):
        pass
 
    def xy2hash(self, x, y):
        """
        Compute the hash value from x/y.
 
        Parameters
        ----------
        x : `float`
            x coordinate along the horizontal axis.
            Must fit within the 32-bit integer range.
        y : `float`
            y coordinate along the vertical axis.
            Must fit within the 32-bit integer range.
 
        Returns
        -------
        hash : `numpy.int64`
            The space-filling hash value associated with the
            given coordinates.
        """
        return self.ij2hash(np.int32(x), np.int32(y))
 
    def ij2hash(self, i, j):
        """
        Compute the hash value from discretized i, j.
 
        Parameters
        ----------
        i : `int`
            i discretized coordinate along the horizontal axis.
            Must fit within the 32-bit integer range.
        j : `int`
            j discretized coordinate along the vertical axis.
            Must fit within the 32-bit integer range.
 
        Returns
        -------
        hash : `numpy.int64`
            The space-filling hash value associated with the
            given coordinates.
        """
        return (self.i02hashi02hash(np.int32(j)) << 1) | self.i02hashi02hash(np.int32(i))
 
    def i02hash(self, i):
        """
        Special case of ij2hash in which the discretized coordinate along
        the vertical axis equals zero.
 
        Parameters
        ----------
        i : `int`
            i discretized coordinate along the horizontal axis.
            Must fit within the 32-bit integer range.
 
        Returns
        -------
        hash : `numpy.int64`
            The space-filling hash value associated with the
            given coordinate.
        """
        val1 = np.int64(self.LUPT_TO_HASHLUPT_TO_HASH[np.uint32(i) >> np.uint32(24)] << np.int64(48))
        val2 = np.int64(self.LUPT_TO_HASHLUPT_TO_HASH[(np.uint32(i) & 0x00FF0000) >> np.uint32(16)] << np.uint64(32))
        val3 = np.int64(self.LUPT_TO_HASHLUPT_TO_HASH[(np.uint32(i) & 0x0000FF00) >> np.uint32(8)] << np.uint64(16))
        val4 = np.int64(self.LUPT_TO_HASHLUPT_TO_HASH[np.uint32(i) & 0x000000FF])
        return val1 | val2 | val3 | val4
 
    def hash2ij(self, h):
        """
        Transforms the given space-filling hash value into a single value
        from which the 2d coordinates can be extracted using ij2i and ij2j.
 
        Parameters
        ----------
        h : `int`
            Space-filling hash value
 
        Returns
        -------
        ij : `np.int64`
            Single value from which 2d coordinates can be extracted.
        """
        val1 = self.LUPT_TO_IJ_INTLUPT_TO_IJ_INT[
            np.int32((np.uint64(h)
                      & np.uint64(0xFF00000000000000)) >> np.uint64(56))] << np.int64(28)
        val2 = self.LUPT_TO_IJ_INTLUPT_TO_IJ_INT[
            np.int32((np.uint64(h)
                      & np.uint64(0x00FF000000000000)) >> np.uint64(48))] << np.int64(24)
        val3 = self.LUPT_TO_IJ_INTLUPT_TO_IJ_INT[
            np.int32((np.uint64(h)
                      & np.uint64(0x0000FF0000000000)) >> np.uint64(40))] << np.int64(20)
        val4 = self.LUPT_TO_IJ_INTLUPT_TO_IJ_INT[
            np.int32((np.uint64(h)
                      & np.uint64(0x000000FF00000000)) >> np.uint64(32))] << np.int64(16)
        val5 = self.LUPT_TO_IJ_INTLUPT_TO_IJ_INT[
            np.int32((np.uint64(h)
                      & np.uint64(0x00000000FF000000)) >> np.uint64(24))] << np.int64(12)
        val6 = self.LUPT_TO_IJ_INTLUPT_TO_IJ_INT[
            np.int32((np.uint64(h)
                      & np.uint64(0x0000000000FF0000)) >> np.uint64(16))] << np.int64(8)
        val7 = self.LUPT_TO_IJ_INTLUPT_TO_IJ_INT[
            np.int32((np.uint64(h)
                      & np.uint64(0x000000000000FF00)) >> np.uint64(8))] << np.int64(4)
        val8 = self.LUPT_TO_IJ_INTLUPT_TO_IJ_INT[
            np.int32((np.uint64(h)
                      & np.uint64(0x00000000000000FF)))]
        return val1 | val2 | val3 | val4 | val5 | val6 | val7 | val8
 
    def hash2i0(self, _hash):
        """
        Special case of hash2ij in which the discretized coordinate along
        the vertical axis is zero.
 
        Parameters
        ----------
        _hash : `int`
            Space-filling hash value.
 
        Returns
        -------
        ij : `np.int64`
            Single value from which 2d coordinates can be extracted.
        """
        assert((0xFFFFFFFF33333333 & np.int64(_hash)) == 0)
        return self.hash2ijhash2ij(_hash)
 
    def ij2i(self, ij):
        """
        Extract the discretized horizontal coordinate from hash2ij.
 
        Parameters
        ----------
        ij : `int`
            The ij result of hash2ij.
 
        Returns
        -------
        i : `np.int32`
            Discretized horizontal coordinate stored in ij.
        """
        return np.int32(ij)
 
    def ij2j(self, ij):
        """
        Extract the discretized vertical coordinate from hash2ij.
 
        Parameters
        ----------
        ij : `int`
            The ij result of hash2ij.
 
        Returns
        -------
        j : `np.int32`
            Discretized vertical coordinate stored in ij.
        """
        return np.int32(np.uint64(ij) >> np.uint64(32))