High-level function: make_sky_image

.travis.yml

@@ -33,6 +33,7 @@ env:
         - SETUP_XVFB=True
         - DEBUG=True
         - ASTROPY_USE_SYSTEM_PYTEST=1
+        - FETCH_HIPS_EXTRA=true
 
     matrix:
         # Make sure that egg_info works without dependencies
@@ -80,6 +81,12 @@ install:
     - git clone git://github.com/astropy/ci-helpers.git
     - source ci-helpers/travis/setup_conda.sh
 
+    - if $FETCH_HIPS_EXTRA; then
+          git clone https://github.com/hipspy/hips-extra.git $HOME/hips-extra;
+          export HIPS_EXTRA=${HOME}/hips-extra;
+      fi
+
+
 script:
    - $MAIN_CMD $SETUP_CMD
 

docs/getting_started.rst

@@ -1,9 +1,35 @@
 .. include:: references.txt
 
+.. doctest-skip-all
+
 .. _gs:
 
 ***************
 Getting started
 ***************
 
-Coming soon ...
+The example below shows a high level use case of the ``hips`` package.
+It fetches a HiPS tile from a remote URL and draws it on a sky image.
+Then it saves it on local disk in FITS file format.
+
+::
+
+    from astropy.io import fits
+    from astropy.coordinates import SkyCoord
+    from astropy.tests.helper import remote_data
+    from hips.utils import WCSGeometry
+    from hips.draw import make_sky_image
+    from hips.tiles import HipsSurveyProperties
+
+    geometry = WCSGeometry.create(
+         skydir=SkyCoord(0, 0, unit='deg', frame='galactic'),
+         shape=(1000, 2000), coordsys='GAL',
+         projection='AIT', cdelt=0.01, crpix=(1000, 500),
+    )
+    url = 'https://raw.githubusercontent.com/hipspy/hips-extra/master/datasets/samples/DSS2Red/properties'
+    hips_survey = HipsSurveyProperties.fetch(url)
+
+    data = make_sky_image(geometry, hips_survey)
+
+    hdu = fits.PrimaryHDU(data=data, header=geometry.fits_header)
+    hdu.writeto('my_image.fits')

hips/draw/simple.py

@@ -1,5 +1,131 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 """HiPS tile drawing -- simple method."""
+from typing import Generator, Any
+import numpy as np
+from skimage import transform as tf
+from ..tiles import HipsSurveyProperties, HipsTile, HipsTileMeta
+from ..utils import WCSGeometry, compute_healpix_pixel_indices
 
 __all__ = [
+    'draw_sky_image',
+    'make_sky_image',
+    'SimpleTilePainter'
 ]
+
+
+def draw_sky_image(geometry: WCSGeometry, tiles: Generator[HipsTile, Any, Any]) -> np.ndarray:
+    """Draw sky image using the simple and quick method.
+
+    Parameters
+    ----------
+    geometry : `~hips.utils.WCSGeometry`
+        An object of WCSGeometry
+    tiles : List[HipsTile]
+        A list of HipsTile
+
+    Returns
+    -------
+    np.ndarray
+        Returns a numpy array containing all HiPS tiles projected onto it
+    """
+    # TODO: Fix type annotation issue
+    sky_sky = np.zeros(geometry.shape)
+    for tile in tiles:
+        painter = SimpleTilePainter(geometry, tile)
+        sky_sky += painter.warp_image()
+    return sky_sky
+
+
+class SimpleTilePainter:
+    """Paint a single tile using a simple projective transformation method.
+
+    The algorithm implemented is described here: :ref:`drawing_algo`.
+
+    Parameters
+    ----------
+    geometry : `~hips.utils.WCSGeometry`
+        An object of WCSGeometry
+    tile : `HipsTile`
+       An object of HipsTile
+    """
+
+    def __init__(self, geometry: WCSGeometry, tile: HipsTile) -> None:
+        self.geometry = geometry
+        self.tile = tile
+
+    @property
+    def projection(self) -> tf.ProjectiveTransform:
+        """Estimate projective transformation on a HiPS tile"""
+        corners = self.tile.meta.skycoord_corners.to_pixel(self.geometry.wcs)
+        src = np.array(corners).T.reshape((4, 2))
+        dst = self.tile.meta.dst
+        pt = tf.ProjectiveTransform()
+        pt.estimate(src, dst)
+        return pt
+
+    def warp_image(self) -> np.ndarray:
+        """Warp a HiPS tile and a sky image"""
+        return tf.warp(
+            self.tile.data,
+            self.projection,
+            output_shape=self.geometry.shape,
+            preserve_range=True,
+        )
+
+
+def fetch_tiles(healpix_pixel_indices: np.ndarray, order: int, hips_survey: HipsSurveyProperties) -> 'HipsTile':
+    """Fetch HiPS tiles from a remote URL.
+
+    Parameters
+    ----------
+    healpix_pixel_indices : np.ndarray
+        A list of HEALPix pixel indices
+    order : int
+        Order of the HEALPix map
+    hips_survey : HipsSurveyProperties
+        An object of HipsSurveyProperties
+
+    Returns
+    -------
+    'HipsTile'
+        Returns an object of  HipsTile
+    """
+    for healpix_pixel_index in healpix_pixel_indices:
+        tile_meta = HipsTileMeta(
+            order=order,
+            ipix=healpix_pixel_index,
+            frame=hips_survey.astropy_frame,
+            file_format='fits',
+        )
+        tile = HipsTile.fetch(tile_meta, hips_survey.tile_access_url + tile_meta.filename)
+        yield tile
+
+
+def make_sky_image(geometry: WCSGeometry, hips_survey: HipsSurveyProperties) -> np.ndarray:
+    """Make sky image: fetch tiles and draw.
+
+    The example for this can be found on the :ref:`gs` page.
+
+    Parameters
+    ----------
+    geometry : `~hips.utils.WCSGeometry`
+        Geometry of the output image
+    hips_survey : `~hips.HipsSurveyProperties`
+        HiPS survey properties
+
+    Returns
+    -------
+    data : `~numpy.ndarray`
+        Output image pixels
+    """
+    healpix_pixel_indices = compute_healpix_pixel_indices(
+        wcs_geometry=geometry,
+        order=hips_survey.hips_order,
+        healpix_frame=hips_survey.astropy_frame,
+    )
+    # TODO: this isn't a good API. Will become better when we have a cache.
+    tiles = fetch_tiles(healpix_pixel_indices, hips_survey.hips_order, hips_survey)
+
+    image_data = draw_sky_image(geometry, tiles)
+
+    return image_data

hips/draw/tests/test_simple.py

@@ -1 +1,55 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
+import numpy as np
+from numpy.testing import assert_allclose
+from astropy.utils.data import get_pkg_data_filename
+from astropy.tests.helper import remote_data
+from ..simple import make_sky_image, draw_sky_image
+from ...tiles import HipsSurveyProperties, HipsTileMeta, HipsTile
+from ...utils.testing import get_hips_extra_file, make_test_wcs_geometry, requires_hips_extra
+
+
+def get_test_tiles():
+    filename = get_pkg_data_filename('../../tiles/tests/data/properties.txt')
+    hips_survey = HipsSurveyProperties.read(filename)
+
+    tile1 = HipsTile.read(
+        meta=HipsTileMeta(order=3, ipix=450, file_format='fits', frame=hips_survey.astropy_frame,
+                          tile_width=512),
+        filename=get_hips_extra_file('datasets/samples/DSS2Red/Norder3/Dir0/Npix450.fits'),
+    )
+
+    tile2 = HipsTile.read(
+        meta=HipsTileMeta(order=3, ipix=451, file_format='fits', frame=hips_survey.astropy_frame,
+                          tile_width=512),
+        filename=get_hips_extra_file('datasets/samples/DSS2Red/Norder3/Dir0/Npix451.fits'),
+    )
+
+    return [tile1, tile2]
+
+
+@requires_hips_extra()
+def test_draw_sky_image():
+    geometry = make_test_wcs_geometry(case=2)
+    tiles = get_test_tiles()
+    data = draw_sky_image(geometry, tiles)
+
+    assert data.shape == geometry.shape
+    assert data.dtype == np.float64
+    assert_allclose(np.sum(data), 4575235421.5126467)
+    assert_allclose(data[400, 500], 2866.0101409848185)
+    assert_allclose(data[400, 501], 2563.6916727348043)
+
+
+@remote_data
+def test_make_sky_image():
+    url = 'https://raw.githubusercontent.com/hipspy/hips-extra/master/datasets/samples/DSS2Red/properties'
+    hips_survey = HipsSurveyProperties.fetch(url)
+    geometry = make_test_wcs_geometry(case=2)
+    data = make_sky_image(geometry, hips_survey)
+    assert data.shape == geometry.shape
+    assert data.dtype == np.float64
+    assert_allclose(data[200, 994], 3717.10091363)
+    assert_allclose(data[200, 995], 3402.55292158)
+
+
+# TODO: add tests for SimpleTilePainter with asserts on the intermediate computed things.

hips/tiles/__init__.py

@@ -1,6 +1,5 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 """Classes and functions to manage HiPS tiles."""
 from .tile import *
-from .tile_meta import *
 from .description import *
 from .surveys import *

hips/tiles/description.py

@@ -24,6 +24,12 @@ class HipsSurveyProperties:
     >>> hips_survey_property.base_url
     'http://alasky.u-strasbg.fr/DSS/DSSColor'
     """
+    hips_to_astropy_frame_mapping = OrderedDict([
+        ('equatorial', 'icrs'),
+        ('galactic', 'galactic'),
+        ('ecliptic', 'ecliptic'),
+    ])
+    """HIPS to Astropy SkyCoord frame string mapping."""
 
     def __init__(self, data: OrderedDict) -> None:
         self.data = data
@@ -99,6 +105,11 @@ def hips_frame(self) -> str:
         """HiPS coordinate frame (`str`)."""
         return self.data['hips_frame']
 
+    @property
+    def astropy_frame(self) -> str:
+        """Astropy coordinate frame (`str`)."""
+        return self.hips_to_astropy_frame_mapping[self.hips_frame]
+
     @property
     def hips_order(self) -> int:
         """HiPS order (`int`)."""
@@ -108,3 +119,13 @@ def hips_order(self) -> int:
     def tile_format(self) -> str:
         """HiPS tile format (`str`)."""
         return self.data['hips_tile_format']
+
+    @property
+    def access_url(self):
+        """HiPS access url"""
+        return self.data['moc_access_url'].rsplit('/', 1)[0]
+
+    @property
+    def tile_access_url(self):
+        """Tile access URL for a HiPS surveys"""
+        return self.access_url + '/Norder' + str(self.hips_order) + '/Dir0/'

hips/tiles/tests/test_description.py

@@ -21,8 +21,17 @@ def test_hips_version(self):
     def test_hips_frame(self):
         assert self.hips_survey_property.hips_frame == 'equatorial'
 
+    def test_astropy_frame(self):
+        assert self.hips_survey_property.astropy_frame == 'icrs'
+
     def test_hips_order(self):
         assert self.hips_survey_property.hips_order == 9
 
     def test_tile_format(self):
         assert self.hips_survey_property.tile_format == 'jpeg'
+
+    def test_access_url(self):
+        assert self.hips_survey_property.access_url == 'http://alasky.u-strasbg.fr/DSS/DSSColor'
+
+    def test_tile_access_url(self):
+        assert self.hips_survey_property.tile_access_url == 'http://alasky.u-strasbg.fr/DSS/DSSColor/Norder9/Dir0/'

hips/tiles/tests/test_tile.py

@@ -1,8 +1,10 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
-from numpy.testing import assert_equal
+import numpy as np
 from astropy.tests.helper import remote_data
-from ..tile import HipsTile
-from ..tile_meta import HipsTileMeta
+from numpy.testing import assert_allclose
+from numpy.testing import assert_equal
+
+from ..tile import HipsTile, HipsTileMeta
 
 
 class TestHipsTile:
@@ -41,3 +43,35 @@ def test_fetch_read_write_fits(self, tmpdir):
         # print(tile2.data.sum())
         # print((tile == tile2).all())
         # assert tile == tile2
+
+
+class TestHipsTileMeta:
+    @classmethod
+    def setup_class(cls):
+        cls.meta = HipsTileMeta(order=3, ipix=450, file_format='fits', frame='icrs', tile_width=512)
+
+    def test_path(self):
+        assert str(self.meta.path) == 'hips/tiles/tests/data'
+
+    def test_filename(self):
+        assert self.meta.filename == 'Npix450.fits'
+
+    def test_full_path(self):
+        assert str(self.meta.full_path) == 'hips/tiles/tests/data/Npix450.fits'
+
+    def test_nside(self):
+        assert self.meta.nside == 8
+
+    def test_dst(self):
+        dst = np.array([[511, 0], [511, 511], [0, 511], [0, 0]])
+        assert_allclose(self.meta.dst, dst)
+
+    def test_skycoord_corners(self):
+        assert_allclose(self.meta.skycoord_corners.data.lat.deg, [-24.624318, -30., -35.685335, -30.])
+        assert_allclose(self.meta.skycoord_corners.data.lon.deg, [264.375, 258.75, 264.375, 270.])
+        assert self.meta.skycoord_corners.frame.name == 'icrs'
+
+        meta = HipsTileMeta(order=3, ipix=450, file_format='fits', frame='galactic', tile_width=512)
+        assert_allclose(meta.skycoord_corners.data.lat.deg, [-24.624318, -30., -35.685335, -30.])
+        assert_allclose(meta.skycoord_corners.data.lon.deg, [264.375, 258.75, 264.375, 270.])
+        assert meta.skycoord_corners.frame.name == 'galactic'