Update drawing code

hips/draw/simple.py

@@ -34,48 +34,52 @@ def draw_sky_image(geometry: WCSGeometry, tiles: List[HipsTile]) -> np.ndarray:
     np.ndarray
         Returns a numpy array containing all HiPS tiles projected onto it
     """
-    # TODO: copy over (and clean up a little) the drawing code you already have in notebooks.
-    # Suggestion: don't debug the code now, just put it here.
-    # First implement simple tile caching so that you can draw without needing to make web queries.
-
-    # Probably this function should just loop over tiles and call a `_draw_sky_image_one_tile`
-    # helper function, and sum up the results.
     all_sky = np.zeros(geometry.shape)
     for tile in tiles:
-        nside = hp.order2nside(tile.meta.order)
-        theta, phi = boundaries(nside, tile.meta.ipix)
-        radec = SkyCoord(ra=phi, dec=np.pi / 2 - theta, unit='radian', frame='icrs')
-        corners = []
-        for i in range(len(radec.ra.deg)):
-            corners.append([radec.ra.deg[i], radec.dec.deg[i]])
-        all_sky += _draw_sky_image_one_tile(corners, WCS(tile.header), tile, all_sky.shape)
+        draw_tile = DrawSkyImageOneTile(geometry, tile, all_sky.shape)
+        draw_tile.compute_corners()
+        draw_tile.apply_projection()
+        all_sky += draw_tile.warp_image()
     return all_sky
 
 
-def _draw_sky_image_one_tile(corners: list, wcs: WCS, tile: HipsTile, shape: tuple) -> np.ndarray:
+class DrawSkyImageOneTile:
     """A private function for drawing a tile over a sky image.
 
     Parameters
     ----------
-    corners : `~astropy.coordinates.SkyCoord`
-        Four corners of a HiPS tile
-    wcs : `~astropy.wcs.WCS`
-        WCS projection of a HiPS tile
-    tile : 'HipsTile'
-        An object of HipsTile
+    geometry : 'WCSGeometry'
+        An object of WCSGeometry
+    tile : HipsTile
+       An object of HipsTile
     shape ; tuple
         Shape of the all-sky image
-
-    Returns
-    -------
-    np.ndarray
-        Returns a numpy array containing the projected HiPS tile
     """
-    src = wcs.wcs_world2pix(corners, 0)
-    dst = np.array([[511, 0], [511, 511], [0, 511], [0, 0]])
-    pt = tf.ProjectiveTransform()
-    pt.estimate(src, dst)
-    return tf.warp(tile.data, pt, output_shape=shape)
+
+    def __init__(self, geometry: WCSGeometry, tile: HipsTile, shape: tuple) -> None:
+        self.geometry = geometry
+        self.tile = tile
+        self.shape = shape
+        self.wcs = WCS(tile.header)
+        self.corners = None
+        self.pt = None
+
+    def compute_corners(self) -> None:
+        nside = hp.order2nside(self.tile.meta.order)
+        theta, phi = boundaries(nside, self.tile.meta.ipix)
+        radec = SkyCoord(ra=phi, dec=np.pi / 2 - theta, unit='radian', frame='icrs')
+        self.corners = []
+        for i in range(len(radec.ra.deg)):
+            self.corners.append([radec.ra.deg[i], radec.dec.deg[i]])
+
+    def apply_projection(self) -> None:
+        src = self.geometry.wcs.wcs_world2pix(self.corners, 0)
+        dst = np.array([[511, 0], [511, 511], [0, 511], [0, 0]])
+        self.pt = tf.ProjectiveTransform()
+        self.pt.estimate(src, dst)
+
+    def warp_image(self) -> np.ndarray:
+        return tf.warp(self.tile.data, self.pt, output_shape=self.shape)
 
 
 def _fetch_tiles(healpix_pixel_indices: np.ndarray, order: int, hips_survey: HipsSurveyProperties) -> 'HipsTile':
@@ -118,26 +122,21 @@ def make_sky_image(geometry: WCSGeometry, hips_survey: HipsSurveyProperties) ->
 
     Examples
     --------
-    Define which image you want:
-
-    >>> geometry = tbd
-
-    Define which HiPS survey you want:
-
-    >>> hips_survey = tbd
-
-    Compute the image:
-
-    >>> from hips import make_sky_image
-    >>> data = make_sky_image(geometry, hips_survey)
-
-    If you want, you can save the image to file:
     >>> from astropy.io import fits
+    >>> 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.header)
     >>> hdu.writeto('my_image.fits')
     """
-    # Compute list of tiles needed
-
     healpix_pixel_indices = compute_healpix_pixel_indices(geometry, hips_survey.hips_order)
     path = Path(os.environ['HIPS_EXTRA'])
     tiles_path = path / 'datasets' / 'samples' / 'DSS2Red' / 'Norder3' / 'Dir0'
@@ -152,7 +151,6 @@ def make_sky_image(geometry: WCSGeometry, hips_survey: HipsSurveyProperties) ->
     # TODO: this isn't a good API. Will become better when we have a cache.
     # tiles = _fetch_tiles(healpix_pixel_indices, order, hips_survey)
 
-    # Draw the image
     image_data = draw_sky_image(geometry, tiles)
 
     return image_data

hips/draw/tests/test_simple.py

@@ -34,7 +34,7 @@ def test_draw_sky_image():
     assert data.shape == geometry.shape
     assert data.dtype == np.float64
 
-    assert data[0, 0] == 0.060903854985925092
+    assert data[0, 0] == 0.0
 
 
 @remote_data
@@ -47,5 +47,4 @@ def test_make_sky_image():
 
     assert data.shape == geometry.shape
     assert data.dtype == np.float64
-
-    assert data[0, 0] == 0.060903854985925092
+    assert data[800, 1000] == 1794.7673494847763