Add test case for get_order function

.travis.yml

@@ -46,7 +46,7 @@ matrix:
     include:
         # Main build -- used for coverage
         - os: linux
-          env: SETUP_CMD='test --coverage'
+          env: SETUP_CMD='test --coverage --remote-data -V'
 
         # Docs build
         - os: linux
@@ -55,7 +55,7 @@ matrix:
 
         # Check that install / tests work with minimal required dependencies
         - os: linux
-          env: SETUP_CMD='test'
+          env: SETUP_CMD='test -V'
                CONDA_DEPENDENCIES='healpy scikit-image Pillow'
 
         # Try MacOS X

hips/draw/simple.py

@@ -10,9 +10,13 @@
 __all__ = [
     'draw_sky_image',
     'make_sky_image',
-    'SimpleTilePainter'
+    'SimpleTilePainter',
+    'compute_matching_hips_order'
 ]
 
+__doctest_skip__ = [
+    'compute_matching_hips_order',
+]
 
 # TODO: Fix type annotation issue
 def draw_sky_image(geometry: WCSGeometry, tiles: Generator[HipsTile, Any, Any], hips_survey: HipsSurveyProperties) -> np.ndarray:
@@ -117,20 +121,57 @@ def fetch_tiles(healpix_pixel_indices: np.ndarray, order: int, hips_survey: Hips
         yield tile
 
 
-def get_order(geometry: WCSGeometry, hips_survey: HipsSurveyProperties) -> int:
-    """Compute the tile order suited for the given geometry and hips_survey"""
+def compute_matching_hips_order(geometry: WCSGeometry, hips_survey: HipsSurveyProperties) -> int:
+    """Compute HiPS tile order matching a given image pixel size.
+
+    Parameters
+    ----------
+    geometry : WCSGeometry
+        Geometry of the output image
+    hips_survey : HipsSurveyProperties
+        An object of HipsSurveyProperties
 
+    Returns
+    -------
+    'int'
+        Returns HiPS order
+
+    Examples
+    --------
+    >>> from hips.draw import compute_matching_hips_order
+    >>> from astropy.coordinates import SkyCoord
+    >>> url = 'http://alasky.unistra.fr/DSS/DSS2Merged/properties'
+    >>> hips_survey = HipsSurveyProperties.fetch(url)
+    >>> geometry = WCSGeometry.create_simple(
+    ...     skydir=SkyCoord(0, 0, unit='deg', frame='icrs'),
+    ...     width=2000, height=1000, fov="3 deg",
+    ...     coordsys='icrs', projection='AIT'
+    ... )
+    >>> compute_matching_hips_order(geometry, hips_survey)
+    7
+    """
+
+    # Sky image angular resolution (pixel size in degree)
     resolution = np.min(proj_plane_pixel_scales(geometry.wcs))
+    desired_order = _get_hips_order_for_resolution(hips_survey.tile_width, resolution)
+    # Return the desired order, or the highest resolution available.
+    # Note that HiPS never has resolution less than 3,
+    # and that limit is handled in _get_hips_order_for_resolution
+    return np.min([desired_order, hips_survey.hips_order])
+
 
-    tile_order = np.log2(hips_survey.tile_width)
+def _get_hips_order_for_resolution(tile_width, resolution):
+    """Finding the best HiPS order by looping through all possible options."""
+    tile_order = np.log2(tile_width)
     full_sphere_area = 4 * np.pi * np.square(180 / np.pi)
+    # 29 is the maximum order supported by healpy and 3 is the minimum order
     for candidate_tile_order in range(3, 29 + 1):
-        tile_resolution = np.sqrt(full_sphere_area / 12 / 4**(candidate_tile_order + tile_order))
-
+        tile_resolution = np.sqrt(full_sphere_area / 12 / 4 ** (candidate_tile_order + tile_order))
+        # Finding the smaller tile order with a resolution equal or better than geometric resolution
         if tile_resolution <= resolution:
             break
 
-    return np.min([candidate_tile_order, hips_survey.hips_order])
+    return candidate_tile_order
 
 
 def make_sky_image(geometry: WCSGeometry, hips_survey: HipsSurveyProperties) -> np.ndarray:
@@ -150,7 +191,7 @@ def make_sky_image(geometry: WCSGeometry, hips_survey: HipsSurveyProperties) ->
     data : `~numpy.ndarray`
         Output image pixels
     """
-    order = get_order(geometry, hips_survey)
+    order = compute_matching_hips_order(geometry, hips_survey)
     healpix_pixel_indices = compute_healpix_pixel_indices(
         wcs_geometry=geometry,
         order=order,

hips/draw/tests/test_simple.py

@@ -1,13 +1,16 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
+import pytest
 import numpy as np
+import healpy as hp
+from astropy.coordinates import SkyCoord
 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 ..simple import make_sky_image, draw_sky_image, compute_matching_hips_order, _get_hips_order_for_resolution
 from ...tiles import HipsSurveyProperties, HipsTileMeta, HipsTile
+from ...utils import WCSGeometry
 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)
@@ -44,14 +47,46 @@ def test_draw_sky_image():
 
 @remote_data
 def test_make_sky_image():
-    url = 'https://raw.githubusercontent.com/hipspy/hips-extra/master/datasets/samples/DSS2Red/properties'
+    # The same example is used in the high level docs getting started page
+    url = 'http://alasky.unistra.fr/DSS/DSS2Merged/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)
+    assert_allclose(data[200, 994], 2213.30874796)
+    assert_allclose(data[200, 995], 2296.93885940)
+
+hips_order_pars = [
+    dict(order=7, fov="3 deg"),
+    dict(order=5, fov="10 deg"),
+    dict(order=4, fov="15 deg"),
+]
+
+
+@requires_hips_extra()
+@pytest.mark.parametrize('pars', hips_order_pars)
+def test_compute_matching_hips_order(pars):
+    full_path = get_hips_extra_file('datasets/samples/2MASS6XH/properties')
+    hips_survey = HipsSurveyProperties.read(filename=full_path)
+    geometry = WCSGeometry.create_simple(
+        skydir=SkyCoord(0, 0, unit='deg', frame='icrs'),
+        width=2000, height=1000, fov=pars['fov'],
+        coordsys='icrs', projection='AIT'
+    )
+    assert compute_matching_hips_order(geometry, hips_survey) == pars['order']
+
+get_hips_order_for_resolution_pars = [
+    dict(tile_width=512, resolution=0.01232, resolution_res=0.06395791924665553, order=4),
+    dict(tile_width=256, resolution=0.0016022, resolution_res=0.003997369952915971, order=8),
+    dict(tile_width=128, resolution=0.00009032, resolution_res=0.00012491781102862408, order=13),
+]
 
 
+@pytest.mark.parametrize('pars', get_hips_order_for_resolution_pars)
+def test_get_hips_order_for_resolution(pars):
+    hips_order = _get_hips_order_for_resolution(pars['tile_width'], pars['resolution'])
+    assert hips_order == pars['order']
+    hips_resolution = hp.nside2resol(hp.order2nside(hips_order))
+    assert_allclose(hips_resolution, pars['resolution_res'])
 # TODO: add tests for SimpleTilePainter with asserts on the intermediate computed things.