Commit 5ae93f0f authored by perdereau's avatar perdereau
Browse files

python v3

parent 9d32df75
......@@ -2,7 +2,7 @@ import math as math
import sys, platform, os
import numpy as np
from astropy.io import fits
import matplotlib as matplotlib
from matplotlib import pyplot as plt
import itertools as itert
from skimage import data, img_as_float
......@@ -14,7 +14,7 @@ from colorsys import hsv_to_rgb, hls_to_rgb
# Choices: ['auto', 'gtk', 'gtk3', 'inline', 'nbagg', 'notebook', 'osx', 'qt', 'qt4', 'qt5', 'tk', 'wx']
# values of the freqs for the slices
rfreq = [1381.,1390.,1420.5,1430.]
rfreq = [1347.,1390.,1420.5,1430.]
normf=[1325.,1465.] # "normal frequencies"
......@@ -22,6 +22,12 @@ names = ['1H','2H','3H','4H','1V','2V','3V','4V','1Hx2H','1Hx3H','1Hx4H','2Hx3H'
'1Vx2V','1Vx3V','1Vx4V','2Vx3V','2Vx4V','3Vx4V','1Hx1V','1Hx2V','1Hx3V','1Hx4V','2Hx1V','2Hx2V','2Hx3V','2Hx4V',
'3Hx1V','3Hx2V','3Hx3V','3Hx4V','4Hx1V','4Hx2V','4Hx3V','4Hx4V']
font = {'family' : 'normal',
'weight' : 'normal',
'size' : 20}
matplotlib.rc('font', **font)
# viewtfmap reads and returns data stored in the TF maps fits files build by visiavg (code in AnaPAON4)
# all TF maps are stored as hdu in the fits file, with some additional data at the end.
......@@ -56,7 +62,8 @@ names = ['1H','2H','3H','4H','1V','2V','3V','4V','1Hx2H','1Hx3H','1Hx4H','2Hx3H'
# hasvar : to handle the cas where there are no variances
# ramod : to handle the case where the time-freq maps are ra-frequency maps [deprecated]
# mask : if true, also looks for a 'mask' file (same name+'_mask') to mask RFIs and satellites
# yrcm : to use another color map (
# yrcm : to use another color map
# usecm : choice of colormap
# raplot : also shows slices at some frequencies [see array rfreq above] as a fct of RA [deprecated - don't ew
# timeplot : also shows slices at soem frequencies as a fct of time [deprecated - don't use]
......@@ -76,11 +83,11 @@ names = ['1H','2H','3H','4H','1V','2V','3V','4V','1Hx2H','1Hx3H','1Hx4H','2Hx3H'
#
def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=False,mtitle="",clean=True,mod=False,
arg=False,cmplx=False,noplo=False,verb=False,yrcm=False,cor=[],cimag=False,creal=False,hasvar=True,
ramod=False,plras=[],nmras=[],mask=False,logsc=False,snratio=False):
arg=False,cmplx=False,noplo=False,verb=False,yrcm=False,cor=[],cimag=False,creal=False,hasvar=True,submed=False,
ramod=False,plras=[],nmras=[],mask=False,logsc=False,snratio=False,usecm="",extaxes=False,choose_cmap=False,nau=True) :
if verb :
print len(cor)
print (len(cor))
filename = folder
if noplo :
clean=False
......@@ -113,14 +120,14 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
i=num # devrait etre juste la pour les sorties des voies th et rfi
if (verb) :
print hdulist[i].header
print( hdulist[i].header)
# get the bin indexes for the slices & the "normal range"
bfreq = [ (np.abs(freqs-rf)).argmin() for rf in rfreq ]
bfnor = [(np.abs(freqs-rf)).argmin() for rf in normf ]
if (verb) :
print ' number of hdus ' + str(len(hdulist))
print (' number of hdus ' + str(len(hdulist)))
nfre = np.size(freqs)
......@@ -129,7 +136,7 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
if timin ==0. :
timin=times[0]/3600.
timax = max(times)/3600.
print ' timin/max = '+str(timin)+', ' +str(timax)
print (' timin/max = '+str(timin)+', ' +str(timax))
if ramod :
timin = min(ras)
timax= max(ras)
......@@ -148,13 +155,26 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
if(cimag):
img = hdulist[ i+1].data
if (verb) :
print 'Image shape :'
print np.shape(img)
# module
print ('Image shape :')
print (np.shape(img))
# module
meds=np.zeros(0)
if submed :
meds=np.median(img,axis=1)
print (len(meds))
print (np.shape(img))
for idx in np.arange(np.shape(img)[0]):
img[idx]=img[idx]-meds[idx]
if mod :
if num>7 :
imgi= hdulist[ i+1].data
if submed :
medsi = np.median(imgi,axis=1)
for idx in np.arange(np.shape(imgi)[0]):
imgi[idx]=imgi[idx]-medsi[idx]
img = np.sqrt(img**2+imgi**2)
if len(cor) ==2 : # for correlation coeff
imca = hdulist[cor[0]].data
......@@ -164,17 +184,25 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
if arg :
if num>7 :
imgi= hdulist[ i+1].data
if submed :
medsi = np.median(imgi,axis=1)
for idx in np.arange(np.shape(imgi)[0]):
imgi[idx]=imgi[idx]-medsi[idx]
img = np.arctan2(imgi,img) * 180. / np.pi
mvmin=-180.
mvmax=180.
mtitle= names[num]+" arg.(deg)"
else :
print "arg demande mais auto selectionne : ERREUR "
print ("arg demande mais auto selectionne : ERREUR ")
exit
mskdat=img*0.
# complex number representation
if cmplx :
imgi= hdulist[ i+1].data
if submed :
medsi = np.median(imgi,axis=1)
for idx in np.arange(np.shape(imgi)[0]):
imgi[idx]=imgi[idx]-medsi[idx]
zmod = np.sqrt(img**2+imgi**2)
arg = np.arctan2(imgi,img)
......@@ -182,9 +210,9 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
aa = (aa+0.5)%1
bb = 1./(abs(zmod)**0.3 +1.)
bb = -bb+1.
print np.shape(aa)
cc = np.ndarray(np.shape(aa))
print np.shape(cc)
cc.fill(1.)
iimg =[]
for ii in np.arange(np.shape(aa)[0]):
......@@ -197,15 +225,16 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
if mask : # reads mask file
if cmplx :
print "mask & cmplx not yet compatible "
print ("mask & cmplx not yet compatible ")
exit
msknam = filename +'_masks.fits'
print msknam+ " " + str(num)
hdmskl = fits.open(msknam)
mskdat = np.fabs(hdmskl[num+1].data-1.)
img = img*mskdat
filename=filename+"_masked"
if logsc :
filename=filename+" (log)"
if mod :
......@@ -224,48 +253,55 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
if mvmax !=0. :
vmax = mvmax
if verb :
print vmin, np.asarray(img[300:900,]).min()
print vmax,np.asarray(img[300:900,]).max()
print (vmin, np.asarray(img[300:900,]).min())
print (vmax,np.asarray(img[300:900,]).max())
if (not noplo) :
# plots TF image - dual axes with bin no and values
im = ax1.imshow(img,aspect='auto' , interpolation='none',vmin=vmin,vmax=vmax)
im = ax1.imshow(img,aspect='auto' ,extent=(timin,timax,frmax,frmin), interpolation='none',vmin=vmin,vmax=vmax)
#m = ax1.imshow(img,aspect='auto' ,extent=(timin,timax,frmax,frmin), interpolation='none',vmin=vmin,vmax=vmax)
if (yrcm) :
im.set_cmap('YlOrBr_r')
im.set_cmap(usecm)
cbar = plt.colorbar( im ) #,fraction=0.021,pad=.05)
cbar.set_norm(mynorm.MyNormalize(vmin=vmin,vmax=vmax,stretch='linear'))
cm_cycle = sorted([i for i in dir(plt.cm) if hasattr(getattr(plt.cm,i),'N')])
cm_index = cm_cycle.index(cbar.get_cmap().name)
if nau :
cbar.ax.set_ylabel('signal (NAU) ', rotation=270,labelpad=20)
else:
cbar.ax.set_ylabel('signal ', rotation=270,labelpad=20)
plt.show()
if mtitle =="":
plt.title(filename + ' ' + names[num],y=1.08)
else :
plt.title(filename + ' ' +mtitle,y=1.08)
plt.title(mtitle+ " - " + names[num] ,y=1.08)
plt.ylabel("Frequency (Mhz)")
if ramod:
plt.xlabel("RA (h)")
else:
plt.xlabel("Temps (h)")
plt.xlabel("Time (h)")
if ramod :
for rap in plras:
ax1.plot([rap,rap],[frmax,frmin],color='red',linestyle='--')
#ax1.annotate(nmras[plras.index(rap)],(rap,frmin),horizontalalignment='center',verticalalignment='bottom')
ax1.annotate(nmras[plras.index(rap)],(rap,frmin),horizontalalignment='left',verticalalignment='bottom',rotation=30)
ax1.yaxis.grid(which="major", color='gray', linestyle='-', linewidth=.5)
ax1.xaxis.grid(which="major", color='gray', linestyle='-', linewidth=.5)
ax2=ax1.twinx()
ax2.yaxis.set_view_interval(frmax,frmin ,ignore=True)
ax2.set_ylabel("Freq. bin",rotation=270.,y=0.9)
if not ramod :
ax3=ax1.twiny()
ax3.xaxis.set_view_interval( timin,timax ,ignore=True)
ax3.set_xlabel("time bin",x=0.85)
if extaxes :
ax2=ax1.twinx()
ax2.yaxis.set_view_interval(frmax,frmin ,ignore=True)
#ax2.set_ylabel("Freq. bin",rotation=270.,y=0.9)
if not ramod :
ax3=ax1.twiny()
ax3.xaxis.set_view_interval( timin,timax ,ignore=True)
#ax3.set_xlabel("time bin",x=0.85)
plt.tight_layout()
if save: # saves png with some name
ext=""
if creal :
......@@ -280,13 +316,8 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
else :
plt.savefig(filename+"_"+ names[num]+"_TFM"+ext+".png")
# ------------ THIS IS TO SWITCH FROM ONE CMAP TO THE OTHER
axnext = fig.add_axes([0.6, 0.04, 0.05, 0.05])
bnext = Button(axnext,cm_cycle[cm_index+1] )
axprev = fig.add_axes([0.7, 0.04, 0.05, 0.05])
bprev = Button(axprev,cm_cycle[cm_index-1] )
axcur = fig.add_axes([0.65, 0.04, 0.05, 0.05], )
bcur = Button(axcur,cm_cycle[cm_index])
# ------------ THIS IS TO SWITCH FROM ONE CMAP TO THE OTHER
def cm_next(val):
cm_cycle = sorted([i for i in dir(plt.cm) if hasattr(getattr(plt.cm,i),'N')])
cm_index = cm_cycle.index(cbar.get_cmap().name)
......@@ -327,9 +358,16 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
indnext=0
bnext.label.set_text(cm_cycle[indnext])
bprev.label.set_text(cm_cycle[indprev])
if choose_cmap :
axnext = fig.add_axes([0.6, 0.04, 0.05, 0.05])
bnext = Button(axnext,cm_cycle[cm_index+1] )
axprev = fig.add_axes([0.7, 0.04, 0.05, 0.05])
bprev = Button(axprev,cm_cycle[cm_index-1] )
axcur = fig.add_axes([0.65, 0.04, 0.05, 0.05], )
bcur = Button(axcur,cm_cycle[cm_index])
bnext.on_clicked(cm_next)
bprev.on_clicked(cm_prev)
bnext.on_clicked(cm_next)
bprev.on_clicked(cm_prev)
#----------------- END CMAP SWITCH
#plt.tight_layout()
......@@ -373,20 +411,20 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
legobj.set_linewidth(2.0)
plt.ylim((mvmin,mvmax))
if mtitle =="":
plt.title(filename + ' ' + names[num],y=1.08)
else :
plt.title(filename + ' ' +mtitle,y=1.08)
plt.title(mtitle+ " - " + names[num] ,y=1.08)
plt.tight_layout()
if save:
if mask :
plt.savefig(filename+"_"+ names[num]+"_TFM_mask_freq.png")
else:
plt.savefig(filename+"_"+ names[num]+"_TFM_freq.png")
if (verb) :
print times.size
if (verb) :
print np.shape(img)
print (times.size())
print (np.shape(img))
#return values
......@@ -395,7 +433,7 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
elif cmplx:
return times/3600.,img,freqs,ras
else :
return times/3600.,img,freqs,ras,[ img[b,] for b in bfreq],ax1,bnext,bprev,mskdat
return times/3600.,img,freqs,ras,[ img[b,] for b in bfreq],ax1,mskdat
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment