Skip to content

Commit d58a628

Browse files
author
Tom Barclay
committed
initial commit
1 parent a6e3d2b commit d58a628

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

111 files changed

+27472
-0
lines changed

‎PyKE.tar

1.15 MB
Binary file not shown.

‎addpath.py

+15
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,15 @@
1+
# Inserts the directory containing links
2+
# to all Python-based tasks in the STSDAS tree to
3+
# the default Python path.
4+
5+
import iraf,sys
6+
7+
# define path to top level Python directory in STSDAS
8+
9+
_path = iraf.osfn('kepler$')
10+
11+
# if directory not already in PYTHONPATH,...
12+
13+
if _path not in sys.path:
14+
# Add the directory to path.
15+
sys.path.insert(1,_path)

‎keparith.par

+12
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,12 @@
1+
infile,s,a,'kepler.fits',,,'Input file (string)'
2+
outfile,s,a,'keparith.fits',,,'Name of output FITS file (string)'
3+
datacol,s,a,'SAP_FLUX',,,'Name of the column containing the flux time series (string)'
4+
constantfunc,s,a,'None','None|median|mean|MAD|std|max|range','A value calculated from a function to be used in operation (string)'
5+
constantval,r,a,None,,,'Use a given number in operation (float)'
6+
operation,s,a,'Add','add|subtract|multiply|divide','Divide flux by a constant (float)'
7+
clobber,b,h,'no',,,'Overwrite output file? (boolean)'
8+
verbose,b,h,'no',,,'Write to log file? (boolean)'
9+
logfile,s,h,'kepcotrend.log',,,'Name of ascii log file? (string)'
10+
status,i,h,0,,,'Exit status (0=good)'
11+
mode,s,h,"al"
12+

‎keparith.py

+218
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,218 @@
1+
"""
2+
Name: keparith.py
3+
Written by: Tom Barclay
4+
Date released: Nov 10 2011
5+
6+
Changelog:
7+
1.0 released
8+
2.0 made it possible to run from the command line if there are arguements given to the call to keparith
9+
2.1 fixed a bug when dividing a number and using PDCSAP_DATA - thanks to Joseph Long for spotting this
10+
"""
11+
12+
13+
14+
import sys
15+
import pyfits
16+
import matplotlib.pyplot as plt
17+
import kepio, kepmsg, kepkey, kepfit, kepstat
18+
from numpy import median,subtract,maximum,ones,multiply,float32,shape,absolute,mean,std,isfinite,where,nan
19+
20+
__svnid__ = "$Id: keparith.py 4326 2012-07-02 21:49:01Z tsbarcl2 $"
21+
__url__ = "$URL: svn+ssh://mstill@murzim.amn.nasa.gov/data-repo/trunk/data/flight/go/PyKE/kepler/keparith.py $"
22+
23+
24+
def MAD(xx,minSd=1E-16):
25+
"""Median Absolute Deviation
26+
"""
27+
med=median(xx,0)
28+
absdev=absolute(subtract(xx,med))
29+
mad=median(absdev,0)
30+
mad=maximum(mad,multiply(ones(mad.shape,float32),
31+
(minSd/1.48)))
32+
return mad
33+
34+
35+
def kepaddconstant(infile,outfile,datacol,constant,constantval,sign
36+
,clobber,verbose,logfile,status):
37+
# log the call
38+
39+
hashline = '----------------------------------------------------------------------------'
40+
kepmsg.log(logfile,hashline,verbose)
41+
call = 'KEPARITH -- '
42+
call += 'infile='+infile+' '
43+
call += 'outfile='+outfile+' '
44+
call += 'datacol='+str(datacol)+' '
45+
call += 'constant='+str(constant)+' '
46+
call += 'constantval='+str(constantval)+' '
47+
call += 'sign'+str(sign)+' '
48+
overwrite = 'n'
49+
if (clobber): overwrite = 'y'
50+
call += 'clobber='+overwrite+ ' '
51+
chatter = 'n'
52+
if (verbose): chatter = 'y'
53+
call += 'verbose='+chatter+' '
54+
call += 'logfile='+logfile
55+
kepmsg.log(logfile,call+'\n',verbose)
56+
57+
# start time
58+
59+
kepmsg.clock('KEPARITH started at',logfile,verbose)
60+
61+
# test log file
62+
63+
logfile = kepmsg.test(logfile)
64+
65+
# clobber output file
66+
67+
if clobber: status = kepio.clobber(outfile,logfile,verbose)
68+
if kepio.fileexists(outfile):
69+
message = 'ERROR -- KEPARITH: ' + outfile + ' exists. Use clobber=yes'
70+
status = kepmsg.err(logfile,message,verbose)
71+
72+
## open input file
73+
74+
instr, status = kepio.openfits(infile,'readonly',logfile,verbose)
75+
76+
if status == 0:
77+
try:
78+
test = str(instr[0].header['FILEVER'])
79+
version = 2
80+
except KeyError:
81+
version = 1
82+
83+
84+
# if version == 1:
85+
# #lc_flux = instr[1].data.field(datacol) + constant
86+
# lc_flux = instr[1].data.field(datacol)
87+
# instr[1].data.field(datacol)[:] = (lc_flux - median(lc_flux)) / MAD(lc_flux)
88+
# elif version == 2:
89+
# #lc_flux = instr[1].data.field(datacol) + constant
90+
# lc_flux = instr[1].data.field(datacol)
91+
# instr[1].data.field(datacol)[:] = (lc_flux - median(lc_flux)) / MAD(lc_flux)
92+
93+
lc_flux = instr[1].data.field(datacol)
94+
95+
if datacol == 'SAP_FLUX':
96+
errcol = 'SAP_FLUX_ERR'
97+
try:
98+
lc_err = instr[1].data.field(errcol)
99+
haveerr = True
100+
except:
101+
haveerr = False
102+
elif datacol == 'PDCSAP_FLUX':
103+
errcol = 'PDCSAP_FLUX_ERR'
104+
try:
105+
lc_err = instr[1].data.field(errcol)
106+
haveerr = True
107+
except:
108+
haveerr = False
109+
else:
110+
errcol = datacol + '_ERR'
111+
try:
112+
lc_err = instr[1].data.field(errcol)
113+
haveerr = True
114+
except:
115+
try:
116+
errcol = datacol + '_err'
117+
lc_err = instr[1].data.field(errcol)
118+
haveerr = True
119+
except:
120+
haveerr = False
121+
122+
#subtractor he just refers to the number that will be added/subtracted
123+
#divided or multiplied
124+
if isinstance(constantval,(long,int,float)) and constant == 'None':
125+
subtractor = float(constantval)
126+
elif constant.lower() == 'median':
127+
subtractor = float(median(lc_flux[isfinite(lc_flux)]))
128+
elif constant.lower() == 'mean':
129+
subtractor = float(mean(lc_flux[isfinite(lc_flux)]))
130+
elif constant.lower() == 'std':
131+
subtractor = float(std(lc_flux[isfinite(lc_flux)]))
132+
elif constant.lower() == 'mad':
133+
subtractor = float(MAD(lc_flux[isfinite(lc_flux)]))
134+
elif constant.lower() == 'max':
135+
subtractor = float(max(lc_flux[isfinite(lc_flux)]))
136+
elif constant.lower() == 'range':
137+
subtractor = float(max(lc_flux[isfinite(lc_flux)]) - min(lc_flux[isfinite(lc_flux)]))
138+
elif str(constant).lower() == 'none':
139+
subtractor = 0.
140+
message = 'No operation will be performed if you select None for the function and do not give a constant value'
141+
status = kepmsg.err(logfile,message,verbose)
142+
else:
143+
message = 'Your constant term is not in the list of possible functions'
144+
status = kepmsg.err(logfile,message,verbose)
145+
146+
if subtractor == 0. and sign == 'divide' and status == 0:
147+
message = 'You are trying to divide by zero: not a good idea.'
148+
status = kepmsg.err(logfile,message,verbose)
149+
150+
151+
152+
if status == 0:
153+
if sign.lower() == 'add':
154+
instr[1].data.field(datacol)[:] = where(isfinite(instr[1].data.field(datacol)[:]),(lc_flux + subtractor),nan)
155+
elif sign.lower() == 'subtract':
156+
instr[1].data.field(datacol)[:] = where(isfinite(instr[1].data.field(datacol)[:]),(lc_flux - subtractor),nan)
157+
elif sign.lower() == 'divide':
158+
instr[1].data.field(datacol)[:] = where(isfinite(instr[1].data.field(datacol)[:]),(lc_flux / subtractor),nan)
159+
if haveerr:
160+
instr[1].data.field(errcol)[:] = where(isfinite(instr[1].data.field(errcol)[:]),(lc_err / subtractor),nan)
161+
elif sign.lower() == 'multiply':
162+
instr[1].data.field(datacol)[:] = where(isfinite(instr[1].data.field(datacol)[:]),(lc_flux * subtractor),nan)
163+
if haveerr:
164+
instr[1].data.field(errcol)[:] = where(isfinite(instr[1].data.field(errcol)[:]),(lc_err * subtractor),nan)
165+
else:
166+
message = 'Your operation need to be one of: add, subtract, divide or multiply'
167+
status = kepmsg.err(logfile,message,verbose)
168+
169+
if status == 0:
170+
instr.writeto(outfile)
171+
172+
173+
if status == 0:
174+
status = kepio.closefits(instr,logfile,verbose)
175+
176+
## end time
177+
178+
if (status == 0):
179+
message = 'KEPARITH completed at'
180+
else:
181+
message = '\nKEPARITH aborted at'
182+
kepmsg.clock(message,logfile,verbose)
183+
184+
# main
185+
186+
if '--shell' in sys.argv:
187+
import argparse
188+
189+
parser = argparse.ArgumentParser(description='Time invariant algebra on light curve data')
190+
parser.add_argument('--shell', action='store_true', help='Are we running from the shell?')
191+
192+
parser.add_argument('infile', help='Name of input file', type=str)
193+
parser.add_argument('outfile', help='Name of FITS file to output', type=str)
194+
195+
parser.add_argument('--datacol', '-d', default='SAP_FLUX', dest='datacol', help='Name of the column containing the flux time series', type=str)
196+
parser.add_argument('--constantfunc', '-f', default='None', dest='constantfunc', help='A value calculated from a function to be used in operation', type=str,
197+
choices=['None','median','mean','MAD','std','max','range'])
198+
parser.add_argument('--constantval', '-v', default='None', dest='constantval', help='Use a given number in operation', type=str)
199+
parser.add_argument('--operation', '-o', default='add', dest='operation', help='Add, subtract, multiply, divide flux by a constant', type=str,
200+
choices=['add','','subtract','multiply','divide'])
201+
202+
parser.add_argument('--clobber', action='store_true', help='Overwrite output file?')
203+
parser.add_argument('--verbose', action='store_true', help='Write to a log file?')
204+
parser.add_argument('--logfile', '-l', help='Name of ascii log file', default='kepcotrend.log', dest='logfile', type=str)
205+
parser.add_argument('--status', '-e', help='Exit status (0=good)', default=0, dest='status', type=int)
206+
207+
208+
args = parser.parse_args()
209+
210+
kepaddconstant(args.infile,args.outfile,args.datacol,args.constantfunc,args.constantval,args.operation,args.clobber,args.verbose,args.logfile,args.status)
211+
212+
213+
else:
214+
from pyraf import iraf
215+
parfile = iraf.osfn("kepler$keparith.par")
216+
t = iraf.IrafTaskFactory(taskname="keparith", value=parfile, function=kepaddconstant)
217+
218+

‎keparray.py

+105
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,105 @@
1+
import numpy as n
2+
import scipy.interpolate
3+
import scipy.ndimage
4+
5+
def rebin2D(a, newdims, method='linear', centre=False, minusone=False):
6+
'''Arbitrary resampling of source array to new dimension sizes.
7+
Currently only supports maintaining the same number of dimensions.
8+
To use 1-D arrays, first promote them to shape (x,1).
9+
10+
Uses the same parameters and creates the same co-ordinate lookup points
11+
as IDL''s congrid routine, which apparently originally came from a VAX/VMS
12+
routine of the same name.
13+
14+
method:
15+
neighbour - closest value from original data
16+
nearest and linear - uses n x 1-D interpolations using
17+
scipy.interpolate.interp1d
18+
(see Numerical Recipes for validity of use of n 1-D interpolations)
19+
spline - uses ndimage.map_coordinates
20+
21+
centre:
22+
True - interpolation points are at the centres of the bins
23+
False - points are at the front edge of the bin
24+
25+
minusone:
26+
For example- inarray.shape = (i,j) & new dimensions = (x,y)
27+
False - inarray is resampled by factors of (i/x) * (j/y)
28+
True - inarray is resampled by(i-1)/(x-1) * (j-1)/(y-1)
29+
This prevents extrapolation one element beyond bounds of input array.
30+
'''
31+
if not a.dtype in [n.float64, n.float32]:
32+
a = n.cast[float](a)
33+
34+
m1 = n.cast[int](minusone)
35+
ofs = n.cast[int](centre) * 0.5
36+
old = n.array( a.shape )
37+
ndims = len( a.shape )
38+
if len( newdims ) != ndims:
39+
print "[congrid] dimensions error. " \
40+
"This routine currently only support " \
41+
"rebinning to the same number of dimensions."
42+
return None
43+
newdims = n.asarray( newdims, dtype=float )
44+
dimlist = []
45+
46+
if method == 'neighbour':
47+
for i in range( ndims ):
48+
base = n.indices(newdims)[i]
49+
dimlist.append( (old[i] - m1) / (newdims[i] - m1) \
50+
* (base + ofs) - ofs )
51+
cd = n.array( dimlist ).round().astype(int)
52+
newa = a[list( cd )]
53+
return newa
54+
55+
elif method in ['nearest','linear']:
56+
# calculate new dims
57+
for i in range( ndims ):
58+
base = n.arange( newdims[i] )
59+
dimlist.append( (old[i] - m1) / (newdims[i] - m1) \
60+
* (base + ofs) - ofs )
61+
# specify old dims
62+
olddims = [n.arange(i, dtype = n.float) for i in list( a.shape )]
63+
64+
# first interpolation - for ndims = any
65+
mint = scipy.interpolate.interp1d( olddims[-1], a, kind=method )
66+
newa = mint( dimlist[-1] )
67+
68+
trorder = [ndims - 1] + range( ndims - 1 )
69+
for i in range( ndims - 2, -1, -1 ):
70+
newa = newa.transpose( trorder )
71+
72+
mint = scipy.interpolate.interp1d( olddims[i], newa, kind=method )
73+
newa = mint( dimlist[i] )
74+
75+
if ndims > 1:
76+
# need one more transpose to return to original dimensions
77+
newa = newa.transpose( trorder )
78+
79+
return newa
80+
elif method in ['spline']:
81+
oslices = [ slice(0,j) for j in old ]
82+
oldcoords = n.ogrid[oslices]
83+
nslices = [ slice(0,j) for j in list(newdims) ]
84+
newcoords = n.mgrid[nslices]
85+
86+
newcoords_dims = range(n.rank(newcoords))
87+
#make first index last
88+
newcoords_dims.append(newcoords_dims.pop(0))
89+
newcoords_tr = newcoords.transpose(newcoords_dims)
90+
# makes a view that affects newcoords
91+
92+
newcoords_tr += ofs
93+
94+
deltas = (n.asarray(old) - m1) / (newdims - m1)
95+
newcoords_tr *= deltas
96+
97+
newcoords_tr -= ofs
98+
99+
newa = scipy.ndimage.map_coordinates(a, newcoords)
100+
return newa
101+
else:
102+
print "Congrid error: Unrecognized interpolation type.\n", \
103+
"Currently only \'neighbour\', \'nearest\',\'linear\',", \
104+
"and \'spline\' are supported."
105+
return None

‎kepbin.par

+18
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,18 @@
1+
infile,s,a,'kepler.fits',,,'Input file (string)'
2+
outfile,s,a,'kepbin.fits',,,'Name of output FITS file (string)'
3+
fluxcol,s,a,'SAP_FLUX',,,'Column containing the flux (string)'
4+
do_nbin,b,a,'no',,,'Specify number of bins? (boolean)'
5+
nbins,i,h,None,,,'Number of bins (integer)'
6+
do_binwidth,b,a,'no',,,'Specify binwidth? (boolean)'
7+
binwidth,r,h,None,,,'bin width (float)'
8+
do_ownbin,b,a,'no',,,'Specify ownbins? (boolean)'
9+
ownbins,s,h,'',,,'own bins file (string)'
10+
method,s,h,'linear','linear|nearest|zero|slinear|quadratic|cubic','method of interpolation (string)'
11+
interpm,s,h,'quad','quad|romberg','method of integration (string)'
12+
plot,b,a,'no',,,'Plot result? (boolean)'
13+
clobber,b,h,'no',,,'Overwrite output file? (boolean)'
14+
verbose,b,h,'no',,,'Write to log file? (boolean)'
15+
logfile,s,h,'kepclip.log',,,'Name of ascii log file? (string)'
16+
status,i,h,0,,,'Exit status (0=good)'
17+
mode,s,h,"al"
18+

0 commit comments

Comments
 (0)