_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 HEALPix `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 HEALPix `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.i02hash(np.int32(j)) << 1) | self.i02hash(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_HASH[np.uint32(i) >> np.uint32(24)] << np.int64(48))
        val2 = np.int64(self.LUPT_TO_HASH[(np.uint32(i) & 0x00FF0000) >> np.uint32(16)] << np.uint64(32))
        val3 = np.int64(self.LUPT_TO_HASH[(np.uint32(i) & 0x0000FF00) >> np.uint32(8)] << np.uint64(16))
        val4 = np.int64(self.LUPT_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_INT[
            np.int32((np.uint64(h)
                      & np.uint64(0xFF00000000000000)) >> np.uint64(56))] << np.int64(28)
        val2 = self.LUPT_TO_IJ_INT[
            np.int32((np.uint64(h)
                      & np.uint64(0x00FF000000000000)) >> np.uint64(48))] << np.int64(24)
        val3 = self.LUPT_TO_IJ_INT[
            np.int32((np.uint64(h)
                      & np.uint64(0x0000FF0000000000)) >> np.uint64(40))] << np.int64(20)
        val4 = self.LUPT_TO_IJ_INT[
            np.int32((np.uint64(h)
                      & np.uint64(0x000000FF00000000)) >> np.uint64(32))] << np.int64(16)
        val5 = self.LUPT_TO_IJ_INT[
            np.int32((np.uint64(h)
                      & np.uint64(0x00000000FF000000)) >> np.uint64(24))] << np.int64(12)
        val6 = self.LUPT_TO_IJ_INT[
            np.int32((np.uint64(h)
                      & np.uint64(0x0000000000FF0000)) >> np.uint64(16))] << np.int64(8)
        val7 = self.LUPT_TO_IJ_INT[
            np.int32((np.uint64(h)
                      & np.uint64(0x000000000000FF00)) >> np.uint64(8))] << np.int64(4)
        val8 = self.LUPT_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.hash2ij(_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))