Add extract_sbprofile.py: Extract SBP from binned image

To replace 'ciao_sbp.sh' and 'chandra_ccdgap_rect.py'

1 files changed, 182 insertions, 0 deletions

diff --git a/bin/extract_sbprofile.py b/bin/extract_sbprofile.py
new file mode 100755
index 0000000..2f55002
--- /dev/null
+++ b/bin/extract_sbprofile.py
@@ -0,0 +1,182 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2016-2017 Weitian LI <liweitianux@live.com>
+# MIT license
+
+"""
+Extract surface brightness profile specified by the region file
+from a binned image file.
+
+NOTE
+----
+The CIAO official science thread "Obtain and Fit a Radial Profile" uses
+the (energy-filtered) event file to extract the SBP.  However, we found
+that it is more convenient and accurate to extract the SBP from a binned
+image instead of an event file (both are energy filtered of course),
+due to the region area calculation.
+
+To bin a image from the event file, the skyfov is usually provided to
+specify the chip boundaries, therefore, the created image has a size just
+enclosing the chips of interest (e.g., ACIS-7, or ACIS-0123).
+And internally, the tool ``dmcopy`` records the skyfov regions in the
+header keywords ``DSVAL?``, which are also used by ``dmcopy`` when
+excluding the detected point sources.
+
+when we provide the binned image to ``dmextract`` to extract SBP,
+the removed point source areas and chip boundaries will be considered,
+and the output SBP will have correct *effective* areas for each SBP region.
+(Also note that the given *exposure map* is also an image, which is
+similar to the input image.)
+
+Therefore, we can just draw a series of annuli (instead of complex pie
+regions), and ignore the chip boundaries as well.
+
+In addition, using an image as the input is also much faster, especially
+with large regions, or needing to consider excluding the outside-FoV
+region for using event file.
+
+
+References
+----------
+* CIAO: Ahelp: dmextract
+  http://cxc.harvard.edu/ciao/ahelp/dmextract.html
+* CIAO: Obtain and Fit a Radial Profile
+  http://cxc.harvard.edu/ciao/threads/radial_profile/
+"""
+
+import os
+import argparse
+import subprocess
+import logging
+
+import numpy as np
+from astropy.io import fits
+
+from _context import acispy
+from acispy.manifest import get_manifest
+from acispy.pfiles import setup_pfiles
+
+
+logging.basicConfig(level=logging.INFO)
+logger = logging.getLogger(__name__)
+
+
+def extract_sbp(outfile, infile, expmap, region, clobber=False):
+    """
+    Extract the surface brightness profile
+
+    If `bkg` is provided, then background subtraction is considered.
+    """
+    clobber = "yes" if clobber else "no"
+    logger.info("Extract SBP: %s[bin sky=@%s] ..." % (infile, region))
+    subprocess.check_call(["punlearn", "dmextract"])
+    subprocess.check_call([
+        "dmextract",
+        "infile=%s[bin sky=@%s]" % (infile, region),
+        "outfile=%s" % outfile, "exp=%s" % expmap,
+        "opt=generic", "clobber=%s" % clobber
+    ])
+    # add ``RMID`` and ``R_ERR`` columns
+    logger.info("Add 'RMID' and 'R_ERR' columns to SBP ...")
+    sbp_rmid = os.path.splitext(outfile)[0] + "_rmid.fits"
+    subprocess.check_call(["punlearn", "dmtcalc"])
+    subprocess.check_call([
+        "dmtcalc",
+        "infile=%s" % outfile, "outfile=%s" % sbp_rmid,
+        "expression=RMID=(R[0]+R[1])/2,R_ERR=(R[1]-R[0])/2",
+        "clobber=%s" % clobber
+    ])
+    os.rename(sbp_rmid, outfile)
+
+
+def make_sbp_txt(outfile, sbp, clobber=False):
+    """
+    Dump SBP data from FITS file to make a TXT data file,
+    which can be used to further create the QDP file, and be used
+    by the SBP fitting tools.
+    """
+    logger.info("Dump SBP data to text file: %s" % outfile)
+    with fits.open(sbp) as sbpfits:
+        rmid = sbpfits["HISTOGRAM"].data["RMID"]
+        r_err = sbpfits["HISTOGRAM"].data["R_ERR"]
+        sur_flux = sbpfits["HISTOGRAM"].data["SUR_FLUX"]
+        sur_flux_err = sbpfits["HISTOGRAM"].data["SUR_FLUX_ERR"]
+    sbpdata = np.column_stack([rmid, r_err, sur_flux, sur_flux_err])
+    np.savetxt(outfile, sbpdata,
+               header="RMID  R_ERR  SUR_FLUX  SUR_FLUX_ERR")
+
+
+def make_sbp_qdp(outfile, sbp_data, clobber=False):
+    """
+    Create a QDP file for the extracted SBP, for easier visualization.
+    """
+    logger.info("Create a QDP for SBP: %s" % outfile)
+    lines = [line.replace("#", "!") for line in open(sbp_data).readlines()]
+    lines = [
+        "READ SERR 1 2\n",
+        'LABEL T "Surface Brightness Profile"\n',
+        'LABEL X "Radius (pixel)"\n',
+        'LABEL Y "Surface Flux (photons/cm\\u2\\d/pixel\\u2\\d/s)"\n',
+        "LOG X Y ON\n"
+    ] + lines
+    with open(outfile, "w") as f:
+        f.writelines(lines)
+
+
+def main():
+    parser = argparse.ArgumentParser(
+            description="Extract surface brightness profile (SBP)")
+    parser.add_argument("-C", "--clobber", dest="clobber", action="store_true",
+                        help="overwrite existing file")
+    parser.add_argument("-r", "--region", dest="region",
+                        help="region file specifying the SBP " +
+                             "(default: 'sbp_reg' from manifest)")
+    parser.add_argument("-i", "--infile", dest="infile", required=True,
+                        help="input binned image ")
+    parser.add_argument("-e", "--expmap", dest="expmap",
+                        help="exposure map (default: 'expmap' from manifest)")
+    parser.add_argument("-o", "--outfile", dest="outfile",
+                        help="output SBP filename (default: same " +
+                        "basename as the input region file)")
+    args = parser.parse_args()
+
+    setup_pfiles(["dmextract", "dmtcalc"])
+
+    manifest = get_manifest()
+    if args.region:
+        region = args.region
+    else:
+        region = manifest.getpath("sbp_reg", relative=True)
+    if args.expmap:
+        expmap = args.expmap
+    else:
+        expmap = manifest.getpath("expmap", relative=True)
+    if args.outfile:
+        outfile = args.outfile
+    else:
+        outfile = os.path.splitext(region)[0] + ".fits"
+    sbp_data = os.path.splitext(outfile)[0] + ".txt"
+    sbp_qdp = os.path.splitext(outfile)[0] + ".qdp"
+
+    extract_sbp(outfile=outfile, infile=args.infile, expmap=expmap,
+                region=region, clobber=args.clobber)
+    make_sbp_txt(outfile=sbp_data, sbp=outfile, clobber=args.clobber)
+    make_sbp_qdp(outfile=sbp_qdp, sbp_data=sbp_data, clobber=args.clobber)
+
+    logger.info("Add SBP files to manifest ...")
+    key = "sbp"
+    manifest.setpath(key, outfile)
+    logger.info("Added '%s' to manifest: %s" % (key, manifest.get(key)))
+    key = "sbp_reg"
+    manifest.setpath(key, region)
+    logger.info("Added '%s' to manifest: %s" % (key, manifest.get(key)))
+    key = "sbp_data"
+    manifest.setpath(key, sbp_data)
+    logger.info("Added '%s' to manifest: %s" % (key, manifest.get(key)))
+    key = "sbp_qdp"
+    manifest.setpath(key, sbp_qdp)
+    logger.info("Added '%s' to manifest: %s" % (key, manifest.get(key)))
+
+
+if __name__ == "__main__":
+    main()