In [12]:
import glob, os, sys
from astropy.io import fits as pyfits
import matplotlib.pyplot as plt
from astropy.io import ascii
import numpy as np
In [43]:
#The new data location
if sys.platform == 'darwin': BASEDIR='/Volumes/PLUME/MWA_DATA/'
if sys.platform == 'linux2': BASEDIR='/mnt/MWA_DATA/'
CHANNELS= ['113-114','125-126','139-140','153-154','169-170','187-188']#'062-063',
#'139-140'#'125-126'#'113-114' #'069-070'#'093-094' #'084-085'#'076-077'#'062-063'
#CHANNELS= ['153-154']
polarization='XX'
INDIR=BASEDIR
date='2015/11/04 '
OBSIDS=['1130642640','1130642936','1130643240','1130643536','1130643840','1130644136','1130644440','1130644736','1130645040']
#OBSIDS=['1130642936']
FIRSTOBSID=OBSIDS[0]
varydr=1 #1 to vary the dynamic range, 0 to keep it fixed
force=1 #Overwrite files if present
oplot_maxima=0
bdiff=0
rdiff=0
#Set this to 1 to create images integrated over time
summed=1
nsummed=60
dpi=40 #dots per inch for the plotting
finchan=[8,14] #the fine channel indices
#The [y0,y1,x0,x1] coordinates of the subplot for each channel
peak_box={'062-063':{'x':[280,629],'y':[300,649]},
'069-070':{'x':[280,629],'y':[300,649]},
'076-077':{'x':[240,589],'y':[350,699]},
'084-085':{'x':[280,629],'y':[300,649]},
'093-094':{'x':[230,579],'y':[340,689]},
'103-104':{'x':[320,669],'y':[290,639]},
'113-114':{'x':[230,579],'y':[340,689]},
'125-126':{'x':[230,579],'y':[340,689]},
'139-140':{'x':[230,579],'y':[340,689]},
'153-154':{'x':[290,639],'y':[300,649]},
'169-170':{'x':[250,599],'y':[320,669]},
'187-188':{'x':[315,664],'y':[285,634]}}
In [45]:
for CHANNEL in CHANNELS:
if bdiff > 0: base_loaded=0
#GET the Peak and RMS for the Dynamic Range of the image.
if varydr > 0:
drfile='DR_'+CHANNEL+'_'+polarization+'.txt'
datadir=BASEDIR+CHANNEL+'/'
print datadir+drfile
if os.path.exists(datadir+drfile):
drdata=ascii.read(datadir+drfile)
tmp=drdata['col1']
fnames=[]
dr_timestrings=[]
for fname in tmp:
#fnames.append(os.path.basename(fname.split('.image')[0]))
filebasename=os.path.basename(fname.split('.image')[0])
obsid=filebasename.split('_')[0]
fnames.append(filebasename+'.fits')
timebase=filebasename.split('_t')[1].split('_')[0]
time=timebase[0:2]+':'+timebase[2:4]+':'+timebase[4:6]
dr_timestrings.append(timebase)
peak=drdata['col2']
rms=drdata['col3']
dr=drdata['col4']
else:
varydr=0
#print fnames
#CHECK if Maxima files exist, record the data
if oplot_maxima > 0:
maxindices=['1','2']
tmpmaxinfo={}
maxinfo={}
for maxindex in maxindices:
maxfile='Max'+maxindex+'_info_'+CHANNEL+'_'+polarization+'.txt'
if os.path.exists(datadir+maxfile):
maxdata=ascii.read(datadir+maxfile)
max_timestrings=[]
for time in maxdata['times']:
max_timestrings.append(''.join(time.split(' ')[1].split(':')))
maxdata['timestring']=max_timestrings
tmpmaxinfo[maxindex]=maxdata
maxinfo[CHANNEL]=tmpmaxinfo
#Calculate the frequencies
tmp=CHANNEL.split('-')
basefreq=int(tmp[0])*1.28 #base frequency in MHz
startfreq=basefreq+finchan[0]*0.04 #Starting frequency
endfreq=basefreq+finchan[1]*0.04 #Starting frequency
#The sub-FOV
sub=[peak_box[CHANNEL]['y'][0],peak_box[CHANNEL]['y'][1],peak_box[CHANNEL]['x'][0],peak_box[CHANNEL]['x'][1]]
plt.ioff() # Turn interactive plotting off: NOT to show each frame
for OBSID in OBSIDS:
INDIR=BASEDIR+CHANNEL+'/'+OBSID+'/'
print 'Working in '+INDIR
img_list = glob.glob(INDIR+'*'+polarization+'*_d001.fits')
img_list.sort()
if len(img_list)==0:
print "#### Found No image matching the search criterion ####";
if bdiff > 0 and OBSID == FIRSTOBSID and base_loaded == 0:
BASEIMG=img_list[0]
hd=pyfits.open(BASEIMG)
basedata=hd[0].data[0,0,sub[0]:sub[1],sub[2]:sub[3]]
hd.close()
base_loaded=1
sumc=0
for IMG in img_list:
NAME = IMG.split('.fits')[0]
OUT_NAME = NAME+'.png'
if bdiff:
OUT_NAME = NAME+'_bdiff.png'
elif rdiff:
OUT_NAME = NAME+'_rdiff.png'
elif summed:
OUT_NAME = NAME+'_summed'+str(summed)+'s.png'
if force > 0:
os.system('rm -f '+OUT_NAME)
#dynrange=[-2,300]
#print os.path.basename(NAME+'.fits')
#print fnames
if varydr > 0 and os.path.basename(NAME+'.fits') in fnames:
imgind=fnames.index(os.path.basename(NAME+'.fits'))
dynrange=[rms[imgind]*0.9,1.0*peak[imgind]]
if not os.path.exists(OUT_NAME):
timestring=IMG.split('_t')[1].split('_')[0]
time=timestring[0:2]+':'+timestring[2:4]+':'+timestring[4:6]
hd=pyfits.open(IMG)
#print IMG
#print hd[0].data.shape
#print imgind
origdata=hd[0].data[0,0,sub[0]:sub[1],sub[2]:sub[3]]
if bdiff:
#Run this to create base difference images
#SUBTRACT THE TWO IMAGES
plotdata=np.subtract(origdata,basedata)
#np.min(data)*0.9,1.0*np.max(data)
dynrange=[0.,1.0*np.max(plotdata)]
elif rdiff and imgind > 0:
#Run this to create running difference images
#SUBTRACT THE TWO IMAGES
plotdata=np.subtract(origdata,rundata)
#np.min(data)*0.9,1.0*np.max(data)
dynrange=[0.,1.0*np.max(plotdata)]
elif summed:
#Run this to get time-integrated images
if sumc > 1 and sumc % nsummed == 0:
plotdata = sumdata
dynrange=[0.,1.0*np.max(plotdata)]
sumdata = origdata
sumc=sumc+1
else:
if sumc == 0: sumdata = origdata
else: sumdata += origdata
sumc=sumc+1
continue
else:
plotdata=origdata
dynrange=[0.,1.0*np.max(plotdata)]
ypixels,xpixels=plotdata.shape
fig=plt.figure(figsize=(xpixels/dpi,ypixels/dpi))
ax=plt.axes([0.08, 0.08, 0.87, 0.87])#frameon=False, xticks=[],yticks=[]
ax.imshow(plotdata,vmin=dynrange[0],vmax=dynrange[1], origin='lower', interpolation='none')
print imgind
printstr = 'MWA '+date + time +' | C' + CHANNEL +' ({:.2f}-{:.2f} MHz)'.format(startfreq,endfreq)
printstr+=' | DR={:.2f}'.format(dr[imgind])
#printstr+=' | DR={:.2f}-{:.2f} Jy'.format(dynrange[0],dynrange[1])
plt.text(xpixels*0.02, ypixels*0.97, printstr,ha='left',va='center',color='white',fontsize=11)
#drstr='DR:['+str(dynrange[0].format('%.2f'))+'-'+str(dynrange[1].format('%.2f'))+']'
#PLOT THE MAX POSITIONS
if oplot_maxima > 0:
ax.autoscale(False)
maxindices=['1','2']
for maxindex in maxindices:
if timestring in maxinfo[CHANNEL][maxindex]['timestring']:
imgind=maxinfo[CHANNEL][maxindex]['timestring'].tolist().index(timestring)
if maxindex == '1':
plotsym='go'
if maxindex == '2':
plotsym='ro'
xpos=maxinfo[CHANNEL][maxindex]['maxlocx_px'][imgind] - peak_box[CHANNEL]['x'][0]
ypos=maxinfo[CHANNEL][maxindex]['maxlocy_px'][imgind] - peak_box[CHANNEL]['y'][0]
ax.plot(xpos,ypos,plotsym,markersize=7)
#plt.text(xpixels*0.9,ypixels*0.96, drstr,
# ha='center',va='center',color='white')
print OUT_NAME
plt.savefig(OUT_NAME)
plt.close()
hd.close()
if rdiff: rundata = origdata
if bdiff > 0: base_loaded=0
plt.ion() # Turn interactive plotting back on
In [30]:
print imgind
In [42]:
322 % 30
Out[42]:
In [ ]: