diff --git a/hips/draw/simple.py b/hips/draw/simple.py index f82df43..9ab55e9 100644 --- a/hips/draw/simple.py +++ b/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 diff --git a/hips/draw/tests/test_simple.py b/hips/draw/tests/test_simple.py index 8f293a7..8c7aa97 100644 --- a/hips/draw/tests/test_simple.py +++ b/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