aboutsummaryrefslogtreecommitdiffstats
path: root/astro/lc_clean.py
blob: a0fdc7b4ad7c8343f45aee6d59e85cbe65c818b1 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
#
# Aaron LI
# Created: 2016-01-16
# Updated: 2016-01-16
#

"""
Clean the lightcurve by fitting the RATE data with a Gaussian model,
and discard the time bins with RATE beyond [mean-n*sigma, mean+n*sigma].
"""

__version__ = "0.1.0"
__date__ = "2016-01-16"

import sys
import argparse

from astropy.io import fits
import numpy as np


class LightCurve:
    """
    X-ray data light curve class
    """
    def __init__(self, lcfile):
        f = fits.open(lcfile)
        self.lc_data = f[1].data
        self.lc_header = f[1].header
        self.time     = self.lc_data['TIME']
        self.rate     = self.lc_data['RATE']
        self.rate_err = self.lc_data['ERROR']
        self.TSTART   = self.lc_header['TSTART']
        self.TSTOP    = self.lc_header['TSTOP']
        self.TIMEDEL  = self.lc_header['TIMEDEL']
        self.TIMEPIXR = self.lc_header['TIMEPIXR']
        f.close()

    def sigma_clip(self, nsigma=3, maxiter=10):
        """
        Iteratively clip the time bins whose value lie beyond the
        range [mean-n*sigma, mean+n*sigma].
        """
        rate = self.rate
        keep_idx = np.ones(rate.shape, dtype=bool)  # all True's
        keep_num = np.sum(keep_idx)
        keep_num0 = np.inf
        i = 0
        while (keep_num < keep_num0):
            if (i >= maxiter):
                print("WARNING: maximum iteration limit reached",
                        file=sys.stderr)
                break
            keep_num0 = keep_num
            i += 1
            mean  = np.mean(rate[keep_idx])
            sigma = np.std(rate[keep_idx])
            cut_low  = mean - nsigma * sigma
            cut_high = mean + nsigma * sigma
            keep_idx = np.logical_and((rate >= cut_low), (rate <= cut_high))
            keep_num = np.sum(keep_idx)
        # save clip results
        self.niter = i
        self.keep_idx = keep_idx
        self.time_clipped = self.time[keep_idx]
        self.rate_clipped = self.rate[keep_idx]

    def make_gti(self, apply_header=True):
        """
        Make new GTIs (good time intervals) according to the clipped
        time bins.
        """
        frac = 0.01  # TIMEDEL fraction to distingush two time bins
        gti_start = []
        gti_stop  = []
        time_start = self.time_clipped
        time_stop  = time_start + self.TIMEDEL
        # first GTI start time
        gti_start.append(time_start[0])
        for tstart, tstop in zip(time_start[1:], time_stop[:-1]):
            if (np.abs(tstart-tstop) <= frac * self.TIMEDEL):
                # time bin continues
                continue
            else:
                # a new GTI start
                gti_start.append(tstart)
                gti_stop.append(tstop)
        # last GTI stop time
        gti_stop.append(time_stop[-1])
        # convert to numpy array
        gti_start = np.array(gti_start)
        gti_stop  = np.array(gti_stop)
        if apply_header:
            # add TSTART to the time
            gti_start += self.TSTART
            gti_stop  += self.TSTART
        # save results
        self.gti_start = gti_start
        self.gti_stop  = gti_stop

    def write_gti(self, filename=None, header=True):
        """
        Write generated GTIs to file or screen (default)
        """
        if isinstance(filename, str):
            outfile = open(filename, 'w')
        else:
            outfile = sys.stdout
        #
        if header:
            outfile.write('# TSTART\tTSTOP\n')
        outfile.write('\n'.join([ '%s\t%s' % (tstart, tstop) \
                for tstart, tstop in zip(self.gti_start, self.gti_stop) ]))
        #
        if isinstance(filename, str):
            outfile.close()


def main():
    parser = argparse.ArgumentParser(
            description="Clean light curve by sigma clipping")
    parser.add_argument("-V", "--version", action="version",
            version="%(prog)s " + "%s (%s)" % (__version__, __date__))
    parser.add_argument("infile",
            help="input lightcurve file; contains [TIME, RATE] columns")
    parser.add_argument("outfile", nargs='?', default=None,
            help="output text-format GTI file; for XSELECT filter time")
    parser.add_argument("-s", "--nsigma", dest="nsigma", type=float,
            default=2.0, help="sigma clipping significant level")
    parser.add_argument("-H", "--no-header", dest="noheader",
            action="store_true", help="not write header to the output file")
    parser.add_argument("-v", "--verbose", dest="verbose",
            action="store_true", help="show verbose information")
    args = parser.parse_args()

    lc = LightCurve(args.infile)
    lc.sigma_clip(nsigma=args.nsigma)
    lc.make_gti(apply_header=True)
    lc.write_gti(filename=args.outfile, header=(not args.noheader))
    if args.verbose:
        exposure = np.sum(lc.gti_stop - lc.gti_start)
        print("# Total GTI: %.2f (s)" % exposure)


if __name__ == "__main__":
    main()


#  vim: set ts=4 sw=4 tw=0 fenc= ft=python: #