From 4ce3c0b5ce2bd4b2d2bb89eaae669a4e3eba6662 Mon Sep 17 00:00:00 2001
From: Aaron LI <aaronly.me@outlook.com>
Date: Sun, 9 Oct 2016 22:31:36 +0800
Subject: utils: Add function "read_fits_healpix()"

This function wraps on the `healpy.read_map()`, but reset the data array
to its original dtype in FITS file, as well as return the FITS header in
`astropy.io.fits.Header` object.
---
 fg21sim/utils/__init__.py |  2 +-
 fg21sim/utils/fits.py     | 29 +++++++++++++++++++++++++++++
 2 files changed, 30 insertions(+), 1 deletion(-)

diff --git a/fg21sim/utils/__init__.py b/fg21sim/utils/__init__.py
index b6cf3fd..97243a7 100644
--- a/fg21sim/utils/__init__.py
+++ b/fg21sim/utils/__init__.py
@@ -1,6 +1,6 @@
 # Copyright (c) 2016 Weitian LI <liweitianux@live.com>
 # MIT license
 
-from .fits import write_fits_healpix
+from .fits import read_fits_healpix, write_fits_healpix
 from .healpix import healpix2hpx, hpx2healpix
 from .logging import setup_logging
diff --git a/fg21sim/utils/fits.py b/fg21sim/utils/fits.py
index e865528..be5dc68 100644
--- a/fg21sim/utils/fits.py
+++ b/fg21sim/utils/fits.py
@@ -9,6 +9,7 @@ from datetime import datetime, timezone
 
 import numpy as np
 from astropy.io import fits
+import healpy as hp
 
 
 # Reference:
@@ -26,6 +27,34 @@ FITS_COLUMN_FORMATS = {
 }
 
 
+def read_fits_healpix(filename):
+    """Read the HEALPix map from a FITS file to 1D array in *RING* ordering.
+
+    Parameters
+    ----------
+    filename : str
+        Filename of the HEALPix FITS file
+
+    Returns
+    -------
+    data : 1D `~numpy.ndarray`
+        HEALPix data in *RING* ordering
+    header : `~astropy.io.fits.Header`
+        Header of the input FITS file
+
+    NOTE
+    ----
+    This function wraps on `healpy.read_map()`, but set the data type of
+    data array to its original value as in FITS file, as well as return
+    FITS header as `~astropy.io.fits.Header` instance.
+    """
+    hdu = fits.open(filename)[0]
+    dtype = hdu.data.dtype
+    header = hdu.header
+    data = hp.read_map(hdu, nest=False, verbose=False)
+    return (data.astype(dtype), header)
+
+
 def write_fits_healpix(filename, hpmap, header=None, clobber=False):
     """Write the HEALPix map to a FITS file with proper header as well
     as the user-provided header.
-- 
cgit v1.2.2