#!/usr/bin/env python
#
version_num = "1.10" # Version number to be updated
"""
VERSION = '1.10';
author = 'Jose Vicente Perea';
date = '2021-02-12';
ChangeLog
=========
Version 1.10 - 2021-02-12 (JVP)
------------
+ FITS file not correctly closed
"fits.open()" replaced by content manager "with fits.open()"
Version 1.9 - 2017-05-10 (JVP)
------------
+ Time bin labels in panel plots corrected
Version 1.8 - 2017-05-09 (JVP)
------------
+ Major task upgrade:
- Re-binning of each exposure TimeSeries to get an optimal panel plot points density: pointsDen. Default: 512
- Final TS plots with the optimal time binning and the corresponding error bars (weight by FRACEXP)
- Re-binning function 'TSrebin' included in module 'ModuleUtil.py'
- Header documentation updated
Version 1.7 - 2017-04-03 (JVP)
------------
+ Correction: If output file name already has a valid extension (png or pdf): keep it as output file name
Version 1.6 - 2017-03-31 (JVP)
------------
+ Env shebang included '#!/usr/bin/env python'
+ If input file name already has a valid extension (png or pdf): keep it as output file name
+ Output figure plot size changed to 8.0 x 10.0 inches
Version 1.5 - 2017-03-24 (JVP)
------------
+ New time binning figured out to get a new_S/N = 3
Dt' = Dt * ( new_S/N / S/N )^1/2 # Assuming Counts_error ~ Counts ^ 1/2
new_S/N = 3
S/N = stddev(rate)/mean(rate_error)
Re-binned TS: Blue line for EPIC
Re-binned TS: Red and Black lines for RGS
+ Mean error rate is shown as a single error bar for each exposure
+ Reference time of all TimeSeries converted to UTC (astropy.time import Time)
+ Number of plots is always 4. If any TimeSeries is missing that plot will be empty
+ PN plot labeled as 'EPN Timing' if PN TS is Timing mode
+ Observation Target : OBJECT keyword removed from the plot
+ Ignoring warnings by using: warnings.catch_warnings()
+ Time binning removed as input option
Version 1.4 - 2016-12-21 (JVP)
------------
+ Continue if some of the keywords SRC_RA, SRC_DEC, SRCNUM do not exist
Version 1.3 - 2016-12-19 (JVP)
------------
+ Small correction in R1 or R2 conditional
Version 1.2 - 2016-12-15 (JVP)
------------
+ EPIC source id and coordinates labels only if they exist
Version 1.1 - 2016-12-02 (JVP)
------------
+ Include matplotlib.use('Agg') to avoid DISPLAY variable problems on PNG creation
Version 1.0 - 2016-11-21 (JVP)
------------
+ Initial version
Python script to plot a merge of Time Series in the same time frame
Original Time Series come from EPIC and/or RGS
------------------------------------------------------------------------------
"""
#
#
# BEFORE RUNNING:
#
# Point the PATH to the proper Python environment : 'ppsdev' or 'ppsprod':
#
# # export PATH=/lhome/pcms/local/anaconda3/envs/ppsdev/bin:$PATH
# # python plotTimeSeries.py -h
#
# ------------------------------------------------------------------------
#
# usage: plotTimeSeries [-h] [--version] [--infile INFILE] [--outfile OUTFILE] --outformat {png,pdf}
#
# Create a TimeSeries merged plot
#
# optional arguments:
# -h, --help show this help message and exit
# --version, -V show program's version number and exit
#
# --infile INFILE Input plain text file containing the TimeSeries file list. Default: TSFileList.asc
# --outfile OUTFILE Output figure file. Default: figure
# --outformat {png,pdf} Output figure format. Default: png
#### --bin BIN Time binning for the binned TimeSeries. Default: 1000 seconds --> REMOVED
#
# ------------------------------------------------------------------------
#
#
# print(__doc__) # Printout of the program documentation: Previous """ enclosed comments
import sys
import os
import argparse
import warnings
###################
# Inut Options
#
parser = argparse.ArgumentParser(description='Create a TimeSeries merged plot', prog='plotTimeSeries.py')
#parser.add_argument('--version', '-v', action='version', version='%(prog)s (%(prog)s-1.0) [-]')
parser.add_argument('--version', '-v', action='version', version='%(prog)s (%(prog)s-{}) [-]'.format(version_num))
parser.add_argument('--infile', help='Input plain text file containing the TimeSeries file list. Default: TSFileList.asc', default='TSFileList.asc')
parser.add_argument('--outfile', help='Output figure file. Default: figure. If the output filename already has a valid extension, it will be kept', default='figure')
parser.add_argument('--outformat', help='Output figure format. Default: png', choices=['png', 'pdf'], default='png')
parser.add_argument('--pointsDen', help='Points density. Default: 512 points in the whole plot', default=512, type=int)
args = parser.parse_args()
#
###################
#
import collections
#
#import re # Support for regular expressions (RE). e.g.: findall(pattern, string)
import numpy as np
from astropy import units as u
from astropy.io import ascii, fits
from astropy.coordinates import SkyCoord
from astropy.table import Table, join # Table: read, write, etc FITS files. join: joining Tables
from astropy.time import Time
import ModuleUtil # For function 'TSrebin'
import matplotlib
matplotlib.use('Agg') # To avoid problems with DISPLAY variable
# raise RuntimeError('Invalid DISPLAY variable')
# RuntimeError: Invalid DISPLAY variable
import matplotlib.pyplot as plt
plt.rcParams['font.size'] = 10.0
###################
# sys.exit(1)
###################
fileList = args.infile
ext = args.outformat
# If input file name already has a valid extension (png or pdf): keep it.
if args.outfile.lower().endswith(('.png', '.pdf')):
#print('\n ----> outfile name already has a valid format')
figure = args.outfile
else:
figure = args.outfile + '.' + ext
points_density = args.pointsDen
if os.path.isfile(fileList) and os.access(fileList, os.R_OK):
#print("Input file " + fileList + " exists and is readable")
pass
else:
print("Either file " + fileList + " is missing or is not readable\n")
sys.exit(1)
data = ascii.read(fileList, names = ['files'], comment='///' ,format='no_header')
plots_n = 4 # Fixed number of plots. Missing plot if there is no LC
#print()
###########################################################################################################################
#
# Read all TS's from files and get :
# times (all TSTART + all TSTOP)
# instruments, exposures
# tstarts() = 1st 'TIME' value [0]
# tstops() = last 'TIME' value [-1]
#
####################################################
times=[]
instruments=[]
exposures =[]
tstarts =[]
tstops =[]
c = 0
for i in data['files']:
# c += 1
file = i
#t = fits.open(file)
with fits.open(file) as t:
instrument = t[1].header['INSTRUME']
expid = t[0].header['EXPIDSTR']
tstart = t[1].header['TSTART']
tstop = t[1].header['TSTOP']
time_bin = t[1].header['TIMEDEL']
try:
timesys = t[1].header['TIMESYS']
timeref = t[1].header['TIMEREF']
mjdref = t[1].header['MJDREF']
except KeyError:
pass
times.append(tstart)
times.append(tstop)
instruments.append(instrument)
exposures.append(expid)
tstarts.append(t[1].data['TIME'][0])
tstops.append(t[1].data['TIME'][-1])
# Find the Optimail Time Bin for each Exposure
#
# Total Plot must have a points density == points_density. Default: 512
#
# Create a Table with: Instr., ExpID, time(0), time(-1)
#
instr_times = Table([instruments, exposures, tstarts, tstops], names=['instrument', 'exposures', 'tstart', 'tstop'])
instr_by_name = instr_times.group_by('instrument')
count_instr=collections.Counter( instr_times['instrument'] )
instr_Dtimes=[]
Dtimes=[]
for i in count_instr.keys(): # Run over the Instruments
#print(' instrument = ', i)
instr_mask = instr_by_name.groups.keys['instrument'] == i
instr_max_time = max(instr_by_name.groups[instr_mask]['tstop'])
instr_min_time = min(instr_by_name.groups[instr_mask]['tstart'])
instr_Dt = (instr_max_time - instr_min_time) / ( points_density - 1 ) # (Total Exposure time) / ( # points in plot )
#print(' instr_Dt = ', instr_Dt)
instr_Dtimes.append(i)
Dtimes.append(instr_Dt)
DTIMES = Table([instr_Dtimes, Dtimes], names=['instrument', 'Dtimes'])
# Create the final Table with: Instr., ExpID, time(0), time(-1), Optimal Time Bin (Dtimes)
#
# Joinning both tables
instr_deltaTimes = join(instr_times, DTIMES, keys='instrument')
### print(instr_deltaTimes) # New time bins table per Instrument and Exposure
max_DTout = max(instr_deltaTimes['Dtimes'])
###########################################################################################################################
maxtime = max(times)
mintime = min(times)
t_ref = mintime # Refernece Time = Minimum of the TSTART keywords
maxtime_ref = maxtime - t_ref
#
# Reference Time in UTC:
# Reference values from the last EPIC LC read in the previous loop
# timesys, timeref, mjdref = TT LOCAL 50814.0
if 'mjdref' in globals() and 'timesys' in globals():
t_ref_utc = Time(mjdref + t_ref/86400,format='mjd',scale=(timesys).lower(), precision=0).utc.iso # Conversion to UTC
print()
print('Minimum TSTART from all Timeseries: Reference Time (UTC) = ', t_ref_utc)
else:
t_ref_utc = t_ref # No UTC conversion
print('Minimum TSTART from all Timeseries: Reference Time could not be converted to UTC = ', t_ref)
#
#
#print('mintime, maxtime = ', mintime, maxtime)
print()
####################################################
##
## Main loop over the TS files:
##
####################################################
ax = ()
f, ax = plt.subplots(plots_n, 1, sharex=True)
f.suptitle('Merged TimeSeries', fontsize=12)
#f.set_size_inches(10,12) # Figure size [optional]
f.set_size_inches(8.0,10.0)
#f.set_size_inches(8.27,11.7) # Figure size [optional]. A4 = 8.27 x 11.7 inch
##f.set_size_inches(10,14.15) # Figure size [optional]. A4 rate
pn_timing = 0
chk_r1 = chk_r2 = 0
for i in data['files']:
c += 1
#file = dir + i
file = i
if file.find('PN') != -1 and file.find('SRCTSR0000.F') != -1:
pn_timing = 1
#t = fits.open(file)
with fits.open(file) as t:
instrument = t[1].header['INSTRUME']
obs_id = t[1].header['OBS_ID']
objectID = t[1].header['OBJECT']
if instrument == 'EPN' or instrument == 'EMOS1' or instrument == 'EMOS2':
try:
ra_deg = t[1].header['SRC_RA']
dec_deg = t[1].header['SRC_DEC']
srcnum = t[1].header['SRCNUM']
coord = SkyCoord(ra=ra_deg*u.degree, dec=dec_deg*u.degree)
coords_hmsdms = coord.to_string('hmsdms')
except KeyError:
pass
expid = t[0].header['EXPIDSTR']
tstart = t[1].header['TSTART']
tstop = t[1].header['TSTOP']
time_bin = t[1].header['TIMEDEL']
print('File: ' , c , ' ' + file + ' ' + instrument )
if 'coords_hmsdms' in globals():
print(' RA,Dec: ', coords_hmsdms)
instruments.append(instrument)
n_instr = len(set(instruments)) # set ==> unique
#----------------------------------------------------------------------------
#
# Get the Optimail Time Bin for each Exposure from Table: 'instr_deltaTimes'
#
# If the Number of TimeSeries points is < , then the original time bin is kept
#
obs_by_name = instr_deltaTimes.group_by('instrument')
mask = obs_by_name.groups.keys['instrument'] == instrument
DT = obs_by_name.groups[mask]['Dtimes'][0]
new_time_bin = DT
#t_bin_frac = int(np.trunc( new_time_bin / time_bin )) # DTout / DTin ('TIMEDEL'). Re-binning factor. Trunc = round to the minor integer
t_bin_frac = int(np.round( new_time_bin / time_bin )) # DTout / DTin ('TIMEDEL'). Re-binning factor. Round = round to the nearest integer
# Read Input TS from file
with warnings.catch_warnings():
warnings.simplefilter("ignore", category=Warning) # To ignore: WARNING: UnitsWarning: 'cts' did not parse as fits unit
tsr=Table.read(file,hdu='RATE')
# DEBUG . Calculate the TS init time bin from the data
#print('==========================================')
#delta_times = []
#timeCol = tsr['TIME']
#for i in range(2, len(timeCol)):
# delta_times.append(timeCol[i] - timeCol[i-1])
#counter=collections.Counter(delta_times)
#print(counter)
#print('==========================================')
# DEBUG
if t_bin_frac <= 1: # ==> No re-binning
ntsr = tsr
ntsr['TIME'] = ntsr['TIME'] - t_ref
new_time_bin = time_bin
print(' No re-binning needed for this Light Curve. Input / Output time bins ratio : ', t_bin_frac)
else:
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
# Re-binning function call
#
ntsr = ModuleUtil.TSrebin(tsr, time_bin, new_time_bin) # import module 'ModuleUtil' required
ntsr['TIME'] = ntsr['TIME'] - t_ref
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#------------------------------------------------------------------------------------------------------------------------ #
# PLOTTING
#
if instrument == 'EPN':
pos = 0
if instrument == 'EMOS1':
pos = 1
if instrument == 'EMOS2':
pos = 2
if instrument == 'RGS1':
pos = 3
chk_r1 += 1
if instrument == 'RGS2':
pos = 3
chk_r2 += 1
xlimits = [0 - 10* max_DTout, maxtime_ref + 10 * max_DTout]
ymin = min( ntsr['RATE'] )
ymax = max( ntsr['RATE'] )
ax[pos].set_xlim( xlimits )
if instrument == 'RGS2':
color='b'
else:
color='r'
ax[pos].errorbar( ntsr['TIME'], ntsr['RATE'], ntsr['ERROR'], fmt='.', c=color, ms=3, ecolor='k')
if instrument == 'EPN' or instrument == 'EMOS1' or instrument == 'EMOS2':
ax[pos].annotate('time bin: {0:.2f} sec'.format(new_time_bin), xy=(0.01, 0.03), xycoords='axes fraction')
if instrument == 'RGS1':
if chk_r1 == 1:
ax[pos].annotate('time bin: {0:.2f} sec'.format(new_time_bin), xy=(0.01, 0.03), xycoords='axes fraction')
if instrument == 'RGS2':
if chk_r2 == 1:
ax[pos].annotate('time bin: {0:.2f} sec'.format(new_time_bin), xy=(0.75, 0.03), xycoords='axes fraction', color=color)
init_ylim = ax[pos].get_ylim()
if ymax > init_ylim[1]:
ax[pos].set_ylim(ymax=ymax)
if ymin < init_ylim[0]:
ax[pos].set_ylim(ymin=ymin)
if pn_timing == 1:
ax[0].set_title('EPN Timing', fontsize=10)
else:
ax[0].set_title('EPN', fontsize=10)
ax[1].set_title('EMOS1', fontsize=10)
ax[2].set_title('EMOS2', fontsize=10)
ax[3].set_title('RGS1+RGS2', fontsize=10)
plt.annotate('OBS_ID: {0}'.format(obs_id), xy=(0.13, 0.93), xycoords='figure fraction')
if 'srcnum' in globals():
plt.annotate('EPIC SRCNUM: {0}'.format(srcnum), xy=(0.58, 0.93), xycoords='figure fraction')
plt.annotate('RA,Dec: {0}'.format(coords_hmsdms), xy=(0.58, 0.91), xycoords='figure fraction')
plt.annotate("rate (count/s)", xy=(0.05, 0.6), xycoords='figure fraction', rotation=90)
plt.xlabel('time (s) after {0} (UTC)'.format(t_ref_utc))
#f.subplots_adjust(hspace=0.0)
plt.savefig(figure)