From fc5a870bf4d8215a8818d3068749ea97e7029932 Mon Sep 17 00:00:00 2001 From: Christoph Deil Date: Thu, 20 Jul 2017 15:35:11 +0200 Subject: [PATCH] Add all-sky image class and major update to existing tile code --- hips/draw/simple.py | 33 ++----- hips/draw/tests/test_simple.py | 20 ++-- hips/tiles/__init__.py | 1 + hips/tiles/allsky.py | 163 +++++++++++++++++++++++++++++++ hips/tiles/tests/test_allsky.py | 127 ++++++++++++++++++++++++ hips/tiles/tests/test_tile.py | 19 +++- hips/tiles/tile.py | 93 ++++++++++++++---- hips/utils/healpix.py | 8 ++ hips/utils/tests/test_healpix.py | 11 ++- 9 files changed, 419 insertions(+), 56 deletions(-) create mode 100644 hips/tiles/allsky.py create mode 100644 hips/tiles/tests/test_allsky.py diff --git a/hips/draw/simple.py b/hips/draw/simple.py index 727f6c0..e9237e0 100644 --- a/hips/draw/simple.py +++ b/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 diff --git a/hips/draw/tests/test_simple.py b/hips/draw/tests/test_simple.py index 1850a69..d3368f7 100644 --- a/hips/draw/tests/test_simple.py +++ b/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) diff --git a/hips/tiles/__init__.py b/hips/tiles/__init__.py index 7a93265..a96dcf3 100644 --- a/hips/tiles/__init__.py +++ b/hips/tiles/__init__.py @@ -2,3 +2,4 @@ """Classes and functions to manage HiPS tiles.""" from .tile import * from .surveys import * +from .allsky import * diff --git a/hips/tiles/allsky.py b/hips/tiles/allsky.py new file mode 100644 index 0000000..ef385db --- /dev/null +++ b/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), + ] diff --git a/hips/tiles/tests/test_allsky.py b/hips/tiles/tests/test_allsky.py new file mode 100644 index 0000000..80bfe2e --- /dev/null +++ b/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) diff --git a/hips/tiles/tests/test_tile.py b/hips/tiles/tests/test_tile.py index 2957e40..d5b087b 100644 --- a/hips/tiles/tests/test_tile.py +++ b/hips/tiles/tests/test_tile.py @@ -10,7 +10,9 @@ class TestHipsTileMeta: @classmethod def setup_class(cls): - cls.meta = HipsTileMeta(order=3, ipix=450, file_format='fits', frame='icrs') + # Note: we're intentionally using positional args + # to make sure tests fail if the order of args is changed + cls.meta = HipsTileMeta(3, 450, 'fits', 'icrs', 512) def test_order(self): assert self.meta.order == 3 @@ -24,10 +26,23 @@ def test_file_format(self): def test_frame(self): assert self.meta.frame == 'icrs' + def test_width(self): + assert self.meta.width == 512 + def test_repr(self): - expected = "HipsTileMeta(order=3, ipix=450, file_format='fits', frame='icrs')" + expected = "HipsTileMeta(order=3, ipix=450, file_format='fits', frame='icrs', width=512)" assert repr(self.meta) == expected + def test_eq(self): + assert self.meta == HipsTileMeta(3, 450, 'fits', 'icrs', 512) + assert self.meta != HipsTileMeta(4, 450, 'fits', 'icrs', 512) + + def test_copy(self): + meta = self.meta.copy() + meta.ipix = 42 + assert self.meta.ipix == 450 + assert meta.ipix == 42 + def test_skycoord_corners(self): coord = self.meta.skycoord_corners assert_allclose(coord.data.lat.deg, [-24.624318, -30., -35.685335, -30.]) diff --git a/hips/tiles/tile.py b/hips/tiles/tile.py index 5a110c7..f19e120 100644 --- a/hips/tiles/tile.py +++ b/hips/tiles/tile.py @@ -1,4 +1,6 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst +from typing import List, Tuple +from copy import deepcopy import warnings import urllib.request from io import BytesIO @@ -11,7 +13,6 @@ from astropy.io import fits from ..utils import healpix_pixel_corners from .io import tile_default_url, tile_default_path -from typing import List __all__ = [ 'HipsTileMeta', @@ -23,6 +24,38 @@ ] +# TODO: this could be a dict. Would that be better? +def compute_image_shape(width: int, height: int, fmt: str) -> Tuple[int]: + """Compute numpy array shape for a given 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 + + Parameters + ---------- + width, height : int + Width and height of the image + fmt : {'fits', 'jpg', 'png'} + Image format + + Returns + ------- + shape : tuple + Numpy array shape + """ + if fmt == 'fits': + return height, width + elif fmt == 'jpg': + return height, width, 3 + elif fmt == 'png': + return height, width, 4 + else: + return ValueError(f'Invalid format: {fmt}') + + class HipsTileMeta: """HiPS tile metadata. @@ -36,13 +69,15 @@ class HipsTileMeta: File format frame : {'icrs', 'galactic', 'ecliptic'} Sky coordinate frame + width : int + Tile width (tiles always have width = height) Examples -------- >>> from hips import HipsTileMeta - >>> tile_meta = HipsTileMeta(order=3, ipix=450, file_format='fits', frame='icrs') + >>> tile_meta = HipsTileMeta(order=3, ipix=450, file_format='fits') >>> tile_meta - HipsTileMeta(order=3, ipix=450, file_format='fits', frame='icrs') + HipsTileMeta(order=3, ipix=450, file_format='fits', frame='icrs', width=512) >>> tile_meta.skycoord_corners None: + def __init__(self, order: int, ipix: int, file_format: str, + frame: str = 'icrs', width: int = 512) -> None: self.order = order self.ipix = ipix self.file_format = file_format self.frame = frame + self.width = width def __repr__(self): return ( - f'HipsTileMeta(order={self.order}, ipix={self.ipix}, ' - f'file_format={self.file_format!r}, frame={self.frame!r})' + 'HipsTileMeta(' + f'order={self.order}, ipix={self.ipix}, ' + f'file_format={self.file_format!r}, frame={self.frame!r}, ' + f'width={self.width}' + ')' ) def __eq__(self, other: 'HipsTileMeta') -> bool: return ( self.order == other.order and self.ipix == other.ipix and - self.file_format == other.file_format + self.file_format == other.file_format and + self.frame == other.frame and + self.width == other.width ) + def copy(self): + """An independent copy.""" + return deepcopy(self) + @property def skycoord_corners(self) -> SkyCoord: """Tile corner sky coordinates (`~astropy.coordinates.SkyCoord`).""" @@ -123,6 +169,7 @@ class HipsTile: def __init__(self, meta: HipsTileMeta, raw_data: BytesIO) -> None: self.meta = meta self.raw_data = raw_data + self._data = None def __eq__(self, other: 'HipsTile') -> bool: return ( @@ -170,31 +217,44 @@ def from_numpy(cls, meta: HipsTileMeta, data: np.ndarray) -> 'HipsTile': @property def children(self) -> List['HipsTile']: """Create four children tiles from parent tile.""" - children_tiles = [] w = self.data.shape[0] // 2 - child_data = [ + data = [ self.data[0: w, 0: w], self.data[0: w, w: w * 2], self.data[w: w * 2, 0: w], self.data[w: w * 2, w: w * 2] ] - for index, data in enumerate(child_data): + tiles = [] + for idx in range(4): meta = HipsTileMeta( self.meta.order + 1, - self.meta.ipix * 4 + index, + self.meta.ipix * 4 + idx, self.meta.file_format, self.meta.frame ) - children_tiles.append(self.from_numpy(meta, data)) + tile = self.from_numpy(meta, data[idx]) + tiles.append(tile) - return children_tiles + return tiles @property def data(self) -> np.ndarray: - """Tile pixel data (`~numpy.ndarray`).""" - fmt = self.meta.file_format - bio = BytesIO(self.raw_data) + """Tile pixel data (`~numpy.ndarray`). + + This is a cached property, it will only be computed once. + """ + if self._data is None: + self._data = self._to_numpy( + self.raw_data, + self.meta.file_format, + ) + + return self._data + + @staticmethod + def _to_numpy(raw_data, fmt): + bio = BytesIO(raw_data) if fmt == 'fits': # At the moment CDS is serving FITS tiles in non-standard FITS files @@ -206,7 +266,6 @@ def data(self) -> np.ndarray: warnings.simplefilter('ignore', VerifyWarning) with fits.open(bio) as hdu_list: data = hdu_list[0].data - # header = hdu_list[0].header elif fmt in {'jpg', 'png'}: with Image.open(bio) as image: data = np.array(image) diff --git a/hips/utils/healpix.py b/hips/utils/healpix.py index 66cd4ff..f29dd6d 100644 --- a/hips/utils/healpix.py +++ b/hips/utils/healpix.py @@ -18,10 +18,13 @@ from .wcs import WCSGeometry __all__ = [ + 'healpix_order_to_npix', 'healpix_skycoord_to_theta_phi', 'healpix_theta_phi_to_skycoord', + 'healpix_pixel_corners', 'healpix_pixels_in_sky_image', + 'hips_order_for_pixel_resolution', ] @@ -34,6 +37,11 @@ """HiPS always uses the nested HEALPix pixel numbering scheme.""" +def healpix_order_to_npix(order: int) -> int: + """HEALPix order to npix.""" + return hp.nside2npix(hp.order2nside(order)) + + def healpix_skycoord_to_theta_phi(skycoord: SkyCoord) -> Tuple[float, float]: """Convert SkyCoord to theta / phi as used in healpy.""" theta = np.pi / 2 - skycoord.data.lat.radian diff --git a/hips/utils/tests/test_healpix.py b/hips/utils/tests/test_healpix.py index 7824f9a..abb2f0c 100644 --- a/hips/utils/tests/test_healpix.py +++ b/hips/utils/tests/test_healpix.py @@ -2,8 +2,17 @@ import pytest from numpy.testing import assert_allclose import healpy as hp -from ..healpix import healpix_pixel_corners, healpix_pixels_in_sky_image, hips_order_for_pixel_resolution from ..testing import make_test_wcs_geometry +from ..healpix import ( + healpix_order_to_npix, + healpix_pixel_corners, + healpix_pixels_in_sky_image, + hips_order_for_pixel_resolution, +) + + +def test_order_to_npix(): + assert healpix_order_to_npix(3) == 768 def test_healpix_pixel_corners():