From 2aaf7d8106e7d596654c9ae1cf0fa75fd642f3eb Mon Sep 17 00:00:00 2001 From: Aaron LI Date: Sun, 19 Feb 2017 16:26:40 +0800 Subject: Setup module 'acispy' and clean 'scripts' directory Setup a module 'acispy' to hold some generic Python modules for better/easier reuse, and clean up the 'scripts' directory, which will be used to hold the CLI tools. --- acispy/__init__.py | 0 acispy/acis.py | 91 +++++++++++++++++++ acispy/ds9.py | 30 +++++++ acispy/pfiles.py | 38 ++++++++ acispy/region.py | 251 +++++++++++++++++++++++++++++++++++++++++++++++++++++ acispy/spectrum.py | 44 ++++++++++ 6 files changed, 454 insertions(+) create mode 100644 acispy/__init__.py create mode 100644 acispy/acis.py create mode 100644 acispy/ds9.py create mode 100644 acispy/pfiles.py create mode 100644 acispy/region.py create mode 100644 acispy/spectrum.py (limited to 'acispy') diff --git a/acispy/__init__.py b/acispy/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/acispy/acis.py b/acispy/acis.py new file mode 100644 index 0000000..91d9fcf --- /dev/null +++ b/acispy/acis.py @@ -0,0 +1,91 @@ +# Copyright (c) 2017 Weitian LI +# MIT license + +""" +Chandra ACIS utilities +""" + +import math +import subprocess +import re + + +class ACIS: + """ + Chandra ACIS detector properties and utilities. + + References + ---------- + [1] CIAO - dictionary - ACIS (advanced camera for imaging and spectroscopy) + http://cxc.harvard.edu/ciao/dictionary/acis.html + [2] CIAO - Dictionary - PI (pulse invariant) + http://cxc.harvard.edu/ciao/dictionary/pi.html + """ + # Pixel size + pixel2arcsec = 0.492 # [arcsec] + # Number of channels + nchannel = 1024 + # Channel energy width + echannel = 14.6 # [eV] + + @classmethod + def energy2channel(self, energy): + """ + Convert energy [eV] to channel number. + """ + return math.floor(energy/self.echannel + 1) + + @classmethod + def get_type(self, filepath): + """ + Determine the Chandra ACIS type (``I`` or ``S``) according the + active ACIS chips. + + Parameters + ---------- + filepath : str + Path to the input FITS file + + Returns + ------- + acis_type : str + ``I`` if ACIS-I, ``S`` if ACIS-S; + otherwise, ``ValueError`` raised. + """ + subprocess.check_call(["punlearn", "dmkeypar"]) + detnam = subprocess.check_output([ + "dmkeypar", "infile=%s" % filepath, "keyword=DETNAM", "echo=yes" + ]).decode("utf-8").strip() + if re.match(r"^ACIS-0123", detnam): + return "I" + elif re.match(r"^ACIS-[0-6]*7", detnam): + return "S" + else: + raise ValueError("unknown chip combination: %s" % detnam) + + @classmethod + def get_chips_str(self, filepath, sep=":"): + """ + Return a string of the chips of interest according to the + active ACIS type. + + Parameters + ---------- + filepath : str + Path to the input FITS file + sep : str, optional + Separator to join the chip ranges, e.g., 0:3, 0-3 + + Returns + ------- + chips : str + ``0:3`` if ACIS-I, ``7`` if ACIS-S; + otherwise, ``ValueError`` raised. + """ + acis_type = self.get_type(filepath) + if acis_type == "I": + return sep.join(["0", "3"]) + elif acis_type == "S": + return "7" + else: + raise ValueError("unknown ACIS type") diff --git a/acispy/ds9.py b/acispy/ds9.py new file mode 100644 index 0000000..e8c8e88 --- /dev/null +++ b/acispy/ds9.py @@ -0,0 +1,30 @@ +# Copyright (c) 2017 Weitian LI +# MIT license + +""" +Wrapper function to view FITS files using DS9. +""" + +import subprocess + + +def ds9_view(filename, regfile=None, regformat="ciao", regsystem="physical", + cmap="he", binfactor=2, scale="linear", smooth=None): + """ + Wrapper function to view FITS files using DS9. + """ + cmd = [ + "ds9", filename, + "-cmap", cmap, + "-bin", "factor", str(binfactor), + "-scale", scale, + ] + if regfile: + cmd += [ + "-regions", "format", regformat, + "-regions", "system", regsystem, + "-regions", regfile, + ] + if smooth: + cmd += ["-smooth", "yes", "-smooth", "radius", str(smooth)] + subprocess.check_call(cmd) diff --git a/acispy/pfiles.py b/acispy/pfiles.py new file mode 100644 index 0000000..8a80a14 --- /dev/null +++ b/acispy/pfiles.py @@ -0,0 +1,38 @@ +# Copyright (c) 2017 Weitian LI +# MIT license +# +# Weitian LI +# 2017-02-06 + +""" +Prepare the CIAO parameter files and setup the PFILES environment +variable to keep the pfiles locally, in order to avoid the conflicts +between multiple instance of the same CIAO tools. +""" + +import os +import subprocess +import shutil + + +def setup_pfiles(tools): + """ + Copy the parameter files of the specified tools to the current + working directory, and setup the ``PFILES`` environment variable. + + Parameters + ---------- + tools : list[str] + Name list of the tools to be set up + """ + for tool in tools: + pfile = subprocess.check_output([ + "paccess", tool + ]).decode("utf-8").strip() + subprocess.check_call(["punlearn", tool]) + try: + shutil.copy(pfile, ".") + except shutil.SameFileError: + pass + # Setup the ``PFILES`` environment variable + os.environ["PFILES"] = "./:" + os.environ["PFILES"] diff --git a/acispy/region.py b/acispy/region.py new file mode 100644 index 0000000..09e42fa --- /dev/null +++ b/acispy/region.py @@ -0,0 +1,251 @@ +# Copyright (c) 2017 Weitian LI +# MIT license + +""" +Region utilities. + +NOTE: +Only support the commonly used region shapes in CIAO format. +While more complex shapes and DS9 format are not considered +at the moment. +""" + +import os +import re + + +class Point: + """ + Point region: point(xc,yc) + """ + shape = "point" + + def __init__(self, xc=None, yc=None): + if isinstance(xc, str): + self.parse(regstr=xc) + else: + self.xc = xc + self.yc = yc + + def parse(self, regstr): + g = re.split(r"[(),]", regstr) + self.xc = float(g[1]) + self.yc = float(g[2]) + + def __str__(self): + return "%s(%.4f,%.4f)" % (self.shape, self.xc, self.yc) + + +class Circle: + """ + Circle region: circle(xc,yc,radius) + """ + shape = "circle" + + def __init__(self, xc=None, yc=None, radius=None): + if isinstance(xc, str): + self.parse(regstr=xc) + else: + self.xc = xc + self.yc = yc + self.radius = radius + + def parse(self, regstr): + g = re.split(r"[(),]", regstr) + self.xc = float(g[1]) + self.yc = float(g[2]) + self.radius = float(g[3]) + + def __str__(self): + return "%s(%.4f,%.4f,%.4f)" % (self.shape, self.xc, self.yc, + self.radius) + + +class Annulus: + """ + Annulus region: + annulus(xc,yc,radius,radius2) + + xc, yc : center + + radius, radius2 : inner, outer radius + """ + shape = "annulus" + + def __init__(self, xc=None, yc=None, radius=None, radius2=None): + if isinstance(xc, str): + self.parse(regstr=xc) + else: + self.xc = xc + self.yc = yc + self.radius = radius + self.radius2 = radius2 + + def parse(self, regstr): + g = re.split(r"[(),]", regstr) + self.xc = float(g[1]) + self.yc = float(g[2]) + self.radius = float(g[3]) + self.radius2 = float(g[4]) + + def __str__(self): + return "%s(%.4f,%.4f,%.4f,%.4f)" % (self.shape, self.xc, self.yc, + self.radius, self.radius2) + + +class Pie: + """ + Pie region: + pie(xc,yc,radius,radius2,angle1,angle2) + + xc, yc : center + + radius, radius2 : inner, outer radius + + angle1, angle2: start, end angle [0, 360) + """ + shape = "pie" + + def __init__(self, xc=None, yc=None, radius=None, radius2=None, + angle1=0.0, angle2=360.0): + if isinstance(xc, str): + self.parse(regstr=xc) + else: + self.xc = xc + self.yc = yc + self.radius = radius + self.radius2 = radius2 + self.angle1 = angle1 + self.angle2 = angle2 + + def parse(self, regstr): + g = re.split(r"[(),]", regstr) + self.xc = float(g[1]) + self.yc = float(g[2]) + self.radius = float(g[3]) + self.radius2 = float(g[4]) + self.angle1 = float(g[5]) + self.angle2 = float(g[6]) + + def __str__(self): + return "%s(%.4f,%.4f,%.4f,%.4f,%.4f,%.4f)" % ( + self.shape, self.xc, self.yc, self.radius, self.radius2, + self.angle1, self.angle2) + + +class Ellipse: + """ + Ellipse region: + ellipse(xc,yc,radius,radius2,rotation) + + xc, yc : center + + radius, radius2 : semi-major / semi-minor axis + + rotation: rotation angle [0, 360) + """ + shape = "ellipse" + + def __init__(self, xc=None, yc=None, radius=None, radius2=None, + rotation=0.0): + if isinstance(xc, str): + self.parse(regstr=xc) + else: + self.xc = xc + self.yc = yc + self.radius = radius + self.radius2 = radius2 + self.rotation = rotation + + def parse(self, regstr): + g = re.split(r"[(),]", regstr) + self.xc = float(g[1]) + self.yc = float(g[2]) + self.radius = float(g[3]) + self.radius2 = float(g[4]) + self.rotation = float(g[5]) + + def __str__(self): + return "%s(%.4f,%.4f,%.4f,%.4f,%.4f)" % ( + self.shape, self.xc, self.yc, self.radius, + self.radius2, self.rotation) + + +class Box: + """ + Box region: box(xc,yc,width,height,rotation) + """ + shape = "box" + + def __init__(self, xc=None, yc=None, width=None, height=None, + rotation=0.0): + if isinstance(xc, str): + self.parse(regstr=xc) + else: + self.xc = xc + self.yc = yc + self.width = width + self.height = height + self.rotation = rotation + + def parse(self, regstr): + g = re.split(r"[(),]", regstr) + self.xc = float(g[1]) + self.yc = float(g[2]) + self.width = float(g[3]) + self.height = float(g[4]) + self.rotation = float(g[5]) + + def __str__(self): + return "%s(%.4f,%.4f,%.4f,%.4f,%.4f)" % ( + self.shape, self.xc, self.yc, self.width, + self.height, self.rotation) + + +class Regions: + """ + Manipulate the region files (CIAO format), as well as parse regions + from strings. + """ + REGION_SHAPES = { + "point": Point, + "circle": Circle, + "annulus": Annulus, + "pie": Pie, + "ellipse": Ellipse, + "box": Box, + } + + def __init__(self, regfile=None): + if regfile: + self.regions = self.load(regfile) + else: + self.regions = [] + + def load(self, infile): + regstr = [] + for line in open(infile): + if not (re.match(r"^\s*#.*$", line) or re.match(r"^\s*$", line)): + regstr.append(line.strip()) + self.regions = self.parse(regstr) + + def save(self, outfile, clobber=False): + if (not clobber) and os.path.exists(outfile): + raise OSError("output file already exists: %s" % outfile) + regstr = [str(reg) for reg in self.regions] + open(outfile, "w").write("\n".join(regstr) + "\n") + + @classmethod + def parse(self, regstr): + """ + Parse the given (list of) region string(s), and return the + corresponding parsed region objects. + """ + if isinstance(regstr, list): + return [self.parse_single(reg) for reg in regstr] + else: + return self.parse_single(regstr) + + @classmethod + def parse_single(self, regstr): + """ + Parse the given single region string to its corresponding + region object. + """ + shape = regstr.strip().split('(')[0].lower() + if shape in self.REGION_SHAPES: + return self.REGION_SHAPES[shape](regstr) + else: + raise ValueError("unknown region shape: %s" % shape) diff --git a/acispy/spectrum.py b/acispy/spectrum.py new file mode 100644 index 0000000..ad38b0d --- /dev/null +++ b/acispy/spectrum.py @@ -0,0 +1,44 @@ +# Copyright (c) 2017 Weitian LI +# MIT license + +""" +Chandra ACIS spectrum. +""" + + +from astropy.io import fits + +from acis import ACIS + + +class Spectrum: + """ + Chandra ACIS spectrum + """ + def __init__(self, filepath): + self.filepath = filepath + self.fitsobj = fits.open(filepath) + ext_spec = self.fitsobj["SPECTRUM"] + self.header = ext_spec.header + # spectral data + self.channel = ext_spec.data.columns["CHANNEL"].array + self.counts = ext_spec.data.columns["COUNTS"].array + # spectral keywords + self.EXPOSURE = self.header.get("EXPOSURE") + self.BACKSCAL = self.header.get("BACKSCAL") + + def calc_pb_flux(self, elow=9500, ehigh=12000): + """ + Calculate the particle background flux: + flux = counts / exposure / area + + Parameters + ---------- + elow, ehigh : float, optional + Lower and upper energy limit for the particle background. + """ + chlow = ACIS.energy2channel(elow) + chhigh = ACIS.energy2channel(ehigh) + counts = self.counts[(chlow-1):chhigh].sum() + flux = counts / self.EXPOSURE / self.BACKSCAL + return flux -- cgit v1.2.2