Support RGB in HEALpix -> HiPS

docs/hipsgen.py

@@ -0,0 +1,52 @@
+"""Example how to generate HiPS from HEALPix.
+
+We will use
+"""
+import logging
+import numpy as np
+import hips
+from astropy.coordinates import SkyCoord
+from astropy_healpix import HEALPix
+
+
+def make_healpix_data(hips_order, tile_width, file_format):
+    """Silly example of HEALPix data.
+
+    Angular distance in deg to (0, 0) as pixel values.
+    """
+    nside = tile_width * (2 ** hips_order)
+    healpix = HEALPix(nside=nside)
+    ipix = np.arange(healpix.npix)
+    lon, lat = healpix.healpix_to_lonlat(ipix)
+    coord = SkyCoord(lon, lat)
+    center = SkyCoord(0, 0, unit="deg")
+    data = coord.separation(center).deg
+
+    if file_format == "fits":
+        return data
+    elif file_format == "jpg":
+        data = data.astype("uint8")
+        return np.moveaxis([data, data + 1, data + 2], 0, -1)
+    elif file_format == "png":
+        data = data.astype("uint8")
+        return np.moveaxis([data, data + 1, data + 2, data + 3], 0, -1)
+    else:
+        raise ValueError()
+
+
+if __name__ == "__main__":
+    logging.basicConfig(level=logging.INFO)
+
+    file_format = "png"
+    hips_order = 3
+    tile_width = 4
+
+    hpx_data = make_healpix_data(hips_order, tile_width, file_format)
+
+    hips.healpix_to_hips(
+        hpx_data=hpx_data,
+        hips_order=hips_order,
+        tile_width=tile_width,
+        base_path="test123",
+        file_format=file_format,
+    )

docs/hipsgen.rst

@@ -0,0 +1,72 @@
+.. include:: references.txt
+
+.. _hipsgen:
+
+*************
+Generate HiPS
+*************
+
+The `~hips.healpix_to_hips` function can be used to generate HiPS data from all-sky HEALPix images.
+We give an example below.
+
+Note that this functionality is very limited, and for most applications you will want to use the
+Java ``hipsgen`` tool or the graphical user interface to ``hipsgen`` from the Aladin desktop application.
+
+To use `~hips.healpix_to_hips` you have to pass in the all-sky HEALPix data as a Numpy array
+(so it has to completely fit into memory). For grayscale images, the HEALPix data array is
+one-dimensional, RGB images with array shape ``(npix, 3)`` can be converted to JPEG HiPS tiles,
+and RGBA (where A is the transparency channel) images with array shape ``(npix, 4)`` can be
+converted to PNG HiPS tiles.
+
+Some data is distributed directly in HEALPix format, e.g. from all-sky surveys like Planck.
+Some data from high-energy telescopes like Fermi-LAT consists of event lists can be used
+to compute HEALPix images. And WCS-based images can also be reprojected to HEALPix images
+using e.g. `reproject <https://reproject.readthedocs.io/en/stable/healpix.html>`__.
+For these use cases `~hips.healpix_to_hips` can be used for the final all-sky HEALPix
+map to HiPS conversion step, giving you full flexibility over the processing and pixel
+color scales etc from Python in the pre-processing steps.
+
+Limitations
+-----------
+
+Currently the following limitations exist:
+
+- No mosaic functionality from WCS images
+  (a lot of work to implement)
+- Only one HiPS resolution is generated at the moment
+  (not hard to write a script to make multiple resolutions)
+- No ``Allsky.fits`` generated yet
+  (shouldn't be hard to generate)
+
+Dataset
+=======
+
+bla bla
+
+Example
+=======
+
+Generate HiPS:
+
+.. code-block:: bash
+
+    $ python docs/hipsgen.py
+
+Open HiPS in Aladin Lite:
+
+.. code-block:: bash
+
+    $ cd test123
+    $ python -m http.server
+    # Then open http://localhost:8000/ in your webbrowser
+
+Or open HiPS in Aladin Desktop:
+
+.. code-block:: bash
+
+    $ TODO: how to load the folder in Aladin
+
+Here's the ``hipsgen.py`` script:
+
+.. literalinclude:: hipsgen.py
+

docs/index.rst

@@ -15,6 +15,7 @@ A Python astronomy package for HiPS : Hierarchical Progressive Surveys.
    about
    installation
    getting_started
+   hipsgen
    api
    changelog
 

hips/draw/__init__.py

@@ -3,3 +3,4 @@
 from .paint import *
 from .ui import *
 from .healpix import *
+from .hipsgen import *

hips/draw/healpix.py

@@ -1,25 +1,31 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
+from typing import List
 import logging
-from pathlib import Path
 import numpy as np
-from astropy_healpix import healpy as hp
-from ..tiles import HipsTile, HipsTileMeta, HipsSurveyProperties
+from ..tiles import HipsTile, HipsTileMeta
 from ..utils.healpix import hips_tile_healpix_ipix_array
 
-__all__ = ["healpix_to_hips_tile", "healpix_to_hips"]
+__all__ = ["healpix_to_hips_tile", "healpix_to_hips_tiles"]
 
 log = logging.getLogger(__name__)
 
 
 def healpix_to_hips_tile(
-    hpx_data: np.ndarray, tile_width: int, tile_idx: int, file_format: str, frame: str
+    hpx_data: np.ndarray,
+    hips_order: int,
+    tile_width: int,
+    tile_idx: int,
+    file_format: str,
+    frame: str = "icrs",
 ) -> HipsTile:
     """Create single HiPS tile from HEALPix data.
 
     Parameters
     ----------
     hpx_data : `~numpy.ndarray`
         Healpix data stored in the "nested" scheme.
+    hips_order : int
+        HEALPix order for the HiPS tiles
     tile_width : int
         Width of the hips tile.
     tile_idx : int
@@ -39,17 +45,38 @@ def healpix_to_hips_tile(
 
     offset_ipix = tile_idx * tile_width ** 2
     ipix = hpx_ipix + offset_ipix
-    data = hpx_data[ipix]
+
+    hpx_data = np.asarray(hpx_data)
+    if file_format == "fits":
+        if hpx_data.ndim != 1:
+            raise ValueError(
+                f"Invalid hpx_data.ndim = {hpx_data.ndim}."
+                " Must be ndim = 1 for file_format='fits'."
+            )
+        data = hpx_data[ipix]
+    elif file_format == "jpg":
+        if hpx_data.ndim != 2 or hpx_data.shape[1] != 3:
+            raise ValueError(
+                f"Invalid hpx_data.shape = {hpx_data.shape}."
+                " Must be shape = (npix, 3) to represent RGB color images."
+            )
+        data = hpx_data[ipix, :]
+    elif file_format == "png":
+        if hpx_data.ndim != 2 or hpx_data.shape[1] != 4:
+            raise ValueError(
+                f"Invalid hpx_data.shape = {hpx_data.shape}."
+                " Must be shape = (npix, 4) to represent RGBA color images."
+            )
+        data = hpx_data[ipix, :]
+    else:
+        raise ValueError(f"Invalid file_format: {file_format!r}")
 
     # np.rot90 returns a rotated view so we make a copy here
     # because the view information is lost on fits io
     data = np.rot90(data).copy()
 
-    hpx_nside = hp.npix2nside(hpx_data.size / tile_width ** 2)
-    hpx_order = int(np.log2(hpx_nside))
-
     meta = HipsTileMeta(
-        order=hpx_order,
+        order=hips_order,
         ipix=tile_idx,
         file_format=file_format,
         frame=frame,
@@ -59,50 +86,41 @@ def healpix_to_hips_tile(
     return HipsTile.from_numpy(meta=meta, data=data)
 
 
-def healpix_to_hips(hpx_data, tile_width, base_path, file_format, frame):
-    """Convert HEALPix image to HiPS.
-
-    This function writes the HiPS to disk.
-    If you don't want that, use `healpix_to_hips_tile` directly.
+def healpix_to_hips_tiles(
+    hpx_data: np.ndarray,
+    hips_order: int,
+    tile_width: int,
+    file_format: str,
+    frame: str = "icrs",
+) -> List[HipsTile]:
+    """Convert HEALPix image to HiPS tiles,
 
     Parameters
     ----------
     hpx_data : `~numpy.ndarray`
-        Healpix data stored in the "nested" scheme.
+        HEALPix data stored in the "nested" scheme.
+    hips_order : int
+        HEALPix order for the HiPS tiles
     tile_width : int
         Width of the hips tiles.
-    base_path : str or `~pathlib.Path`
-        Base path.
     file_format : {'fits', 'jpg', 'png'}
         HiPS tile file format
     frame : {'icrs', 'galactic', 'ecliptic'}
         Sky coordinate frame
-    """
-    base_path = Path(base_path)
-    base_path.mkdir(exist_ok=True, parents=True)
 
-    path = base_path / "properties"
-    log.info(f"Writing {path}")
-    HipsSurveyProperties(
-        {
-            "hips_tile_format": file_format,
-            "hips_tile_width": tile_width,
-            "hips_frame": frame,
-        }
-    ).write(path)
-
-    n_tiles = hpx_data.size // tile_width ** 2
+    Returns
+    -------
+    tiles : Generator of HipsTile
+        HiPS tiles
+    """
+    n_tiles = 12 * (2 ** hips_order)
 
     for tile_idx in range(n_tiles):
-        tile = healpix_to_hips_tile(
+        yield healpix_to_hips_tile(
             hpx_data=hpx_data,
+            hips_order=hips_order,
             tile_width=tile_width,
             tile_idx=tile_idx,
             file_format=file_format,
             frame=frame,
         )
-
-        path = base_path / tile.meta.tile_default_path
-        log.info(f"Writing {path}")
-        path.parent.mkdir(exist_ok=True, parents=True)
-        tile.write(path)