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

4 files changed, 17 insertions, 24 deletions

diff --git a/bin/healpix2hpx b/bin/healpix2hpx
index c684e09..670e27a 100755
--- a/bin/healpix2hpx
+++ b/bin/healpix2hpx
@@ -28,8 +28,6 @@ def main():
     parser.add_argument("outfile", help="output FITS image in HPX projection")
     parser.add_argument("-C", "--clobber", action="store_true",
                         help="overwrite the existing output file")
-    parser.add_argument("-F", "--float", action="store_true",
-                        help="use float (single precision) instead of double")
     parser.add_argument("-l", "--log", dest="loglevel", default=None,
                         choices=["DEBUG", "INFO", "WARNING",
                                  "ERROR", "CRITICAL"],
@@ -68,9 +66,6 @@ def main():
     hpx_data, hpx_header = healpix2hpx(args.infile,
                                        append_history=history,
                                        append_comment=comments)
-    if args.float:
-        logger.info("HPX FITS images: use single-precision float numbers")
-        hpx_data = hpx_data.astype(np.float32)
     hdu = fits.PrimaryHDU(data=hpx_data, header=hpx_header)
     hdu.writeto(args.outfile, clobber=args.clobber, checksum=True)
     logger.info("HPX FITS images write to: %s" % args.outfile)
diff --git a/bin/hpx2healpix b/bin/hpx2healpix
index ac65fbd..ca485ac 100755
--- a/bin/hpx2healpix
+++ b/bin/hpx2healpix
@@ -13,9 +13,6 @@ import sys
 import argparse
 import logging
 
-import numpy as np
-from astropy.io import fits
-
 import fg21sim
 from fg21sim.configs import configs
 from fg21sim.utils import hpx2healpix, write_fits_healpix, setup_logging
@@ -28,8 +25,6 @@ def main():
     parser.add_argument("outfile", help="output HEALPix data file")
     parser.add_argument("-C", "--clobber", action="store_true",
                         help="overwrite the existing output file")
-    parser.add_argument("-F", "--float", action="store_true",
-                        help="use float (single precision) instead of double")
     parser.add_argument("-l", "--log", dest="loglevel", default=None,
                         choices=["DEBUG", "INFO", "WARNING",
                                  "ERROR", "CRITICAL"],
@@ -68,9 +63,6 @@ def main():
     hp_data, hp_header = hpx2healpix(args.infile,
                                      append_history=history,
                                      append_comment=comments)
-    if args.float:
-        logger.info("HEALPix data: use single-precision float numbers")
-        hp_data = hp_data.astype(np.float32)
     write_fits_healpix(args.outfile, hpmap=hp_data, header=hp_header,
                        clobber=args.clobber)
     logger.info("HEALPix data write to FITS file: %s" % args.outfile)
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):