Download

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