#!/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)