diff --git a/docs/hipsgen.py b/docs/hipsgen.py new file mode 100644 index 0000000..d3db83f --- /dev/null +++ b/docs/hipsgen.py @@ -0,0 +1,52 @@ +"""Example how to generate HiPS from HEALPix. + +We will use +""" +import logging +import numpy as np +import hips +from astropy.coordinates import SkyCoord +from astropy_healpix import HEALPix + + +def make_healpix_data(hips_order, tile_width, file_format): + """Silly example of HEALPix data. + + Angular distance in deg to (0, 0) as pixel values. + """ + nside = tile_width * (2 ** hips_order) + healpix = HEALPix(nside=nside) + ipix = np.arange(healpix.npix) + lon, lat = healpix.healpix_to_lonlat(ipix) + coord = SkyCoord(lon, lat) + center = SkyCoord(0, 0, unit="deg") + data = coord.separation(center).deg + + if file_format == "fits": + return data + elif file_format == "jpg": + data = data.astype("uint8") + return np.moveaxis([data, data + 1, data + 2], 0, -1) + elif file_format == "png": + data = data.astype("uint8") + return np.moveaxis([data, data + 1, data + 2, data + 3], 0, -1) + else: + raise ValueError() + + +if __name__ == "__main__": + logging.basicConfig(level=logging.INFO) + + file_format = "png" + hips_order = 3 + tile_width = 4 + + hpx_data = make_healpix_data(hips_order, tile_width, file_format) + + hips.healpix_to_hips( + hpx_data=hpx_data, + hips_order=hips_order, + tile_width=tile_width, + base_path="test123", + file_format=file_format, + ) diff --git a/docs/hipsgen.rst b/docs/hipsgen.rst new file mode 100644 index 0000000..663c1f0 --- /dev/null +++ b/docs/hipsgen.rst @@ -0,0 +1,72 @@ +.. include:: references.txt + +.. _hipsgen: + +************* +Generate HiPS +************* + +The `~hips.healpix_to_hips` function can be used to generate HiPS data from all-sky HEALPix images. +We give an example below. + +Note that this functionality is very limited, and for most applications you will want to use the +Java ``hipsgen`` tool or the graphical user interface to ``hipsgen`` from the Aladin desktop application. + +To use `~hips.healpix_to_hips` you have to pass in the all-sky HEALPix data as a Numpy array +(so it has to completely fit into memory). For grayscale images, the HEALPix data array is +one-dimensional, RGB images with array shape ``(npix, 3)`` can be converted to JPEG HiPS tiles, +and RGBA (where A is the transparency channel) images with array shape ``(npix, 4)`` can be +converted to PNG HiPS tiles. + +Some data is distributed directly in HEALPix format, e.g. from all-sky surveys like Planck. +Some data from high-energy telescopes like Fermi-LAT consists of event lists can be used +to compute HEALPix images. And WCS-based images can also be reprojected to HEALPix images +using e.g. `reproject `__. +For these use cases `~hips.healpix_to_hips` can be used for the final all-sky HEALPix +map to HiPS conversion step, giving you full flexibility over the processing and pixel +color scales etc from Python in the pre-processing steps. + +Limitations +----------- + +Currently the following limitations exist: + +- No mosaic functionality from WCS images + (a lot of work to implement) +- Only one HiPS resolution is generated at the moment + (not hard to write a script to make multiple resolutions) +- No ``Allsky.fits`` generated yet + (shouldn't be hard to generate) + +Dataset +======= + +bla bla + +Example +======= + +Generate HiPS: + +.. code-block:: bash + + $ python docs/hipsgen.py + +Open HiPS in Aladin Lite: + +.. code-block:: bash + + $ cd test123 + $ python -m http.server + # Then open http://localhost:8000/ in your webbrowser + +Or open HiPS in Aladin Desktop: + +.. code-block:: bash + + $ TODO: how to load the folder in Aladin + +Here's the ``hipsgen.py`` script: + +.. literalinclude:: hipsgen.py + diff --git a/docs/index.rst b/docs/index.rst index f152c9b..01a7d52 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -15,6 +15,7 @@ A Python astronomy package for HiPS : Hierarchical Progressive Surveys. about installation getting_started + hipsgen api changelog diff --git a/hips/draw/__init__.py b/hips/draw/__init__.py index 14575cd..fcb29b8 100644 --- a/hips/draw/__init__.py +++ b/hips/draw/__init__.py @@ -3,3 +3,4 @@ from .paint import * from .ui import * from .healpix import * +from .hipsgen import * diff --git a/hips/draw/healpix.py b/hips/draw/healpix.py index 00619ee..f4795c1 100644 --- a/hips/draw/healpix.py +++ b/hips/draw/healpix.py @@ -1,18 +1,22 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst +from typing import List import logging -from pathlib import Path import numpy as np -from astropy_healpix import healpy as hp -from ..tiles import HipsTile, HipsTileMeta, HipsSurveyProperties +from ..tiles import HipsTile, HipsTileMeta from ..utils.healpix import hips_tile_healpix_ipix_array -__all__ = ["healpix_to_hips_tile", "healpix_to_hips"] +__all__ = ["healpix_to_hips_tile", "healpix_to_hips_tiles"] log = logging.getLogger(__name__) def healpix_to_hips_tile( - hpx_data: np.ndarray, tile_width: int, tile_idx: int, file_format: str, frame: str + hpx_data: np.ndarray, + hips_order: int, + tile_width: int, + tile_idx: int, + file_format: str, + frame: str = "icrs", ) -> HipsTile: """Create single HiPS tile from HEALPix data. @@ -20,6 +24,8 @@ def healpix_to_hips_tile( ---------- hpx_data : `~numpy.ndarray` Healpix data stored in the "nested" scheme. + hips_order : int + HEALPix order for the HiPS tiles tile_width : int Width of the hips tile. tile_idx : int @@ -39,17 +45,38 @@ def healpix_to_hips_tile( offset_ipix = tile_idx * tile_width ** 2 ipix = hpx_ipix + offset_ipix - data = hpx_data[ipix] + + hpx_data = np.asarray(hpx_data) + if file_format == "fits": + if hpx_data.ndim != 1: + raise ValueError( + f"Invalid hpx_data.ndim = {hpx_data.ndim}." + " Must be ndim = 1 for file_format='fits'." + ) + data = hpx_data[ipix] + elif file_format == "jpg": + if hpx_data.ndim != 2 or hpx_data.shape[1] != 3: + raise ValueError( + f"Invalid hpx_data.shape = {hpx_data.shape}." + " Must be shape = (npix, 3) to represent RGB color images." + ) + data = hpx_data[ipix, :] + elif file_format == "png": + if hpx_data.ndim != 2 or hpx_data.shape[1] != 4: + raise ValueError( + f"Invalid hpx_data.shape = {hpx_data.shape}." + " Must be shape = (npix, 4) to represent RGBA color images." + ) + data = hpx_data[ipix, :] + else: + raise ValueError(f"Invalid file_format: {file_format!r}") # np.rot90 returns a rotated view so we make a copy here # because the view information is lost on fits io data = np.rot90(data).copy() - hpx_nside = hp.npix2nside(hpx_data.size / tile_width ** 2) - hpx_order = int(np.log2(hpx_nside)) - meta = HipsTileMeta( - order=hpx_order, + order=hips_order, ipix=tile_idx, file_format=file_format, frame=frame, @@ -59,50 +86,41 @@ def healpix_to_hips_tile( return HipsTile.from_numpy(meta=meta, data=data) -def healpix_to_hips(hpx_data, tile_width, base_path, file_format, frame): - """Convert HEALPix image to HiPS. - - This function writes the HiPS to disk. - If you don't want that, use `healpix_to_hips_tile` directly. +def healpix_to_hips_tiles( + hpx_data: np.ndarray, + hips_order: int, + tile_width: int, + file_format: str, + frame: str = "icrs", +) -> List[HipsTile]: + """Convert HEALPix image to HiPS tiles, Parameters ---------- hpx_data : `~numpy.ndarray` - Healpix data stored in the "nested" scheme. + HEALPix data stored in the "nested" scheme. + hips_order : int + HEALPix order for the HiPS tiles tile_width : int Width of the hips tiles. - base_path : str or `~pathlib.Path` - Base path. file_format : {'fits', 'jpg', 'png'} HiPS tile file format frame : {'icrs', 'galactic', 'ecliptic'} Sky coordinate frame - """ - base_path = Path(base_path) - base_path.mkdir(exist_ok=True, parents=True) - path = base_path / "properties" - log.info(f"Writing {path}") - HipsSurveyProperties( - { - "hips_tile_format": file_format, - "hips_tile_width": tile_width, - "hips_frame": frame, - } - ).write(path) - - n_tiles = hpx_data.size // tile_width ** 2 + Returns + ------- + tiles : Generator of HipsTile + HiPS tiles + """ + n_tiles = 12 * (2 ** hips_order) for tile_idx in range(n_tiles): - tile = healpix_to_hips_tile( + yield healpix_to_hips_tile( hpx_data=hpx_data, + hips_order=hips_order, tile_width=tile_width, tile_idx=tile_idx, file_format=file_format, frame=frame, ) - - path = base_path / tile.meta.tile_default_path - log.info(f"Writing {path}") - path.parent.mkdir(exist_ok=True, parents=True) - tile.write(path) diff --git a/hips/draw/hipsgen.py b/hips/draw/hipsgen.py new file mode 100644 index 0000000..9413fc2 --- /dev/null +++ b/hips/draw/hipsgen.py @@ -0,0 +1,154 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +import logging +import numpy as np +from ..tiles import HipsSurveyProperties, HipsTileAllskyArray, HipsTile, HipsTileMeta +from .healpix import healpix_to_hips_tiles +from pathlib import Path + +__all__ = ["healpix_to_hips", "HipsSurvey"] + +log = logging.getLogger(__name__) + +HTML_TEMPLATE = r""" + + + + +
+ + +""" + + +# TODO: change to class with each step as a method? +def healpix_to_hips( + hpx_data: np.ndarray, + hips_order: int, + tile_width: int, + base_path: str, + file_format: str, + frame: str = "icrs", +): + """Convert HEALPix image to HiPS. + + Directly writes files to output folder. + + Parameters + ---------- + hpx_data : `~numpy.ndarray` + Healpix data stored in the "nested" scheme. + hips_order : int + HEALPix order for the HiPS tiles + tile_width : int + Width of the hips tiles. + base_path : str or `~pathlib.Path` + Base path. + file_format : {'fits', 'jpg', 'png'} + HiPS tile file format + frame : {'icrs', 'galactic', 'ecliptic'} + Sky coordinate frame + """ + # Make and write properties file + base_path = Path(base_path) + base_path.mkdir(exist_ok=True, parents=True) + path = base_path / "properties" + log.info(f"Writing {path}") + properties = HipsSurveyProperties( + { + "hips_tile_format": file_format, + "hips_tile_width": tile_width, + "hips_order": hips_order, + "hips_frame": frame, + } + ) + properties.write(path) + + # Make and write index.html + txt = HTML_TEMPLATE.format_map({"name": "test123", "imgFormat": file_format}) + path = base_path / "index.html" + log.info(f"Writing {path}") + path.write_text(txt) + + # Make and write tiles + for tile in healpix_to_hips_tiles( + hpx_data, hips_order, tile_width, file_format, frame + ): + path = base_path / tile.meta.tile_default_path + log.info(f"Writing {path}") + path.parent.mkdir(exist_ok=True, parents=True) + tile.write(path) + + # Make and write allsky file + hips_survey = HipsSurvey(base_path) + allsky = hips_survey.make_allsky(order=hips_order, file_format=file_format) + path = base_path / f"Norder{hips_order}/Allsky.{file_format}" + log.info(f"Writing {path}") + allsky.write(path) + + +# TODO: generalise and move to a better location +class HipsSurvey: + """HiPS survey container. + + Represents one HiPS survey, acting as a manager + to interact with the various pieces that make up a HiPS: + + - One `HipsSurveyProperties` + - Many `HipsTile` + - One `HipsTileAllskyArray` per order + + TODO: do we need several classes, for the different cases, + or can they be handled by one `HipsSurvey` class? + + - in memory + - on local disk + - on server + + Functionality should include: + + - locate, read, write files + - scan and print summary of contents (e.g. which orders / tiles are present) + - copy / clone HiPS (possibly only small parts, with selections) locally and from servers + - generate allsky from tiles? put this functionality here or somewhere else? + + Parameters + ---------- + base_path : `pathlib.Path` + Base path (folder where the ``properties`` file is located) + """ + + def __init__(self, base_path): + self.base_path = base_path + + def make_allsky(self, order: int, file_format: str) -> HipsTileAllskyArray: + tiles = self.read_tiles(order, file_format) + # TODO: make from_tiles work with generators. For now we have to make a list + tiles = list(tiles) + return HipsTileAllskyArray.from_tiles(tiles) + + def read_tiles(self, order: int, file_format: str): + properties = HipsSurveyProperties.read(self.base_path / "properties") + tile_path = self.base_path / f"Norder{order}" + + for path in tile_path.glob(f"Dir*/Npix*.{file_format}"): + ipix = int(path.as_posix().split("/")[-1].split(".")[0][4:]) + frame = properties.hips_frame + width = properties.tile_width + meta = HipsTileMeta(order, ipix, file_format, frame, width) + yield HipsTile.read(meta, path) diff --git a/hips/draw/tests/test_healpix.py b/hips/draw/tests/test_healpix.py index 0a16ab4..31882cd 100644 --- a/hips/draw/tests/test_healpix.py +++ b/hips/draw/tests/test_healpix.py @@ -1,61 +1,75 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst -import pytest import numpy as np -from PIL import Image -from astropy.io import fits -from numpy.testing import assert_allclose, assert_equal +from numpy.testing import assert_equal from astropy_healpix import healpy as hp -from ..healpix import healpix_to_hips, healpix_to_hips_tile +from ..healpix import healpix_to_hips_tile -def test_healpix_to_hips_tile(): - nside, tile_width = 4, 2 - npix = hp.nside2npix(nside) - hpx_data = np.arange(npix, dtype="uint8") +def make_hpx_data(file_format): + npix = hp.nside2npix(4) + data = np.arange(npix, dtype="uint8") + + if file_format == "fits": + return data + elif file_format == "jpg": + return np.moveaxis([data, data + 1, data + 2], 0, -1) + elif file_format == "png": + return np.moveaxis([data, data + 1, data + 2, data + 3], 0, -1) + else: + raise ValueError() + + +def test_healpix_to_hips_tile_fits(): + hpx_data = make_hpx_data("fits") + tile = healpix_to_hips_tile( - hpx_data=hpx_data, - tile_width=tile_width, - tile_idx=0, - file_format="fits", - frame="galactic", + hpx_data=hpx_data, tile_width=2, tile_idx=0, file_format="fits" ) - assert_equal(tile.data, [[1, 3], [0, 2]]) assert tile.meta.order == 1 assert tile.meta.ipix == 0 assert tile.meta.file_format == "fits" - assert tile.meta.frame == "galactic" + assert tile.meta.frame == "icrs" assert tile.meta.width == 2 + assert tile.data.dtype == np.uint8 + assert tile.data.shape == (2, 2) + assert_equal(tile.data, [[1, 3], [0, 2]]) + -@pytest.mark.parametrize("file_format", ["fits", "png"]) -def test_healpix_to_hips(tmpdir, file_format): - nside, tile_width = 4, 2 - npix = hp.nside2npix(nside) - hpx_data = np.arange(npix, dtype="uint8") +def test_healpix_to_hips_tile_jpg(): + hpx_data = make_hpx_data("jpg") - healpix_to_hips( - hpx_data=hpx_data, - tile_width=tile_width, - base_path=tmpdir, - file_format=file_format, - frame="galactic", + tile = healpix_to_hips_tile( + hpx_data=hpx_data, tile_width=2, tile_idx=0, file_format="jpg" ) - # The test data is filled with np.arange(), here we reproduce the sum of the - # indices in the nested scheme manually for comparison - desired = hpx_data.reshape((-1, tile_width, tile_width)) - - for idx, val in enumerate(desired): - filename = str(tmpdir / f"Norder1/Dir0/Npix{idx}.{file_format}") - if file_format is "fits": - data = fits.getdata(filename) - data = np.rot90(data, k=-1) - else: - data = np.array(Image.open(filename)) - data = data.T - assert_allclose(val, data) - - properties = (tmpdir / "properties").read_text(encoding=None) - assert file_format in properties - assert "galactic" in properties + assert tile.meta.order == 1 + assert tile.meta.ipix == 0 + assert tile.meta.file_format == "jpg" + assert tile.meta.width == 2 + + assert tile.data.dtype == np.uint8 + assert tile.data.shape == (2, 2, 3) + + # Note: we don't assert on the tile data for JPEG, + # because JPEG encoding noise is large and results + # can vary depending on JPEG lib and machine used. + + +def test_healpix_to_hips_tile_png(): + hpx_data = make_hpx_data("png") + + tile = healpix_to_hips_tile( + hpx_data=hpx_data, tile_width=2, tile_idx=0, file_format="png" + ) + + assert tile.meta.order == 1 + assert tile.meta.ipix == 0 + assert tile.meta.file_format == "png" + assert tile.meta.frame == "icrs" + assert tile.meta.width == 2 + + assert tile.data.dtype == np.uint8 + assert tile.data.shape == (2, 2, 4) + assert_equal(tile.data[..., 0], [[1, 3], [0, 2]]) diff --git a/hips/draw/tests/test_hipsgen.py b/hips/draw/tests/test_hipsgen.py new file mode 100644 index 0000000..224d8e7 --- /dev/null +++ b/hips/draw/tests/test_hipsgen.py @@ -0,0 +1,34 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +import pytest +import numpy as np +from numpy.testing import assert_equal +from ...tiles import HipsTile, HipsTileMeta +from ..hipsgen import healpix_to_hips +from .test_healpix import make_hpx_data + + +@pytest.mark.parametrize("file_format", ["fits", "png"]) +def test_healpix_to_hips(tmpdir, file_format): + hpx_data = make_hpx_data(file_format) + + healpix_to_hips( + hpx_data=hpx_data, tile_width=2, base_path=tmpdir, file_format=file_format + ) + + properties = (tmpdir / "properties").read_text(encoding=None) + assert file_format in properties + assert "icrs" in properties + + # Check one tile + filename = str(tmpdir / f"Norder1/Dir0/Npix2.{file_format}") + meta = HipsTileMeta(order=1, ipix=2, file_format=file_format, width=2) + tile = HipsTile.read(meta, filename) + + if file_format == "fits": + assert tile.data.dtype == np.uint8 + assert tile.data.shape == (2, 2) + assert_equal(tile.data, [[9, 11], [8, 10]]) + elif file_format == "png": + assert tile.data.dtype == np.uint8 + assert tile.data.shape == (2, 2, 4) + assert_equal(tile.data[..., 0], [[9, 11], [8, 10]]) diff --git a/hips/tiles/allsky.py b/hips/tiles/allsky.py index 829bbb5..ef9958b 100644 --- a/hips/tiles/allsky.py +++ b/hips/tiles/allsky.py @@ -92,9 +92,6 @@ 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