<jupyter_start><jupyter_text>AR1 Observation report  File to be processed<jupyter_code>filename =<jupyter_output><empty_output><jupyter_text>Set up environment<jupyter_code>%pylab inline
import katdal
import scape
import pickle
import katpoint
import pysolr
import h5py
import warnings
warnings.simplefilter('ignore')
scape._logger.setLevel('ERROR')
def get_reduction_metadata(filename, reduction_name=None):
    #Get all reduction products from filename  with given reduction_name 
    #(or all reduction products if reduction_name is None)
    mysolr = pysolr.Solr('http://kat-archive.kat.ac.za:8983/solr/kat_core')
    fn_search_result = mysolr.search('Filename:'+filename)
    if fn_search_result.hits < 1:
        return []
    CASProductId = fn_search_result.docs[0]['CAS.ProductId']
    reduction_products = mysolr.search('InputDataProductId:'+CASProductId)
    if reduction_name==None:
        return reduction_products.docs
    else:
        return [product for product in reduction_products.docs if product.get('ReductionName')==reduction_name]
try:
    res = get_reduction_metadata(filename.split('/')[-1],reduction_name = 'AR1 Generate RFI Flags');
    flag_results = res[0]['CAS.ReferenceDatastore']
    for fl in flag_results:
        if fl.endswith('flags.h5'):
            flag_filename = fl.split(':')[-1]
    print flag_filename
    h5_flags = h5py.File(flag_filename,'r')
    h5_flags['flags'].shape
    cal_flags = True
except:
    print "no calibration pipeline flags found"
    cal_flags=False<jupyter_output><empty_output><jupyter_text>Load data<jupyter_code>data = katdal.open(filename)
if cal_flags:
    data._flags = h5_flags['flags']

N = data.shape[0]
ext = 930
step = max(int(floor(N/ext)),1)
print step
print data
data.obs_params<jupyter_output><empty_output><jupyter_text>Load flags<jupyter_code>data.select()
M = 4 * data.shape[1] / 4096

start_chan = 2200*M//4
end_chan = 2800*M//4
if data.receivers[data.receivers.keys()[0]][0] == 'u' : 
    start_chan = 0#2200*M//4
    end_chan = data.shape[1]#2800*M//4    
data.select(channels=slice(start_chan,end_chan))
flags = data.flags[:]

figure(figsize=(10,5))
pcolormesh(data.timestamps-data.timestamps[0],data.channel_freqs/1e6,sum(flags,axis=2).T)
title('Flags in "quiet" part of the band')
xlabel('Time (s), since %s' % (katpoint.Timestamp(data.start_time).local(),))
ylabel('Frequency/[MHz]')
xlim(0,data.timestamps[-1]-data.timestamps[0])
ylim(data.channel_freqs[0]/1e6,data.channel_freqs[-1]/1e6)
dat = data.file['Data/flags'][:,slice(start_chan,end_chan),:].flatten()
labels =  list(data.file['/Data/flags_description'][:,0])
plt.xticks(range(len(labels)), labels, rotation='vertical')
plt.ylabel('Percentage Flagged')
plt.title('Flag Types (%i samples)'%(dat.shape[0]))
plt.plot(np.unpackbits(dat[:,np.newaxis],axis=1).sum(axis=0)/np.float(dat.shape[0])*100, '*')
plt.grid()<jupyter_output><empty_output><jupyter_text>Plot timeseries, spectrograph and spectrum for each antenna<jupyter_code>activity = []
for ant in data.ants:
    activity.append(data.sensor['Antennas/%s/activity'%ant.name])
scan_types = set(hstack(activity))

activity_count = zeros([len(scan_types),len(data.ants)])
for s,st in enumerate(scan_types):
    for i in range(len(data.ants)):
        activity_count[s,i] = sum(vstack(activity)[i,:]==st)

figure(figsize=(12,2))
for s,st in enumerate(scan_types):
    plot(np.arange(len(data.ants)),activity_count[s,:],"o",label=st)
plt.xticks(np.arange(len(data.ants)), [ant.name for ant in data.ants])
legend(loc=5)
grid()
xlim(-1,len(data.ants)+2)
title('Number of dumps in each scan state')
mpl.rc("savefig", dpi=300)
for p,active_pol in enumerate(['H','V']):
    data.select()
    data.select(pol=active_pol,scans=['scan','track'],dumps=slice(0,-1,step),channels=slice(0,-1,M))
    num_ants = len(data.ants)
    baseline_names = [('%s - %s' % (inpA[:-1], inpB[:-1])) for inpA, inpB in data.corr_products[num_ants:]]
    
    
    num_bls = data.shape[2] - num_ants
    
    scan_targets, scan_az, scan_el, scan_timestamps, scan_vis = [], [], [], [], []
    for scan_ind, state, target in data.scans():
        num_ts = data.shape[0]
        ts = data.timestamps
        vis = data.vis[:]
        flags = data.flags[:]
        #vis[:,:,:num_ants][flags[:,:,:num_ants]] = np.nan
        vis[flags] = np.nan
        scan_vis.append(vis)
        scan_az.append(median(data.az,axis=0))
        scan_el.append(median(data.el,axis=0))
        scan_timestamps.append(ts - data.start_time.secs)
    num_chans, sample_period = len(data.channel_freqs), data.dump_period
    scan_freqinds = [np.arange(num_bls * num_chans)] * len(scan_timestamps)
    scan_phase = [np.angle(vis[:,:,num_ants:]).T.reshape(-1, vis.shape[0]) for vis in scan_vis]
    scan_autos = [np.abs(vis[:,:,:num_ants]) for vis in scan_vis]
    x1,x2 = scan_timestamps[0][0],scan_timestamps[-1][-1]
    y1,y2 = scan_freqinds[0][0]-0.5,scan_freqinds[0][-1]+0.5
    f1,f2 = data.channel_freqs[0]/1e6, data.channel_freqs[-1]/1e6
    for i in range(num_ants):
        ant_name = data.corr_products[i][0]
        amp = vstack([scan[:,:,i] for scan in scan_autos])
        az = hstack([a[i] for a in scan_az])
        el = hstack([e[i] for e in scan_el])
        fig, [[ax1, ax2], [ax3, ax4]] = subplots(2, 2, figsize=(12, 8))
        fig.suptitle('Ant: %s'%(ant_name), fontsize=16)
        fig.subplots_adjust(top=0.95)

        ax1.errorbar(hstack(scan_timestamps),nanmean(amp,axis=1),nanstd(amp,axis=1),color='grey')
        ax1.fill_between(hstack(scan_timestamps), np.nanmin(amp,axis=1), np.nanmax(amp,axis=1),color='lightgrey')
        ax1.plot(hstack(scan_timestamps),np.nanmean(amp,axis=1),color='b')
        ax1.set_title('Time series')
        ax1.set_xlabel('Time (s), since %s' % (katpoint.Timestamp(data.start_time).local(),))
        ax1.set_ylabel('Correlator counts')
        ax1.set_ylim(0,4000)
        ax1.set_xlim(x1,x2)
        for t in scan_timestamps:
            ax1.axvline(t[0],linestyle='dotted',color='mistyrose',linewidth=1)
        
        ma_amp = ma.masked_where(np.isnan(amp.T),amp.T)
        ax2.pcolormesh(hstack(scan_timestamps),data.channel_freqs/1e6,10*log10(ma_amp))
        ax2.set_title('Spectrogram')
        ax2.set_xlabel('Time (s), since %s' % (katpoint.Timestamp(data.start_time).local(),))
        ax2.set_xlim(x1,x2)
        ax2.set_ylabel('Frequency/[MHz]')
        ax2.set_ylim(f1,f2)
        
        ax3.errorbar(data.channel_freqs/1e6,nanmean(amp,axis=0),nanstd(amp,axis=0),color='grey')
        ax3.fill_between(data.channel_freqs/1e6, np.nanmin(amp,axis=0), np.nanmax(amp,axis=0),color='lightgrey')
        ax3.plot(data.channel_freqs/1e6,nanmean(amp,axis=0),color='b')
        ax3.set_title('Spectrum')
        ax3.set_xlabel('Frequency/[MHz]')
        ax3.set_xlim(f1,f2)
        ax3.set_ylabel('Correlator counts')
        ax3.set_ylim(0,4000)
        
        ax4.plot(az,el,"*")
        ax4.set_title('Pointing')
        ax4.set_xlabel('Azimuth/[Degrees]')
        ax4.set_ylabel('Elevation/[Degrees]')
        
    #An 'inch' per baseline
    if num_bls > 0 :
        figure(0,figsize=(12,num_bls),tight_layout=True)
    else :
        figure(0,figsize=(12,1),tight_layout=True)
    #Put plots side-by-side if there are <10 scans.     
    if len(scan_vis)<10:
        subplot (1,2,p+1)
    else:
        subplot(2,1,p+1)
    #segms, labels, lines = scape.plots_basic.plot_segments(scan_timestamps, scan_freqinds, scan_phase, labels=scan_targets,add_breaks=False)
    #imshow(hstack(scan_phase),aspect='auto',origin='lower',extent=[x1,x2,y1,y2],interpolation='none')
    sp=hstack(scan_phase)
    imshow(sp[~np.all(np.isnan(sp),axis=1)],aspect='auto',origin='lower',extent=[x1,x2,y1,y2],interpolation='nearest')
    plt.xlabel('Time (s), since %s' % (katpoint.Timestamp(data.start_time).local(),))
    if p == 1 and len(scan_vis) < 10:
        plt.yticks(np.arange(num_chans // 2, num_bls * num_chans, num_chans), repeat('',num_bls))
    else:
        plt.yticks(np.arange(num_chans // 2, num_bls * num_chans, num_chans), baseline_names)
    for yval in range(0, num_bls * num_chans, num_chans):
        plt.axhline(yval, color='k', lw=2)
    plt.title('Raw visibility phase per baseline, Pol %s'%active_pol)
data.timestamps[1] - data.timestamps[0]
for ant in data.ants:
    print ant.description
print "Scape version: %s"%scape.__version__
print "katdal version: %s"%katdal.__version__
print "katpoint version: %s"%katpoint.__version__
print 'Report Ended:'
!date<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup