Skip to content
Snippets Groups Projects
Commit 0e129342 authored by Vlad Kondratiev's avatar Vlad Kondratiev
Browse files

added aditional info fields to the diagnostics plot (dmc, rmc, polc, elevation, azimuth, transit)

parent 45ae96d5
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
l2com/lofar2-psr-quickview.py
+ 60
− 34
View file @ 0e129342
...@@ -27,7 +27,7 @@ import matplotlib.cm as cm ...@@ -27,7 +27,7 @@ import matplotlib.cm as cm
import matplotlib.colors as colors import matplotlib.colors as colors
import matplotlib.ticker as ticker import matplotlib.ticker as ticker
from matplotlib.ticker import * from matplotlib.ticker import *
from datetime import datetime from datetime import datetime, timedelta
import getpass import getpass
import argparse import argparse
...@@ -200,12 +200,13 @@ if __name__ == "__main__": ...@@ -200,12 +200,13 @@ if __name__ == "__main__":
default="undef", type=str) default="undef", type=str)
parser.add_argument("--pipeid", help="pipeline ID, mostly to be used by the pulp2-cwl-folding pipeline, for info purposes", parser.add_argument("--pipeid", help="pipeline ID, mostly to be used by the pulp2-cwl-folding pipeline, for info purposes",
default="undef", type=str) default="undef", type=str)
parser.add_argument("--paz-rfi-fraction", help="RFI fraction in % calculated by 'paz -r', for info purposes", parser.add_argument("--paz-rfi-fraction", help="RFI fraction in %% calculated by 'paz -r', for info purposes",
default="undef", type=str) default="undef", type=str)
parser.add_argument("--clfd-rfi-fraction", help="RFI fraction in % calculated by 'clfd', for info purposes", parser.add_argument("--clfd-rfi-fraction", help="RFI fraction in %% calculated by 'clfd', for info purposes",
default="undef", type=str) default="undef", type=str)
parser.add_argument("-v", "--verbose", help="more verbose output", parser.add_argument("-o", "--outstem", help="the stem for the output png file, by default the input file name is used without extension",
action="store_true") default="", type=str)
parser.add_argument("-v", "--verbose", help="more verbose output", action="store_true")
args = parser.parse_args() args = parser.parse_args()
# input file # input file
...@@ -224,7 +225,13 @@ if __name__ == "__main__": ...@@ -224,7 +225,13 @@ if __name__ == "__main__":
if key=="file": if key=="file":
value=os.path.basename(value) value=os.path.basename(value)
datainfo[key] = [descr, value] datainfo[key] = [descr, value]
datainfo['obsid'] = ['Observation Id', datainfo['file'][1].split("_")[0]] if key=="coord":
coord=datainfo[key][1]
(ra, dec, sign) = coord.split("-") + ["-"] if "-" in coord else coord.split("+") + ["+"]
dec=sign+dec
if key=="dmc" or key=="rmc" or key=="polc":
datainfo[key][1] = "True" if datainfo[key][1]==1 else "False"
datainfo['obsid'] = ['Observation Id', datainfo['file'][1].split("_")[0][1:]]
datainfo['pipeid'] = ['Pipeline Id', args.pipeid] datainfo['pipeid'] = ['Pipeline Id', args.pipeid]
# getting topo folding period # getting topo folding period
cmd="vap -n -c period %s | awk '{print $2}' -" % (infile) cmd="vap -n -c period %s | awk '{print $2}' -" % (infile)
...@@ -237,9 +244,24 @@ if __name__ == "__main__": ...@@ -237,9 +244,24 @@ if __name__ == "__main__":
tstart= str(os.popen(cmd).readlines()[0][:-1]) tstart= str(os.popen(cmd).readlines()[0][:-1])
ttime = tstart.split(" ")[-1] ttime = tstart.split(" ")[-1]
(day, month, year) = tstart.split(" ")[0].split(".") (day, month, year) = tstart.split(" ")[0].split(".")
datainfo['tstart'] = ['Start date/time: %s-%s-%sT%s' % (year, month, day, ttime), mjd] startobs="%s-%s-%sT%s" % (year, month, day, ttime)
datainfo['tstart'] = ['UTC start date/time: %s' % (startobs), mjd]
endobs=(datetime.strptime(startobs, '%Y-%m-%dT%H:%M:%S.%f')+timedelta(seconds=float(datainfo["length"][1]))).strftime('%Y-%m-%dT%H:%M:%S.%f')
midobs=(datetime.strptime(startobs, '%Y-%m-%dT%H:%M:%S.%f')+timedelta(seconds=float(datainfo["length"][1])/2.)).strftime('%Y-%m-%dT%H:%M:%S.%f')
timangles = { "startobs": [startobs], "midobs": [midobs], "endobs": [endobs] }
for key in timangles.keys():
cmd="azza.pl -ra %s -dec %s -t %s | grep EL | cut -d ' ' -f 3" % (ra, dec, timangles[key][0])
elev= str(os.popen(cmd).readlines()[0][:-1])
cmd="azza.pl -ra %s -dec %s -t %s | grep AZ | cut -d ' ' -f 3" % (ra, dec, timangles[key][0])
azim= str(os.popen(cmd).readlines()[0][:-1])
cmd="azza.pl -ra %s -dec %s -t %s | grep HA | cut -d ' ' -f 3-4" % (ra, dec, timangles[key][0])
ha= str(os.popen(cmd).readlines()[0][:-1])
timangles[key].extend([elev, azim, ha])
datainfo['elevation'] = ['Elevation (in deg) at the start / mid-point / end', "%s / %s / %s" % (timangles["startobs"][1], timangles["midobs"][1], timangles["endobs"][1])]
datainfo['azimuth'] = ['Azimuth (in deg) at the start / mid-point / end', "%s / %s / %s" % (timangles["startobs"][2], timangles["midobs"][2], timangles["endobs"][2])]
datainfo['transit'] = ['Hour angle at the start / mid-point / end', "%s / %s / %s" % (timangles["startobs"][3], timangles["midobs"][3], timangles["endobs"][3])]
if args.paz_rfi_fraction != "undef" or args.clfd_rfi_fraction != "undef": if args.paz_rfi_fraction != "undef" or args.clfd_rfi_fraction != "undef":
datainfo['rfi'] = ['RFI fraction in % by paz / clfd', "%s / %s" % (args.paz_rfi_fraction, args.clfd_rfi_fraction)] datainfo['rfi'] = ['RFI fraction (in %) by paz / clfd', "%s / %s" % (args.paz_rfi_fraction, args.clfd_rfi_fraction)]
# reading 1st file # reading 1st file
workraw = pc.Archive_load(infile) workraw = pc.Archive_load(infile)
...@@ -250,20 +272,20 @@ if __name__ == "__main__": ...@@ -250,20 +272,20 @@ if __name__ == "__main__":
bandpass=workraw.clone() bandpass=workraw.clone()
workraw.remove_baseline() workraw.remove_baseline()
workraw.centre_max_bin() workraw.centre_max_bin()
nchan1 = workraw.get_nchan() nchan = workraw.get_nchan()
nsubint1 = workraw.get_nsubint() nsubint = workraw.get_nsubint()
cfreq = workraw.get_centre_frequency() # center freq in MHz cfreq = workraw.get_centre_frequency() # center freq in MHz
bw = abs(workraw.get_bandwidth()) # bandwidth in MHz bw = abs(workraw.get_bandwidth()) # bandwidth in MHz
lowfreq = cfreq - bw/2.0 # lowest freq in MHz lowfreq = cfreq - bw/2.0 # lowest freq in MHz
chan_bw = bw/nchan1 # channel width in MHz chan_bw = bw/nchan # channel width in MHz
tobs = workraw.integration_length() # obs duration (in seconds) tobs = workraw.integration_length() # obs duration (in seconds)
stokesraw=workraw.clone() stokesraw=workraw.clone()
if nchan1 > 1: stokesraw.fscrunch() if nchan > 1: stokesraw.fscrunch()
if nsubint1 > 1: stokesraw.tscrunch() if nsubint > 1: stokesraw.tscrunch()
stokesraw.bscrunch(args.bscrunch) stokesraw.bscrunch(args.bscrunch)
nbins1 = stokesraw.get_nbin() nbins = stokesraw.get_nbin()
stokesraw.convert_state('Stokes') stokesraw.convert_state('Stokes')
r = stokesraw.get_data() r = stokesraw.get_data()
...@@ -283,14 +305,14 @@ if __name__ == "__main__": ...@@ -283,14 +305,14 @@ if __name__ == "__main__":
dynspecon=workraw.clone() dynspecon=workraw.clone()
dynspecon.pscrunch() dynspecon.pscrunch()
dynspecon.bscrunch(args.bscrunch) dynspecon.bscrunch(args.bscrunch)
nbins1 = dynspecon.get_nbin() nbins = dynspecon.get_nbin()
r = dynspecon.get_data() r = dynspecon.get_data()
#time stokes f phase #time stokes f phase
data5 = r[:,0,:,:] data5 = r[:,0,:,:]
dynspecon.tscrunch() dynspecon.tscrunch()
dynspecon.fscrunch() dynspecon.fscrunch()
prof = dynspecon.get_data()[0,0,0,:] prof = dynspecon.get_data()[0,0,0,:]
(prof, rot_bins, off_left, off_right) = auto_find_off_window(prof, 0, nbins1) (prof, rot_bins, off_left, off_right) = auto_find_off_window(prof, 0, nbins)
range_mean = np.mean(prof[off_left:off_right]) range_mean = np.mean(prof[off_left:off_right])
range_rms = np.std(prof[off_left:off_right]) range_rms = np.std(prof[off_left:off_right])
...@@ -300,13 +322,13 @@ if __name__ == "__main__": ...@@ -300,13 +322,13 @@ if __name__ == "__main__":
range_snrmean = np.mean(range_prof) range_snrmean = np.mean(range_prof)
range_weq = abs(np.sum(range_prof)/range_snrpeak) range_weq = abs(np.sum(range_prof)/range_snrpeak)
binshift = int(5*range_weq) binshift = int(5*range_weq)
if binshift > nbins1//2: binshift=nbins1//2 if binshift > nbins//2: binshift=nbins//2
lbin=nbins1//2-binshift lbin=nbins//2-binshift
rbin=nbins1//2+binshift rbin=nbins//2+binshift
phases=np.array([bin/nbins1 for bin in range(lbin,rbin)]) phases=np.array([bin/nbins for bin in range(lbin,rbin)])
range_crit=(range_prof>4.5) range_crit=(range_prof>4.5)
range_bins=np.arange(0,nbins1)[range_crit] range_bins=np.arange(0,nbins)[range_crit]
range_bins=trim_bins(range_bins) # trimming bins range_bins=trim_bins(range_bins) # trimming bins
range_profcrit=range_prof[range_bins] range_profcrit=range_prof[range_bins]
range_oncount=np.size(range_bins) range_oncount=np.size(range_bins)
...@@ -316,7 +338,7 @@ if __name__ == "__main__": ...@@ -316,7 +338,7 @@ if __name__ == "__main__":
psrstat_snr= float(os.popen(cmd).readlines()[0][:-1]) psrstat_snr= float(os.popen(cmd).readlines()[0][:-1])
datainfo['snr'] = ['Signal-to-noise ratio from psrstat / snr.py', "%.2f / %.2f" % (psrstat_snr, range_snrpeak)] datainfo['snr'] = ['Signal-to-noise ratio from psrstat / snr.py', "%.2f / %.2f" % (psrstat_snr, range_snrpeak)]
sp2=normalize_dynspec(data5, nchan1, nsubint1) sp2=normalize_dynspec(data5, nchan, nsubint)
dataon=np.mean(sp2[:,:,range_bins[0]:range_bins[-1]], axis=2) dataon=np.mean(sp2[:,:,range_bins[0]:range_bins[-1]], axis=2)
dataoff=np.mean(sp2[:,:,off_left:off_right], axis=2) dataoff=np.mean(sp2[:,:,off_left:off_right], axis=2)
v=dataon v=dataon
...@@ -340,7 +362,7 @@ if __name__ == "__main__": ...@@ -340,7 +362,7 @@ if __name__ == "__main__":
ax11 = fig.add_subplot(gs0[0, 0]) # profile ax11 = fig.add_subplot(gs0[0, 0]) # profile
ax11.plot(phases, (data1[0,lbin:rbin]-mean)/rms, '-', color='black', lw=1, label='I') ax11.plot(phases, (data1[0,lbin:rbin]-mean)/rms, '-', color='black', lw=1, label='I')
linpol=np.array([np.sqrt(data1[1][bin]*data1[1][bin]+data1[2][bin]*data1[2][bin]) for bin in range(nbins1)]) linpol=np.array([np.sqrt(data1[1][bin]*data1[1][bin]+data1[2][bin]*data1[2][bin]) for bin in range(nbins)])
ax11.plot(phases, (linpol[lbin:rbin]-mean)/rms, '-', color='red', lw=1, label='L') ax11.plot(phases, (linpol[lbin:rbin]-mean)/rms, '-', color='red', lw=1, label='L')
ax11.plot(phases, (data1[3,lbin:rbin]-mean)/rms, '-', color='blue', lw=1, label='V') ax11.plot(phases, (data1[3,lbin:rbin]-mean)/rms, '-', color='blue', lw=1, label='V')
ax11.legend() ax11.legend()
...@@ -352,15 +374,15 @@ if __name__ == "__main__": ...@@ -352,15 +374,15 @@ if __name__ == "__main__":
ax21 = fig.add_subplot(gs1[0, 0]) # stokes I ax21 = fig.add_subplot(gs1[0, 0]) # stokes I
stokes_iquv=workraw.clone() stokes_iquv=workraw.clone()
if nsubint1 > 1: stokes_iquv.tscrunch() if nsubint > 1: stokes_iquv.tscrunch()
stokes_iquv.bscrunch(args.bscrunch) stokes_iquv.bscrunch(args.bscrunch)
nbins1 = stokes_iquv.get_nbin() nbins = stokes_iquv.get_nbin()
stokes_iquv.convert_state('Stokes') stokes_iquv.convert_state('Stokes')
r = stokes_iquv.get_data() r = stokes_iquv.get_data()
#time stokes f phase #time stokes f phase
data2 = r[0,:,:,:] data2 = r[0,:,:,:]
sp=normalize(data2[0,:,:], nchan1) sp=normalize(data2[0,:,:], nchan)
v=sp[:,lbin:rbin] v=sp[:,lbin:rbin]
gmin, gmax = np.nanpercentile(v[(v!=-np.inf) & (v!=np.nan)], (0.05, 99.95)) gmin, gmax = np.nanpercentile(v[(v!=-np.inf) & (v!=np.nan)], (0.05, 99.95))
gmax = np.max([np.abs(gmin), np.abs(gmax)]) gmax = np.max([np.abs(gmin), np.abs(gmax)])
...@@ -370,7 +392,7 @@ if __name__ == "__main__": ...@@ -370,7 +392,7 @@ if __name__ == "__main__":
plt.ylabel("Frequency (MHz)") plt.ylabel("Frequency (MHz)")
ax22 = fig.add_subplot(gs1[0, 1]) # stokes Q ax22 = fig.add_subplot(gs1[0, 1]) # stokes Q
sp=normalize(data2[1,:,:], nchan1) sp=normalize(data2[1,:,:], nchan)
v=sp[:,lbin:rbin] v=sp[:,lbin:rbin]
gmin, gmax = np.nanpercentile(v[(v!=-np.inf) & (v!=np.nan)], (0.05, 99.95)) gmin, gmax = np.nanpercentile(v[(v!=-np.inf) & (v!=np.nan)], (0.05, 99.95))
gmax = np.max([np.abs(gmin), np.abs(gmax)]) gmax = np.max([np.abs(gmin), np.abs(gmax)])
...@@ -384,7 +406,7 @@ if __name__ == "__main__": ...@@ -384,7 +406,7 @@ if __name__ == "__main__":
plt.xlabel("Pulse phase", labelpad=0) plt.xlabel("Pulse phase", labelpad=0)
ax23 = fig.add_subplot(gs1[0, 2]) # stokes U ax23 = fig.add_subplot(gs1[0, 2]) # stokes U
sp=normalize(data2[2,:,:], nchan1) sp=normalize(data2[2,:,:], nchan)
v=sp[:,lbin:rbin] v=sp[:,lbin:rbin]
gmin, gmax = np.nanpercentile(v[(v!=-np.inf) & (v!=np.nan)], (0.05, 99.95)) gmin, gmax = np.nanpercentile(v[(v!=-np.inf) & (v!=np.nan)], (0.05, 99.95))
gmax = np.max([np.abs(gmin), np.abs(gmax)]) gmax = np.max([np.abs(gmin), np.abs(gmax)])
...@@ -398,7 +420,7 @@ if __name__ == "__main__": ...@@ -398,7 +420,7 @@ if __name__ == "__main__":
plt.xlabel("Pulse phase", labelpad=0) plt.xlabel("Pulse phase", labelpad=0)
ax24 = fig.add_subplot(gs1[0, 3]) # stokes V ax24 = fig.add_subplot(gs1[0, 3]) # stokes V
sp=normalize(data2[3,:,:], nchan1) sp=normalize(data2[3,:,:], nchan)
v=sp[:,lbin:rbin] v=sp[:,lbin:rbin]
gmin, gmax = np.nanpercentile(v[(v!=-np.inf) & (v!=np.nan)], (0.05, 99.95)) gmin, gmax = np.nanpercentile(v[(v!=-np.inf) & (v!=np.nan)], (0.05, 99.95))
gmax = np.max([np.abs(gmin), np.abs(gmax)]) gmax = np.max([np.abs(gmin), np.abs(gmax)])
...@@ -413,9 +435,9 @@ if __name__ == "__main__": ...@@ -413,9 +435,9 @@ if __name__ == "__main__":
ax31 = fig.add_subplot(gs2[0, 0]) # time-phase ax31 = fig.add_subplot(gs2[0, 0]) # time-phase
timephase=workraw.clone() timephase=workraw.clone()
if nchan1 > 1: timephase.fscrunch() if nchan > 1: timephase.fscrunch()
timephase.bscrunch(args.bscrunch) timephase.bscrunch(args.bscrunch)
nbins1 = timephase.get_nbin() nbins = timephase.get_nbin()
timephase.pscrunch() timephase.pscrunch()
r = timephase.get_data() r = timephase.get_data()
#time stokes f phase #time stokes f phase
...@@ -432,7 +454,7 @@ if __name__ == "__main__": ...@@ -432,7 +454,7 @@ if __name__ == "__main__":
#time stokes f phase #time stokes f phase
data4 = np.mean(np.sum(r[:,:,:,:], axis=3), axis=0) data4 = np.mean(np.sum(r[:,:,:,:], axis=3), axis=0)
dB = lambda x: 10 * np.log10(x) dB = lambda x: 10 * np.log10(x)
freqs=np.array([lowfreq + chan * chan_bw for chan in range(nchan1)]) freqs=np.array([lowfreq + chan * chan_bw for chan in range(nchan)])
ax32.plot(freqs, dB(data4[0,:]), color="k", label="PP", lw=0.5) ax32.plot(freqs, dB(data4[0,:]), color="k", label="PP", lw=0.5)
ax32.plot(freqs, dB(data4[1,:]), color="b", label="QQ", lw=0.5) ax32.plot(freqs, dB(data4[1,:]), color="b", label="QQ", lw=0.5)
#ax32.plot(freqs, dB(data4[2,:]), color="g", label="Re[PQ]") #ax32.plot(freqs, dB(data4[2,:]), color="g", label="Re[PQ]")
...@@ -452,7 +474,8 @@ if __name__ == "__main__": ...@@ -452,7 +474,8 @@ if __name__ == "__main__":
ax12.yaxis.set_major_formatter(FormatStrFormatter("")) ax12.yaxis.set_major_formatter(FormatStrFormatter(""))
ax12.yaxis.set_minor_formatter(FormatStrFormatter("")) ax12.yaxis.set_minor_formatter(FormatStrFormatter(""))
ordered_fields_to_show=['obsid', 'pipeid', 'file', 'name', 'coord', 'tstart', 'length', 'nbin', 'nchan', 'npol', 'nsubint', 'period', 'dm', 'rm', 'freq', 'bw', 'site', 'rcvr:name', 'be:name', 'snr', 'rfi'] ordered_fields_to_show=['obsid', 'pipeid', 'file', 'name', 'coord', 'tstart', 'length', 'nbin', 'nchan', 'npol', 'nsubint', \
'period', 'dm', 'rm', 'freq', 'bw', 'dmc', 'rmc', 'polc', 'site', 'rcvr:name', 'be:name', 'snr', 'rfi', 'elevation', 'azimuth', 'transit']
key_str=val_str=descr_str="" key_str=val_str=descr_str=""
for key in ordered_fields_to_show: for key in ordered_fields_to_show:
if key not in datainfo: continue if key not in datainfo: continue
...@@ -483,6 +506,9 @@ if __name__ == "__main__": ...@@ -483,6 +506,9 @@ if __name__ == "__main__":
gs.tight_layout(fig) gs.tight_layout(fig)
plt.tight_layout() plt.tight_layout()
plt.savefig(infile.split(".ar")[0] + ".quickview.png", bbox_inches="tight") if args.outstem == "":
args.outstem = infile.split(".ar")[0]
args.outstem += ".quickview.png"
plt.savefig(args.outstem, bbox_inches="tight")
if not args.saveonly: if not args.saveonly:
plt.show() plt.show()
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment