aboutsummaryrefslogtreecommitdiffstats
path: root/astro/21cm/extract_slice.py
diff options
context:
space:
mode:
Diffstat (limited to 'astro/21cm/extract_slice.py')
-rwxr-xr-xastro/21cm/extract_slice.py84
1 files changed, 0 insertions, 84 deletions
diff --git a/astro/21cm/extract_slice.py b/astro/21cm/extract_slice.py
deleted file mode 100755
index 576b7ee..0000000
--- a/astro/21cm/extract_slice.py
+++ /dev/null
@@ -1,84 +0,0 @@
-#!/usr/bin/env python3
-#
-# Copyright (c) 2017 Weitian LI <weitian@aaronly.me>
-# MIT License
-#
-
-"""
-Extract a slice from the 21cm cube and save as FITS image.
-"""
-
-import os
-import sys
-import argparse
-
-import numpy as np
-import astropy.io.fits as fits
-
-
-def main():
- outfile_default = "{prefix}_z{z:05.2f}_N{Nside}_L{Lside}_s{sidx}.fits"
-
- parser = argparse.ArgumentParser(
- description="Extract a slice from cube and save as FITS image")
- parser.add_argument("-C", "--clobber", dest="clobber",
- action="store_true",
- help="overwrite existing files")
- parser.add_argument("-d", "--dtype", dest="dtype", default="float32",
- help="NumPy dtype of data cubes (default: float32)")
- parser.add_argument("-z", "--redshift", dest="redshift",
- type=float, required=True,
- help="redshift of the input data cube")
- parser.add_argument("-L", "--len-side", dest="Lside",
- type=float, required=True,
- help="Side length of the cube [comoving Mpc]")
- parser.add_argument("-s", "--slice-idx", dest="sidx",
- type=int, default=None,
- help="slice index to be extracted (default: "
- "the central slice)")
- parser.add_argument("-u", "--unit", dest="unit",
- help="data unit (e.g., K, mK)")
- parser.add_argument("-i", "--infile", dest="infile", required=True,
- help="input data cube")
- parser.add_argument("-o", "--outfile", dest="outfile",
- default=outfile_default,
- help="output FITS image slice (default: %s)" %
- outfile_default)
- parser.add_argument("-p", "--prefix", dest="prefix", required=True,
- help="prefix for the output FITS image")
- args = parser.parse_args()
-
- cube = np.fromfile(open(args.infile, "rb"), dtype=args.dtype)
- Nside = round(cube.shape[0] ** (1.0/3))
- print("Read cube: %s (Nside=%d)" % (args.infile, Nside))
- if args.sidx is None:
- sidx = int(Nside / 2.0)
- elif args.idx >= 0 and args.idx < Nside:
- sidx = args.idx
- else:
- raise ValueError("invalid slice index: %s" % args.sidx)
- outfile = args.outfile.format(prefix=args.prefix, z=args.redshift,
- Nside=Nside, Lside=args.Lside, sidx=sidx)
- if os.path.exists(outfile) and not args.clobber:
- raise OSError("output file already exists: %s" % outfile)
-
- cube = cube.reshape((Nside, Nside, Nside))
- simg = cube[:, :, sidx]
- header = fits.Header()
- header["REDSHIFT"] = args.redshift
- header["Lside"] = (args.Lside, "Cube side length [comoving Mpc]")
- header["Nside"] = (Nside, "Number of pixels on each cube side")
- header["SliceIdx"] = (sidx, "Index of this extracted slice")
- if args.unit:
- header["BUNIT"] = (args.unit, "Data unit")
- header.add_history(" ".join(sys.argv))
- hdu = fits.PrimaryHDU(data=simg, header=header)
- try:
- hdu.writeto(outfile, overwrite=args.clobber)
- except TypeError:
- hdu.writeto(outfile, clobber=args.clobber)
- print("Extracted #%d slice: %s" % (sidx, outfile))
-
-
-if __name__ == "__main__":
- main()