From c1c5b3d5e899f8e385c90aae4ff881ecad3c7422 Mon Sep 17 00:00:00 2001 From: Aaron LI 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