diff options
| -rwxr-xr-x | bin/healpix2hpx | 5 | ||||
| -rwxr-xr-x | bin/hpx2healpix | 8 | ||||
| -rw-r--r-- | fg21sim/utils/fits.py | 19 | ||||
| -rw-r--r-- | fg21sim/utils/healpix.py | 9 |
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): |