From c1c5b3d5e899f8e385c90aae4ff881ecad3c7422 Mon Sep 17 00:00:00 2001
From: Aaron LI <aaronly.me@outlook.com>
Date: Mon, 10 Oct 2016 00:29:31 +0800
Subject: utils: Preseve the dtype when read/write FITS files
* utils/fits.py: hack the dtype to ignore the byteorder (FITS data use
big endianness, e.g., dtype(">f4"))
* utils/healpix.py: explicit convert the dtype and log the dtype
* bin/healpix2hpx, bin/hpx2healpix: remove the --float argument
* other minor fixes/updates
---
fg21sim/utils/fits.py | 19 +++++++++++--------
fg21sim/utils/healpix.py | 9 ++++++---
2 files changed, 17 insertions(+), 11 deletions(-)
(limited to 'fg21sim/utils')
diff --git a/fg21sim/utils/fits.py b/fg21sim/utils/fits.py
index 365ed7e..96b5f58 100644
--- a/fg21sim/utils/fits.py
+++ b/fg21sim/utils/fits.py
@@ -53,8 +53,10 @@ def read_fits_healpix(filename):
if isinstance(filename, fits.BinTableHDU):
hdu = filename
else:
- hdu = fits.open(filename)[0]
- dtype = hdu.data.dtype
+ # Read the first extended table
+ hdu = fits.open(filename)[1]
+ # Hack to ignore the dtype byteorder, use native endianness
+ dtype = np.dtype(hdu.data.field(0).dtype.type)
header = hdu.header
data = hp.read_map(hdu, nest=False, verbose=False)
return (data.astype(dtype), header)
@@ -94,9 +96,12 @@ def write_fits_healpix(filename, hpmap, header=None, clobber=False):
- This function (currently) only implement the very basic feature of
the `healpy.write_map()`.
"""
- hpmap = np.array(hpmap)
+ hpmap = np.asarray(hpmap)
if hpmap.ndim != 1:
raise ValueError("Invalid HEALPix data: only support 1D array")
+ # Hack to ignore the dtype byteorder, use native endianness
+ dtype = np.dtype(hpmap.dtype.type)
+ hpmap = hpmap.astype(dtype)
#
npix = hpmap.size
nside = int((npix / 12) ** 0.5)
@@ -120,13 +125,11 @@ def write_fits_healpix(filename, hpmap, header=None, clobber=False):
#
hdr["EXTNAME"] = ("HEALPIX", "Name of the binary table extension")
hdr["CREATOR"] = (__name__, "File creator")
- hdr["DATE"] = (
- datetime.now(timezone.utc).astimezone().isoformat(),
- "File creation date"
- )
+ hdr["DATE"] = (datetime.now(timezone.utc).astimezone().isoformat(),
+ "File creation date")
# merge user-provided header
if header is not None:
- hdr.update(fits.Header(header))
+ hdr.extend(fits.Header(header))
#
hdu = fits.BinTableHDU.from_columns([
fits.Column(name="I", array=hpmap, format=colfmt)
diff --git a/fg21sim/utils/healpix.py b/fg21sim/utils/healpix.py
index 0307f08..ef1bd22 100644
--- a/fg21sim/utils/healpix.py
+++ b/fg21sim/utils/healpix.py
@@ -64,6 +64,8 @@ def healpix2hpx(data, header=None, append_history=None, append_comment=None):
else:
hp_data, hp_header = np.asarray(data), fits.header.Header(header)
logger.info("Read HEALPix data from array and header")
+ dtype = hp_data.dtype
+ logger.info("HEALPix dtype: {0}".format(dtype))
logger.info("HEALPix index ordering: %s" % hp_header["ORDERING"])
if hp_header["ORDERING"] != "RING":
raise ValueError("only 'RING' ordering currently supported")
@@ -72,7 +74,7 @@ def healpix2hpx(data, header=None, append_history=None, append_comment=None):
logger.info("HEALPix data: Npix=%d, Nside=%d" % (npix, nside))
if nside != hp_header["NSIDE"]:
raise ValueError("HEALPix data Nside does not match the header")
- hp_data = np.concatenate([hp_data, [np.nan]])
+ hp_data = np.append(hp_data, np.nan).astype(dtype)
hpx_idx = _calc_hpx_indices(nside)
# fix indices of "-1" to set empty pixels with above appended "nan"
hpx_idx[hpx_idx == -1] = len(hp_data) - 1
@@ -80,7 +82,7 @@ def healpix2hpx(data, header=None, append_history=None, append_comment=None):
hpx_header = _make_hpx_header(hp_header,
append_history=append_history,
append_comment=append_comment)
- return (hpx_data, hpx_header)
+ return (hpx_data.astype(hp_data.dtype), hpx_header)
def hpx2healpix(data, header=None, append_history=None, append_comment=None):
@@ -114,6 +116,7 @@ def hpx2healpix(data, header=None, append_history=None, append_comment=None):
else:
hpx_data, hpx_header = np.asarray(data), fits.header.Header(header)
logger.info("Read HPX image from array and header")
+ logger.info("HPX image dtype: {0}".format(hpx_data.dtype))
logger.info("HPX coordinate system: ({0}, {1})".format(
hpx_header["CTYPE1"], hpx_header["CTYPE2"]))
if ((hpx_header["CTYPE1"], hpx_header["CTYPE2"]) !=
@@ -136,7 +139,7 @@ def hpx2healpix(data, header=None, append_history=None, append_comment=None):
hp_header = _make_healpix_header(hpx_header, nside=nside,
append_history=append_history,
append_comment=append_comment)
- return (hp_data, hp_header)
+ return (hp_data.astype(hpx_data.dtype), hp_header)
def _calc_hpx_indices(nside):
--
cgit v1.2.2