Add all-sky image class and major update to existing tile code

hips/draw/simple.py

@@ -1,10 +1,11 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 """HiPS tile drawing -- simple method."""
 import numpy as np
-from typing import Tuple, List
+from typing import List
 from astropy.wcs.utils import proj_plane_pixel_scales
 from skimage.transform import ProjectiveTransform, warp
 from ..tiles import HipsSurveyProperties, HipsTile, HipsTileMeta
+from ..tiles.tile import compute_image_shape
 from ..utils import WCSGeometry, healpix_pixels_in_sky_image, hips_order_for_pixel_resolution
 
 __all__ = [
@@ -124,29 +125,6 @@ def tiles(self) -> List[HipsTile]:
 
         return self._tiles
 
-    @property
-    def shape(self) -> Tuple[int]:
-        """Shape of the output image.
-
-        The output image shape is 2-dim for grayscale, and 3-dim for color images:
-
-        * ``shape = (height, width)`` for FITS images with one grayscale channel
-        * ``shape = (height, width, 3)`` for JPG images with three RGB channels
-        * ``shape = (height, width, 4)`` for PNG images with four RGBA channels
-        """
-        height = self.geometry.shape.height
-        width = self.geometry.shape.width
-        tile_format = self.tile_format
-
-        if tile_format == 'fits':
-            return height, width
-        elif tile_format == 'jpg':
-            return height, width, 3
-        elif tile_format == 'png':
-            return height, width, 4
-        else:
-            return ValueError(f'Invalid tile format: {tile_format}')
-
     def warp_image(self, tile: HipsTile) -> np.ndarray:
         """Warp a HiPS tile and a sky image."""
         return warp(
@@ -160,7 +138,12 @@ def draw_tiles(self) -> np.ndarray:
         """Draw HiPS tiles onto an empty image."""
         tiles = self.tiles
 
-        image = np.zeros(self.shape, dtype=np.float32)
+        shape = compute_image_shape(
+            width=self.geometry.shape.width,
+            height=self.geometry.shape.height,
+            fmt=self.tile_format
+        )
+        image = np.zeros(shape, dtype=np.float32)
         for tile in tiles:
             tile_image = self.warp_image(tile)
             # TODO: put better algorithm here instead of summing pixels

hips/draw/tests/test_simple.py

@@ -64,16 +64,13 @@ def setup_class(cls):
             width=2000, height=1000, fov="3 deg",
             coordsys='icrs', projection='AIT',
         )
-        cls.simple_tile_painter = SimpleTilePainter(cls.geometry, cls.hips_survey, 'fits')
+        cls.painter = SimpleTilePainter(cls.geometry, cls.hips_survey, 'fits')
 
     def test_draw_hips_order(self):
-        assert self.simple_tile_painter.draw_hips_order == 7
-
-    def test_shape(self):
-        assert self.simple_tile_painter.shape == (1000, 2000)
+        assert self.painter.draw_hips_order == 7
 
     def test_tile_indices(self):
-        assert list(self.simple_tile_painter.tile_indices)[:4] == [69623, 69627, 69628, 69629]
+        assert list(self.painter.tile_indices)[:4] == [69623, 69627, 69628, 69629]
 
     draw_hips_order_pars = [
         dict(order=7, fov="3 deg"),
@@ -93,13 +90,14 @@ def test_compute_matching_hips_order(self, pars):
         assert simple_tile_painter.draw_hips_order == pars['order']
 
     def test_run(self):
-        self.simple_tile_painter.run()
-        assert_allclose(self.simple_tile_painter.image[200, 994], 2120)
+        self.painter.run()
+        assert self.painter.image.shape == (1000, 2000)
+        assert_allclose(self.painter.image[200, 994], 2120)
 
     def test_draw_hips_tile_grid(self):
-        self.simple_tile_painter.plot_mpl_hips_tile_grid()
+        self.painter.plot_mpl_hips_tile_grid()
 
     def test_draw_debug_image(self):
-        tile = self.simple_tile_painter.tiles[3]
-        image = self.simple_tile_painter.image
+        tile = self.painter.tiles[3]
+        image = self.painter.image
         plot_mpl_single_tile(self.geometry, tile, image)

hips/tiles/__init__.py

@@ -2,3 +2,4 @@
 """Classes and functions to manage HiPS tiles."""
 from .tile import *
 from .surveys import *
+from .allsky import *

hips/tiles/allsky.py

@@ -0,0 +1,163 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+from typing import List
+import numpy as np
+from ..utils.healpix import healpix_order_to_npix
+from ..tiles import HipsTile
+from ..tiles.tile import compute_image_shape
+
+__all__ = [
+    'HipsTileAllskyArray',
+]
+
+__doctest_skip__ = [
+    'HipsTileAllskyArray',
+]
+
+
+class HipsTileAllskyArray(HipsTile):
+    """All-sky tile array container.
+
+    To quote from section 4.3.2 "Allsky preview file" of the HiPS IVOA working draft:
+    "The tiles at low orders (0 to 3) may be packaged
+    together into a unique file called Allsky."
+
+    This class implements that all-sky tile array format.
+
+    TODO: We're sub-classing `~hips.HipsTile` here at the moment.
+    This is weird!
+    Probably the current ``HipsTile`` should be renamed ``ImageIO``
+    or be split up into functions that do image I/O in the three supported formats?
+
+    TODO: We re-use the `~hips.HipsTileMeta` class to store ``order``
+    as well as other info like ``file_format`` and ``frame``.
+    Note that ``ipix`` doesn't make sense for an ``AllSkyTileArray``.
+    Maybe there's  a better way to handle this without code duplication?
+
+    Examples
+    --------
+
+    Load an example existing HiPS all-sky image
+    (unfortunately one has to pass a dummy ipix value here):
+
+    >>> from hips import HipsTileAllskyArray, HipsTileMeta
+    >>> meta = HipsTileMeta(order=3, ipix=-1, file_format='jpg', frame='icrs')
+    >>> url = 'http://alasky.unistra.fr/Fermi/Color/Norder3/Allsky.jpg'
+    >>> allsky = HipsTileAllskyArray.fetch(meta, url)
+
+    Now you can extract tiles (e.g. for drawing):
+
+    >>> tile = allsky.tile(ipix=42)
+    >>> tile.meta
+    HipsTileMeta(order=3, ipix=42, file_format='jpg', frame='icrs', width=64)
+
+    TODO: add an example how to go the other way: combine tiles into an allsky image.
+    """
+
+    def __repr__(self):
+        return (
+            'HipsTileAllskyArray('
+            f'format={self.meta.file_format!r}, order={self.meta.order}, '
+            f'width={self.width}, height={self.height}, '
+            f'n_tiles={self.n_tiles}, '
+            f'tile_width={self.tile_width}'
+            ')'
+        )
+
+    @property
+    def width(self) -> int:
+        """Image pixel width (int)."""
+        return self.data.shape[1]
+
+    @property
+    def height(self) -> int:
+        """Image pixel height (int)"""
+        return self.data.shape[0]
+
+    @property
+    def n_tiles(self) -> int:
+        """Number of tiles in the image (int)."""
+        return healpix_order_to_npix(self.meta.order)
+
+    @property
+    def n_tiles_in_row(self) -> int:
+        """Number of tiles per tile row (int)."""
+        return int(np.sqrt(self.n_tiles))
+
+    @property
+    def tile_width(self) -> int:
+        """Pixel width of a single tile (int)."""
+        return self.width // self.n_tiles_in_row
+
+    @classmethod
+    def from_tiles(cls, tiles: List[HipsTile]) -> 'HipsTileAllskyArray':
+        """Create all-sky image from list of tiles."""
+        meta = tiles[0].meta.copy()
+        data = cls.tiles_to_allsky_array(tiles)
+        # TODO: check return type here.
+        # Pycharm warns that a `HipsTile` is returned here, not a `HipsTileAllskyArray`
+        # Is this true or a bug in their static code analysis?
+        return cls.from_numpy(meta, data)
+
+    @staticmethod
+    def tiles_to_allsky_array(tiles: List[HipsTile]) -> np.ndarray:
+        """Combine tiles into an all-sky image."""
+        # Compute all-sky image parameters that we need below
+        n_tiles = len(tiles)
+        n_tiles_in_row = int(np.sqrt(n_tiles))
+        n_tiles_in_col = (n_tiles // n_tiles_in_row) + 1
+        tile_width = tiles[0].meta.width
+
+        # Make an empty all-sky image
+        shape = compute_image_shape(
+            width=tile_width * n_tiles_in_row,
+            height=tile_width * n_tiles_in_col,
+            fmt=tiles[0].meta.file_format,
+        )
+        data = np.empty(shape, tiles[0].data.dtype)
+
+        # Copy over the tile data into the all-sky image
+        for tile in tiles:
+            tile_slice = HipsTileAllskyArray._tile_slice(
+                ipix=tile.meta.ipix,
+                tile_width=tile.meta.width,
+                n_tiles_in_row=n_tiles_in_row,
+            )
+            data[tile_slice] = tile.data
+
+        return data
+
+    @property
+    def tiles(self) -> List[HipsTile]:
+        """Split into a list of `~hips.HipsTile`.
+
+        This is called when using the all-sky image for drawing.
+        """
+        return [self.tile(ipix) for ipix in range(self.n_tiles)]
+
+    def tile(self, ipix: int, copy: bool = True) -> HipsTile:
+        """Extract one of the tiles (`~hips.HipsTile`)
+
+        A copy of the data by default.
+        For drawing we could avoid the copy by passing ``copy=False`` here.
+        """
+        meta = self.meta.copy()
+        meta.ipix = ipix
+        meta.width = self.tile_width
+
+        tile_slice = self._tile_slice(ipix, self.tile_width, self.n_tiles_in_row)
+        data = self.data[tile_slice]
+
+        if copy:
+            data = data.copy()
+
+        return HipsTile.from_numpy(meta, data)
+
+    @staticmethod
+    def _tile_slice(ipix, tile_width, n_tiles_in_row):
+        """Compute the 2-dim slice in the allsky ``data`` for a given tile."""
+        w = tile_width
+        row_idx, col_idx = divmod(ipix, n_tiles_in_row)
+        return [
+            slice(row_idx * w, (row_idx + 1) * w),
+            slice(col_idx * w, (col_idx + 1) * w),
+        ]

hips/tiles/tests/test_allsky.py

@@ -0,0 +1,127 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+from pathlib import Path
+import pytest
+from astropy.tests.helper import remote_data
+from numpy.testing import assert_equal
+from ...utils.testing import get_hips_extra_file, requires_hips_extra
+from ..tile import HipsTileMeta
+from ..allsky import HipsTileAllskyArray
+
+TEST_CASES = [
+    dict(
+        meta=dict(order=3, ipix=463, file_format='jpg'),
+        url='http://alasky.unistra.fr/Fermi/Color/Norder3/Allsky.jpg',
+        filename='datasets/samples/FermiColor/Norder3/Allsky.jpg',
+
+        repr="HipsTileAllskyArray(format='jpg', order=3, width=1728, "
+             "height=1856, n_tiles=768, tile_width=64)",
+        dtype='uint8',
+        shape=(1856, 1728, 3),
+        pix_idx=[[510], [5]],
+        pix_val=[[90, 89, 85]],
+        tile_pix_val=[49, 44, 38],
+    ),
+]
+
+
+def _read_tile(pars):
+    meta = HipsTileMeta(**pars['meta'])
+    filename = get_hips_extra_file(pars['filename'])
+    return HipsTileAllskyArray.read(meta, filename)
+
+
+@requires_hips_extra()
+@pytest.mark.parametrize('pars', TEST_CASES)
+def test_read(pars):
+    # Check that reading tiles in various formats works,
+    # i.e. that pixel data in numpy array format
+    # has correct shape, dtype and values
+    allsky = _read_tile(pars)
+
+    assert repr(allsky) == pars['repr']
+    assert allsky.meta.order == pars['meta']['order']
+    assert isinstance(allsky.raw_data, bytes)
+
+    data = allsky.data
+    assert data.shape == pars['shape']
+    assert data.dtype.name == pars['dtype']
+    assert_equal(allsky.data[pars['pix_idx']], pars['pix_val'])
+
+
+@requires_hips_extra()
+@pytest.mark.parametrize('pars', TEST_CASES)
+def test_from_numpy(pars):
+    tile = _read_tile(pars)
+    tile2 = HipsTileAllskyArray.from_numpy(tile.meta, tile.data)
+
+    # JPEG encoding is lossy. So in that case, output pixel value
+    # aren't exactly the same as input pixel values
+    if tile.meta.file_format != 'jpg':
+        assert_equal(tile.data, tile2.data)
+
+
+@remote_data
+@requires_hips_extra()
+@pytest.mark.parametrize('pars', TEST_CASES)
+def test_fetch(pars):
+    # Check that tile HTTP fetch gives the same result as tile read from disk
+    meta = HipsTileMeta(**pars['meta'])
+    tile_fetch = HipsTileAllskyArray.fetch(meta, url=pars['url'])
+
+    filename = get_hips_extra_file(pars['filename'])
+    tile_read = HipsTileAllskyArray.read(meta, filename)
+
+    assert tile_fetch == tile_read
+
+
+@requires_hips_extra()
+@pytest.mark.parametrize('pars', TEST_CASES)
+def test_write(tmpdir, pars):
+    # Check that tile I/O works, i.e. round-trips on write / read
+    tile = _read_tile(pars)
+
+    filename = str(tmpdir / Path(pars['filename']).name)
+    tile.write(filename)
+    tile2 = HipsTileAllskyArray.read(tile.meta, filename)
+
+    assert tile == tile2
+
+
+@requires_hips_extra()
+@pytest.mark.parametrize('pars', TEST_CASES)
+def test_tile(pars):
+    allsky = _read_tile(pars)
+
+    tile = allsky.tile(0)
+    print(tile.meta)
+
+    assert tile.data.shape == (64, 64, 3)
+    assert_equal(tile.data[0, 0], pars['tile_pix_val'])
+
+
+@requires_hips_extra()
+@pytest.mark.parametrize('pars', TEST_CASES)
+def test_from_tiles(pars):
+    # Check that ``from_tiles`` works properly
+    # For now, we check that ``tiles`` and ``from_tiles`` round-trip
+    # TODO: it would probably be better to test them separately,
+    # asserting on each of the two step. Round-trip can work, if
+    # the same mistake is made in each conversion step.
+    allsky = _read_tile(pars)
+    # print(allsky)
+    # allsky.write('/tmp/allsky.jpg')
+
+    tiles = allsky.tiles
+
+    # allsky2 = HipsTileAllskyArray.from_tiles(tiles)
+    # print(allsky2)
+    # allsky.write('/tmp/allsky2.jpg')
+
+    # TODO: at the moment `HipsTileAllskyArray` always does a JPG encoding,
+    # it can't hold the numpy array data unchanged.
+    # This still doesn't work, because when going ``allsky.tiles`` another
+    # JPG encoding happens.
+    # I did check the JPG files written above. They look the same, so it's working.
+    # Sigh.
+    data2 = HipsTileAllskyArray.tiles_to_allsky_array(tiles)
+    # assert_equal(allsky.data, data2)