From 401f08d355aba8db82318adf88968b45b95fcdf7 Mon Sep 17 00:00:00 2001 From: Adeel Ahmad Date: Sat, 28 Jul 2018 16:06:30 +0500 Subject: [PATCH 1/3] Add function for HiPS to HEALPix conversion This function is replicated from: https://gist.github.com/tboch/f68cd1bb1529d8ac12184b40e54ba692#file-hips2healpix-py --- hips/draw/healpix.py | 120 +++++++++++++++++++++++++++++++- hips/draw/tests/test_healpix.py | 12 +++- 2 files changed, 128 insertions(+), 4 deletions(-) diff --git a/hips/draw/healpix.py b/hips/draw/healpix.py index f61ffbe..d10360a 100644 --- a/hips/draw/healpix.py +++ b/hips/draw/healpix.py @@ -1,13 +1,17 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst -from pathlib import Path import numpy as np -from astropy_healpix import healpy as hp +import healpy as hp +from pathlib import Path +from astropy.io import fits + from ..tiles import HipsTile, HipsTileMeta, HipsSurveyProperties from ..utils.healpix import hips_tile_healpix_ipix_array __all__ = [ 'healpix_to_hips_tile', 'healpix_to_hips', + 'hips_to_healpix_array', + 'hips_to_healpix' ] @@ -88,3 +92,115 @@ def healpix_to_hips(hpx_data, tile_width, base_path, file_format='fits'): properties = HipsSurveyProperties(data=data) properties.write(base_path / 'properties') + +def hips_to_healpix_array(order: int): + """Create array giving mapping between: + + HEALPix index ==> XY index + XY index ==> HEALPix index + + Parameters + ---------- + order : int + Order of HEALPix array. + """ + def fill_array(npix, nsize, pos): + size = nsize ** 2 + npix_idx = [[0 for i in range(size // 4)] for j in range(4)] + nb = [0 for i in range(4)] + + for i in range(size): + if (i % nsize) < (nsize // 2): + dg = 0 + else: + dg = 1 + + if i < (size / 2): + bh = 1 + else: + bh = 0 + + quad = (dg << 1) | bh + + if pos is None: + j = i + else: + j = pos[i] + + npix[j] = npix[j] << 2 | quad + npix_idx[quad][nb[quad]] = j + nb[quad] += 1 + + if size > 4: + for i in range(4): + fill_array(npix, nsize // 2, npix_idx[i]) + + if order == 0: + xy_to_healpix = healpix_to_xy = [0] + return + + nsize = 2 ** order + xy_to_healpix = [0 for i in range(nsize ** 2)] + healpix_to_xy = [0 for i in range(nsize ** 2)] + + fill_array(xy_to_healpix, nsize, None) + + for i in range(len(xy_to_healpix)): + healpix_to_xy[xy_to_healpix[i]] = i + + return healpix_to_xy, xy_to_healpix + +def hips_to_healpix(hips_url: str, npix: int, hpx_output_path: str): + """Given a HiPS survey, generate a HEALPix map. + + Parameters + ---------- + hips_url : str + URL of HiPS survey. + npix : int + HEALPix pixel number. + hpx_output_path : str + HEALPix output path. + """ + # Query the HiPS survey properties. + hips_properties = HipsSurveyProperties.fetch(hips_url + "/properties") + + hips_tile_width = hips_properties.tile_width + order = int(np.log2(hips_tile_width)) + + healpix_to_xy, xy_to_healpix = hips_to_healpix_array(order) + + healpix_map = np.empty(hp.nside2npix(2 ** 12), dtype=np.double) # 2 ** 12? + healpix_map[:] = hp.pixelfunc.UNSEEN + + for ipix in range(npix): + try: + tile = fits.open(f'{hips_url}/Norder3/Dir0/Npix{ipix}.fits') + except: + continue + + # BSCALE and BZERO. + bzero = 0 + bscale = 1 + has_blank_value = False + + if 'BZERO' in tile[0].header: + bzero = tile[0].header['BZERO'] + if 'BSCALE' in tile[0].header: + bscale = tile[0].header['BSCALE'] + if 'BLANK' in tile[0].header: + has_blank_value = True + blank_value_scaled = bscale * tile[0].header['BLANK'] + bzero + + for k in range(hips_tile_width ** 2): + x = k % hips_tile_width + y = k // hips_tile_width + healpix_index = xy_to_healpix[k] + pixel_value = tile[0].data[y][x] + + if has_blank_value and np.fabs(pixel_value - blank_value_scaled) < 1e-6: + pixel_value = hp.pixelfunc.UNSEEN + + healpix_map[ipix * hips_tile_width * hips_tile_width + healpix_index] = pixel_value + + hp.write_map(str(hpx_output_path), healpix_map, coord='C', nest=True) diff --git a/hips/draw/tests/test_healpix.py b/hips/draw/tests/test_healpix.py index eb04767..968e457 100644 --- a/hips/draw/tests/test_healpix.py +++ b/hips/draw/tests/test_healpix.py @@ -1,11 +1,12 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst import pytest import numpy as np +import healpy as hp from PIL import Image from astropy.io import fits from numpy.testing import assert_allclose -from astropy_healpix import healpy as hp -from ..healpix import healpix_to_hips + +from ..healpix import healpix_to_hips, hips_to_healpix @pytest.mark.parametrize('file_format', ['fits', 'png']) @@ -36,3 +37,10 @@ def test_healpix_to_hips(tmpdir, file_format): properties = (tmpdir / 'properties').read_text(encoding=None) assert file_format in properties + +def test_hips_to_healpix(tmpdir): + hips_to_healpix(hips_url='http://alasky.u-strasbg.fr/Pan-STARRS/DR1/g', + npix=768, + hpx_output_path=tmpdir / 'panstarrs-g.hpx') + + healpix_map = hp.read_map(str(tmpdir / 'panstarrs-g.hpx')) From 18f0f6e103c069aafa793d5785e26a1adaadd22f Mon Sep 17 00:00:00 2001 From: Adeel Ahmad Date: Sat, 28 Jul 2018 16:20:21 +0500 Subject: [PATCH 2/3] Add mypy annotations in healpix_to_hips_tile and healpix_to_hips function --- hips/draw/healpix.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/hips/draw/healpix.py b/hips/draw/healpix.py index d10360a..6e74515 100644 --- a/hips/draw/healpix.py +++ b/hips/draw/healpix.py @@ -2,6 +2,7 @@ import numpy as np import healpy as hp from pathlib import Path +from typing import Union from astropy.io import fits from ..tiles import HipsTile, HipsTileMeta, HipsSurveyProperties @@ -15,7 +16,8 @@ ] -def healpix_to_hips_tile(hpx_data, tile_width, tile_idx, file_format) -> HipsTile: +def healpix_to_hips_tile(hpx_data: np.ndarray, tile_width: int, + tile_idx: int, file_format: str) -> HipsTile: """Create single hips tile from healpix data given a tile index. Parameters @@ -59,7 +61,8 @@ def healpix_to_hips_tile(hpx_data, tile_width, tile_idx, file_format) -> HipsTil return HipsTile.from_numpy(meta=meta, data=data) -def healpix_to_hips(hpx_data, tile_width, base_path, file_format='fits'): +def healpix_to_hips(hpx_data: np.ndarray, tile_width: int, + base_path: Union[str, Path], file_format='fits') -> None: """Convert HEALPix image to HiPS. Parameters From b50d4db0e3065ab5834345ec61f08349fbb730e6 Mon Sep 17 00:00:00 2001 From: Adeel Ahmad Date: Wed, 22 Aug 2018 16:46:34 +0500 Subject: [PATCH 3/3] Update hips_to_healpix function to take HiPS tiles as the parameter HEALPix pixels are now passed as a NumPy array. The test case is updated to covert HEALPix to HiPS and then convert the resulting HiPS tiles to a HEALPix map. --- hips/draw/healpix.py | 112 +++++++------------------------- hips/draw/tests/test_healpix.py | 22 +++++-- 2 files changed, 39 insertions(+), 95 deletions(-) diff --git a/hips/draw/healpix.py b/hips/draw/healpix.py index 6e74515..3b9dca4 100644 --- a/hips/draw/healpix.py +++ b/hips/draw/healpix.py @@ -1,9 +1,10 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst import numpy as np -import healpy as hp from pathlib import Path -from typing import Union +from typing import Union, List from astropy.io import fits +from astropy_healpix import healpy as hp +import healpy as hp from ..tiles import HipsTile, HipsTileMeta, HipsSurveyProperties from ..utils.healpix import hips_tile_healpix_ipix_array @@ -11,7 +12,6 @@ __all__ = [ 'healpix_to_hips_tile', 'healpix_to_hips', - 'hips_to_healpix_array', 'hips_to_healpix' ] @@ -96,74 +96,22 @@ def healpix_to_hips(hpx_data: np.ndarray, tile_width: int, properties = HipsSurveyProperties(data=data) properties.write(base_path / 'properties') -def hips_to_healpix_array(order: int): - """Create array giving mapping between: - - HEALPix index ==> XY index - XY index ==> HEALPix index - - Parameters - ---------- - order : int - Order of HEALPix array. - """ - def fill_array(npix, nsize, pos): - size = nsize ** 2 - npix_idx = [[0 for i in range(size // 4)] for j in range(4)] - nb = [0 for i in range(4)] - - for i in range(size): - if (i % nsize) < (nsize // 2): - dg = 0 - else: - dg = 1 - - if i < (size / 2): - bh = 1 - else: - bh = 0 - - quad = (dg << 1) | bh - - if pos is None: - j = i - else: - j = pos[i] - - npix[j] = npix[j] << 2 | quad - npix_idx[quad][nb[quad]] = j - nb[quad] += 1 - - if size > 4: - for i in range(4): - fill_array(npix, nsize // 2, npix_idx[i]) - - if order == 0: - xy_to_healpix = healpix_to_xy = [0] - return - - nsize = 2 ** order - xy_to_healpix = [0 for i in range(nsize ** 2)] - healpix_to_xy = [0 for i in range(nsize ** 2)] - - fill_array(xy_to_healpix, nsize, None) - - for i in range(len(xy_to_healpix)): - healpix_to_xy[xy_to_healpix[i]] = i - - return healpix_to_xy, xy_to_healpix - -def hips_to_healpix(hips_url: str, npix: int, hpx_output_path: str): +def hips_to_healpix(hips_url: str, hips_tiles: List[HipsTile], healpix_pixels: np.ndarray) -> np.ndarray: """Given a HiPS survey, generate a HEALPix map. Parameters ---------- + hips_tiles : List[HipsTile] + List of HiPS tiles. hips_url : str URL of HiPS survey. - npix : int - HEALPix pixel number. - hpx_output_path : str - HEALPix output path. + healpix_pixels : `~numpy.ndarray` + HEALPix pixel numbers. + + Returns + ------- + healpix_map : `~numpy.ndarray` + HEALPix map. """ # Query the HiPS survey properties. hips_properties = HipsSurveyProperties.fetch(hips_url + "/properties") @@ -171,39 +119,25 @@ def hips_to_healpix(hips_url: str, npix: int, hpx_output_path: str): hips_tile_width = hips_properties.tile_width order = int(np.log2(hips_tile_width)) - healpix_to_xy, xy_to_healpix = hips_to_healpix_array(order) + # xy_to_healpix = hips_tile_healpix_ipix_array(shift_order=order).flatten() healpix_map = np.empty(hp.nside2npix(2 ** 12), dtype=np.double) # 2 ** 12? healpix_map[:] = hp.pixelfunc.UNSEEN - for ipix in range(npix): - try: - tile = fits.open(f'{hips_url}/Norder3/Dir0/Npix{ipix}.fits') - except: - continue - - # BSCALE and BZERO. - bzero = 0 - bscale = 1 + for index, tile in enumerate(hips_tiles): has_blank_value = False - if 'BZERO' in tile[0].header: - bzero = tile[0].header['BZERO'] - if 'BSCALE' in tile[0].header: - bscale = tile[0].header['BSCALE'] - if 'BLANK' in tile[0].header: + if 'BLANK' in tile.data: has_blank_value = True blank_value_scaled = bscale * tile[0].header['BLANK'] + bzero - for k in range(hips_tile_width ** 2): - x = k % hips_tile_width - y = k // hips_tile_width - healpix_index = xy_to_healpix[k] - pixel_value = tile[0].data[y][x] + # healpix_index = xy_to_healpix[index] + pixel_value = tile.data - if has_blank_value and np.fabs(pixel_value - blank_value_scaled) < 1e-6: - pixel_value = hp.pixelfunc.UNSEEN + if has_blank_value and np.fabs(pixel_value - blank_value_scaled) < 1e-6: + pixel_value = hp.pixelfunc.UNSEEN - healpix_map[ipix * hips_tile_width * hips_tile_width + healpix_index] = pixel_value + healpix_map_range = index * (hips_tile_width ** 2) + healpix_map[healpix_map_range: healpix_map_range + tile.meta.width] = pixel_value[0] - hp.write_map(str(hpx_output_path), healpix_map, coord='C', nest=True) + return healpix_map diff --git a/hips/draw/tests/test_healpix.py b/hips/draw/tests/test_healpix.py index 968e457..b970d8f 100644 --- a/hips/draw/tests/test_healpix.py +++ b/hips/draw/tests/test_healpix.py @@ -1,12 +1,13 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst import pytest import numpy as np -import healpy as hp from PIL import Image from astropy.io import fits +from astropy_healpix import healpy as hp from numpy.testing import assert_allclose -from ..healpix import healpix_to_hips, hips_to_healpix +from ..healpix import healpix_to_hips, healpix_to_hips_tile +from ..healpix import hips_to_healpix @pytest.mark.parametrize('file_format', ['fits', 'png']) @@ -39,8 +40,17 @@ def test_healpix_to_hips(tmpdir, file_format): assert file_format in properties def test_hips_to_healpix(tmpdir): - hips_to_healpix(hips_url='http://alasky.u-strasbg.fr/Pan-STARRS/DR1/g', - npix=768, - hpx_output_path=tmpdir / 'panstarrs-g.hpx') + nside, tile_width = 4, 2 + npix = hp.nside2npix(nside) + healpix_data = np.arange(npix, dtype='uint8') + + n_tiles = healpix_data.size // tile_width ** 2 + + tiles = [] + for tile_idx in range(n_tiles): + tiles.append(healpix_to_hips_tile(hpx_data=healpix_data, tile_width=tile_width, + tile_idx=tile_idx, file_format='fits')) - healpix_map = hp.read_map(str(tmpdir / 'panstarrs-g.hpx')) + healpix_map = hips_to_healpix(hips_url='http://alasky.u-strasbg.fr/Pan-STARRS/DR1/g', + hips_tiles=tiles, + healpix_pixels=healpix_data)