Python HiPS package

This is the documentation page for hips. A Python astronomy package for HiPS : Hierarchical Progressive Surveys.

About

Description

HiPS (Hierarchical Progressive Surveys) is a way to store large astronomical survey sky image and catalog datasets on servers (such as HiPS at CDS), that allows clients to efficiently fetch only the image tiles or catalog parts for a given region of the sky they are interested in. Similar to Google maps, but for astronomy (see the HiPS paper).

This is a Python package to fetch and draw HiPS data.

It was just started in summer of 2017 and isn’t stable or feature complete yet. Feedback and contributions welcome!

Other resources

(If you have a HiPS-related webpage or tool or resource you’d like mentioned here, let us know!)

Thanks

This package is being developed as part of Google Summer of Code 2017 program by Adeel Ahmad, with Thomas Boch (CDS, Strasbourg) and Christoph Deil (MPIK, Heidelberg) as mentors. We would like to thank Google, CDS, MPIK for their support!

If you’re interested, you should follow Adeel’s blog: https://adl1995.github.io/

Also: thanks to the Astropy team for developing and maintaining the affiliated package-template and the ci-helpers! The recently introduced cookie-cutter makes it even quicker to set up a new package like this one in a good, maintainable way.

Installation

The hips package works with Python 3.6 or later, on Linux, MacOS and Windows.

Installing the latest stable version is possible either using pip or conda.

How to install the latest development version is desribed on the Develop page.

Using pip

To install hips with pip from PyPI, run:

pip install hips --no-deps

Note

The --no-deps flag is optional, but highly recommended if you already have Numpy installed, since otherwise pip will sometimes try to “help” you by upgrading your Numpy installation, which may not always be desired.

Using conda

To install hips with Anaconda from the conda-forge channel on anaconda.org simply run:

conda install -c conda-forge hips

Check installation

To check if you have hips installed, where it was installed and which version you have:

$ python
>>> import hips  # doctest: +SKIP
>>> hips.__version__  # doctest: +SKIP
# -> prints which version you have
>>> hips  # doctest: +SKIP
# -> prints where hips is installed

To see if you have the latest stable, released version of hips, you can find that version here:

Next you could try running the examples at Getting started and see if you get the expected output.

It’s usually not necessary, but if you find that your hips installation gives errors or unexpected results for examples that should work, you can run the hips automated tests via:

python -c 'import hips; hips.test()'

For more information on automated tests, see the Develop page.

Dependencies

The hips package has the following requirements:

  • Python 3.6 or later!
  • Numpy 1.11 or later
  • Astropy 1.2 or later
  • astropy-healpix 0.2 or later
  • scikit-image 0.12 or later. (Older versions could work, but aren’t tested.)
  • Pillow 4.0 or later. (Older versions could work, but aren’t tested.) Pillow is the friendly Python Imaging Library (PIL) fork, for JPEG and PNG tile I/O.

In addition, the following packages are needed for optional functionality:

  • Matplotlib 2.0 or later. Used for plotting in examples.
  • tqdm. Used for showing progress bar either on terminal or in Jupyter notebook.
  • aiohttp. Used for fetching HiPS tiles.

We have some info at Why only Python 3? on why we don’t support legacy Python (Python 2).

Getting started

This is a quick getting started guide for the Python hips package.

Make a sky image

To draw a sky image from HiPS image tiles with the hips package, follow the following three steps:

  1. Specify the sky image geometry you want by creating a WCSGeometry object:

    from astropy.coordinates import SkyCoord
from hips import WCSGeometry

geometry = WCSGeometry.create(
     skydir=SkyCoord(0, 0, unit='deg', frame='galactic'),
     width=2000, height=1000, fov="3 deg",
     coordsys='galactic', projection='AIT',
)

  2. Specify the HiPS survey you want. You just need to provide a valid HiPS survey ID.

    A good address that lists available HiPS data is http://aladin.u-strasbg.fr/hips/list

    hips_survey = 'CDS/P/DSS2/red'

  3. Call the make_sky_image function to fetch the HiPS data and draw it, returning an object of HipsDrawResult. By default a progress bar is shown for fetching and drawing HiPS tiles. For batch processing, this can be turned off by passing a keyword argument: progress_bar=False:

    from hips import make_sky_image
result = make_sky_image(geometry, hips_survey, 'fits')

Of course, you could change the parameters to chose any sky image geometry and available HiPS survey you like!

Work with the result

The HipsDrawResult object from the last section makes it easy for you to plot, save or analyse the sky image. To generate a quick-look plot of the sky image, with rectangles outlining the HiPS tiles that were fetched and drawn to create the sky image:

result.plot()

this will result in the following plot:

(Source code, png, hires.png, pdf)

_images/plot_fits.png

To save the sky image to a file:

result.write_image('my_image.fits')

To analyse the data, or make a publication-quality plot, you can get the sky image pixel data as a numpy.ndarray:

>>> result.image

and the sky image astropy.wcs.WCS mapping pixel to sky coordinates via:

>>> result.geometry.wcs

To print out summary information about the result:

>>> print(result)

The HipsDrawResult object also gives access to the HipsTile objects that were used for drawing the sky image, as well as other things.

Plot using Astropy visualization toolkit

Astropy provides a framework for plotting astronomical images with coordinates. It builds on top of Matplotlib and provides functionalities such as image normalization (scaling and stretching), smart histogram plotting, RGB color image creation from separate images. The framework also allows for customization of plotting styles.

The example below is for the FITS format and controls the stretch of the image through normalization. For FITS tiles, the data type is either int16 or float32:

import matplotlib.pyplot as plt
from astropy.visualization.mpl_normalize import simple_norm

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
im = ax.imshow(result.image, origin='lower', norm=norm, cmap='gray')
fig.colorbar(im)

(Source code, png, hires.png, pdf)

_images/plot_fits_astropy.png

RGB tiles can be plotted in much the same way as above, however, it is uncommon to apply an extra stretch in this case. For jpg and png tiles, the data type is uint8.

Note

For png tiles, there are four channel i.e. RGBA. The alpha channel is used for controlling the transparency of the image.

The example provided here is trivial. Astropy provides numerous other features, customizability options, and in-depth examples. Please see their documentation at: https://docs.astropy.org/en/stable/visualization

Make a color sky image

HiPS supports color images in jpg and png format. Making a color sky image works the same as the grayscale image example above, except that you get back a 3-dim Numpy array with (R, G, B) channels for jpg or (R, G, B, A) channels (A is transparency) for png.

Here’s an example using jpg and http://alasky.u-strasbg.fr/Fermi/Color:

(Source code, png, hires.png, pdf)

_images/plot_jpg.png

HiPS data

We plan to implement functionality to manage HiPS data, i.e. download it and cache it on a local disk. This isn’t available yet, at the moment we simply use Python lists of HipsTile objects, which have a read method for a given filename and a fetch method for a given URL.

More advanced examples

This package is under heavy development, it’s changing quickly.

We’ll add advanced examples and detailed documentation once things have stabilised a bit.

For now, if you know Python, you can look at the code and tests to see what’s available: https://github.com/hipspy/hips

What next?

That’s it, now you’ve seen the main features of the hips package. Note that there is API documentation explaining all available functions, classes and parameters.

If you have any questions, or find something not working or a missing feature, please get in touch by posting on our Github issue tracker.

Reference/API

hips Package

A Python astronomy package for HiPS : Hierarchical Progressive Surveys.

At the moment a client for HiPS images, but other contributions (HiPS catalogs or HiPS image generation) welcome!

Functions

fetch_tiles(tile_metas, hips_survey, …) Fetch a list of HiPS tiles.
healpix_to_hips(hpx_data, tile_width, base_path) Convert HEALPix image to HiPS.
healpix_to_hips_tile(hpx_data, tile_width, …) Create single hips tile from healpix data given a tile index.
make_sky_image(geometry, …) Make sky image: fetch tiles and draw.
test([package, test_path, args, plugins, …]) Run the tests using py.test.

Classes

HipsDrawResult(image, geometry, tile_format, …) HiPS draw result object (sky image and more).
HipsPainter(geometry, …) Paint a sky image from HiPS image tiles.
HipsSurveyProperties(data) HiPS properties container.
HipsSurveyPropertiesList(data) HiPS survey properties list.
HipsTile(meta, raw_data) HiPS tile container.
HipsTileAllskyArray(meta, raw_data) All-sky tile array container.
HipsTileMeta(order, ipix, file_format, …) HiPS tile metadata.
WCSGeometry(wcs, width, height) Sky image geometry: WCS and image shape.

Changelog

0.3 (unreleased)

  • no changes yet

0.2

Version 0.2 of hips was released on October 28, 2017.

  • Change from using healpy to astropy-healpix. This means hips now works on Windows! [#109]
  • Introduce asynchronous fetching of HiPS tiles [#106]
  • Add progress bar support for fetching and drawing HiPS tiles [#105]
  • Add reporting functionality for HipsPainter [#104]
  • Fix tile splitting criterion [#101]

0.1

This first version of the hips package was released on July 28, 2017. It contains a first implementation to fetch and draw tiles.

This is a very early release, to get some users and feedback. Note that the API will change in the coming weeks, and you can also expect new features, fixes and performance and usability enhancements.

The hips package started as a project developed as part of Google summer of code 2017, i.e. planning in early 2017 and coding in June 2017.

Develop

Hello!

Want to contribute to the hips package?

Great! Talk to us by filing a Github issue any time (it doesn’t have to be a concrete feature request or bug report).

This package was created using the Astropy affiliated package template, and everything works pretty much as in Astropy and most affiliated packages.

We didn’t write any developer docs specifically for this package yet. For now, check out the Astropy core package developer docs, or just talk to us if you have any questions.

Install development version

Install the latest development version from https://github.com/hipspy/hips :

git clone https://github.com/hipspy/hips
cd hips
pip install .

Then run the tests with either of these commands:

python -m pytest -v hips
python setup.py test -V

To run all tests and get a coverage report:

python setup.py test -V --remote-data --coverage

To build the documentation, do:

python setup.py build_docs

Get the hips-extra test datasets

To run tests accessing files from hips-extra repository:

git clone https://github.com/hipspy/hips-extra.git

Developers must have it cloned on their system, otherwise some test cases will be skipped. This contains tiles from different HiPS surveys which are used by the drawing module. After this, the HIPS_EXTRA environment variable must be set up on their system. On UNIX operating systems, this can be set using

export HIPS_EXTRA=path/to/hips-extra

Why only Python 3?

This package requires Python 3.6 or later.

It will not work with Python 2.7 or 3.5!

This was a choice we made when starting this package in summer of 2017, at a time when Jupyter had just made their Python 3 only release and other packages we depend on (like Astropy) were about to drop support for legacy Python (Python 2).

Supporting only Python 3 means we e.g. get these benefits:

  • async / await for asynchronous HiPS tile fetching (introduced in Python 3.5)
  • Keyword-only arguments (introduced in Python 3.0)
  • Type annotations (some only introduced in Python 3.6)
  • f-strings (introduced in Python 3.6)

At the moment, the only Python 3.6 only feature we use are f-strings, so if several potential users that are on Python 3.5 and can’t easily upgrade for some reason complain, we might consider supporting Python 3.5 in the next release.

HiPS tile drawing

This section describes the HiPS tile drawing algorithm implemented in this package, to create a sky image for a given WCS.

The original description was for the algorithm implemented in Aladin Lite, written by Thomas Boch. In the meantime, the algorithm implemented in this Python package has deviated a bit, they are no longer the same.

The implementation is based on numpy, astropy, healpy and scikit-image.

Naive algorithm

This is a naive (one could also say: simple and fast) algorithm for drawing HiPS tiles using affine transformations, implemented in the HipsPainter class and usually executed by users via the high-level make_sky_image function.

First we compute and fetch the tiles that are needed for the given sky image:

  1. The user specifies a WCSGeometry, which is a astropy.wcs.WCS as well as a width and height of the sky image to compute.
  2. Compute HiPS order corresponding to the requested image size/resolution. The attributes of HiPS properties needed for this are hips_order (order at the tile level) and hips_tile_width (number of pixels for tile width and height). If hips_tile_width is missing, a default value of 512 is assumed.
  3. Compute the list of tiles corresponding to the image FoV. This is done by computing the HiPS tile HEALPix pixel index for every pixel in the sky image and then computing the unique set.
  4. Fetch (HTTP calls or from a local cache) all tiles in the list.

Then we draw the tiles one by one using these steps:

  1. For each tile, compute the world coordinates of the tile corner vertices, using healpy boundaries function.
  2. For each tile, project vertices in image coordinates according to the requested WCS (performing ICRS to Galactic frame transformation if the HiPS and sky image aren’t in the same frame already).
  3. We extend the tile by 1 pixel in all directions in order to hide “stitches” with other tiles drawing (TODO: not done yet. needed?)
  4. The four corners uniquely define a projective transform between pixel coordinates on the tile and the output sky image. We use scikit-image to compute and apply that transform, which uses cubic spline interpolation under the hood. Thus the output is always float data, even if the input was integer RGB image data.

At the moment, we simply initialise the output sky image with pixel values of zero, and then sum the sky images we get from each projected tile. This is inefficient, and can result in incorrect pixel values at the lines corresponding to tile borders. We plan to implement a better (more efficient and more correct) way to handle that soon.

Note that any algorithm using interpolation is not fully conserving flux or counts. This might be a concern if you use the resulting sky images for data analysis. It’s your responsibility to decide if using this method is appropriate for your application or not!

Tile distortion issue

While the algorithm previously described works fine for HiPS tiles not distorted, it brings some astrometry offsets for distorted tiles. This distortion is strongly visible in the HEALPix scheme for tiles at the boundary between the equatorial zone and the polar cap.

An example of such distortions is shown in the example below (uncheck Activate deformations reduction algorithm to view the astrometry offsets): http://cds.unistra.fr/~boch/AL/test-reduce-deformations2.html

To overcome this problem, Aladin Desktop and Aladin Lite use the following recursive strategy: for tiles either too large (one edge is >300 pixels when projected) or too distorted (ratio of smaller diagonal on larger diagonal is smaller than 0.7 and one of the diagonal is > 150 pixels when projected):

  • We consider 4 children tiles, dynamically generated from the pixels of their father. Each children tile has a width and height equal to half of its father’s width/height.
  • For each children tile, we compute the world coordinates of its vertices, project them and either draw it if not too distorted or repeat the process by splitting again into 4 children tiles.

The recursion is limited by a maximum number of recursive steps (for 512x512 tiles, you are limited to a maximum of 9 steps as 2^9=512) and/or a maximum order (maximum order set arbitrarily at 19 in Aladin Desktop).

Precise algorithm

Note

The precise algorithm isn’t implemented yet.

Contrary the previous algorithm which used affine transformations, the idea here for the drawing step is to scan the result image pixels, and for each of them interpolate (Nearest neighbour or bilinear) the value, ie compute the indexes of nearest neighbour(s), retrieve the values of these pixels and merge them to determine the value of the target pixel. This is very similar to what reproject is doing.

One challenge is that one needs to know how to find the tile and pixel corresponding to a given HEALPix index. The correspondance between a HEALPix index and a pixel in a HiPS tile is given by a hpx2xy array (see method createHpx2xy in class cds.tools.pixtools.Util from Aladin Desktop source code.)

WCS for FITS tiles

It seems that the astrometry of a HiPS tile can be accurately described using a WCS header like this one (example for HiPS in equatorial frame, Norder 3, Npix 448):

NAXIS   =                    2 / number of data axes
NAXIS1  =                  512 / length of data axis 1
NAXIS2  =                  512 / length of data axis 1
CRPIX1  =              -2047.5 / Coordinate reference pixel
CRPIX2  =              -5631.5 / Coordinate reference pixel
CD1_1   = -1.0986328125000E-02 / Transformation matrix (rot + scale)
CD1_2   = -1.0986328125000E-02 / Transformation matrix (rot + scale)
CD2_1   =  1.0986328125000E-02 / Transformation matrix (rot + scale)
CD2_2   = -1.0986328125000E-02 / Transformation matrix (rot + scale)
CTYPE1  = 'RA---HPX'            / Longitude in an HPX projection
CTYPE2  = 'DEC--HPX'            /  Latitude in an HPX projection
CRVAL1  =                   0. / [deg] Longitude at the reference point
CRVAL2  =                   0. / [deg]  Latitude at the reference point
PV2_1   =                   4 / HPX H parameter (longitude)
PV2_2   =                   3 / HPX K parameter  (latitude)

HPX projection is supported by WCSLib. It is understood by DS9. Support in other tools (reproject, Montage, etc) is unclear and has to be tested.

Note

It seems that we can define a WCS for each tile. If so, this would allow us to simply use the reproject package to draw the tiles, which would be an alternative “precise” algorithm.