📣 An issue occured with the embedded container registry on October 25 2021, between 10:30 and 12:10 (UTC+2). Any persisting issues should be reported to CC-IN2P3 Support. 🐛

Commit 4caca2db authored by perdereau's avatar perdereau
Browse files

add comments (viewtfmap) + set pressure to 0 (altaz.py)

parent b96d0cb6
......@@ -13,20 +13,75 @@ import time
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.]
normf=[1325.,1465.]
normf=[1325.,1465.] # "normal frequencies"
names = ['1H','2H','3H','4H','1V','2V','3V','4V','1Hx2H','1Hx3H','1Hx4H','2Hx3H','2Hx4H','3Hx4H',
'1Vx2V','1Vx3V','1Vx4V','2Vx3V','2Vx4V','3Vx4V','1Hx1V','1Hx2V','1Hx3V','1Hx4V','2Hx1V','2Hx2V','2Hx3V','2Hx4V',
'3Hx1V','3Hx2V','3Hx3V','3Hx4V','4Hx1V','4Hx2V','4Hx3V','4Hx4V']
# 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.
# First the autocorrelations for 1H->4V two consecutive HDUs with mean of FFT amplitudes , then variance
# the for the cross-correlations - the TF amplitudes are complex : in each case one HDU with real part, one with imaginary part averages
# and two for their respective variances
# NB the variance computation is optional
# at the end of the fits files, acddition HDUs contain 1D data (in reverse order) : times (center of bins), frequencies and RA
# (omitted in the case od RA-freq maps
#
# --------------------------------
# main arguments of viewtfmap
# --------------------------------
# folder = path to the fits files (w/o extension
# num = indez in the array of auto and cross-correlation (named in 'names' list above) of the quantity we want to plot/read
# mvmin/max : min and max values of the color scale (if not doable from the plot window)
# save : set tu True to save png files of the plots
# cor : to compute the correlation coefficient between 2 auto-cor channels [deprecated]
# clean : cleans previous plot (external window)
# mod : in case of cross-correlation computes the module of the complex amplitude
# (default = real part only)
# arg : in case of cross-correlation computes the argument of the complex amplitude
# cimag : to view the imaginary part
# creal : to view the real part [deprecated, default]
# cmplx :get a complete visu of the complex : color <> argument, hue <> amplitude
# hasvar : nadle varianvces if they have been stored (default)
# mtitle= string to set title to plots
# verb : verbosity level
# noplo : does not make any plot, only read and output
# 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 (
# 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]
#returns :
# times (in hours) [1d array]
# time-freq (or ra-freq) map [2d array] of demanded channel
# frequencies (MHz) [1D] ,
# ras RAs (in h)[ 1D]
# list of 1D array with slices for ewach freq in bfreq (if any)
# only if map plotted, some plot info : 3 items ax1,bnext,bprev
# mskdat : time-freq map, with masked values set to 0
#
#
#
#
#
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):
if verb :
print len(cor)
filename = folder
#print i
#print num
if noplo :
clean=False
if clean:
......@@ -35,15 +90,16 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
hdulist = fits.open(folder+'.fits')
winwid = (18,8)
winwid = (18,8) # window size
times=hdulist[len(hdulist)-1].data
ras=hdulist[len(hdulist)-2].data #41 in 43
freqs= hdulist[len(hdulist)-3].data
times=hdulist[len(hdulist)-1].data # time @ bin centers, in seconds
ras=hdulist[len(hdulist)-2].data # RA @ bin centers (hdu e.g. 41 in 43)
freqs= hdulist[len(hdulist)-3].data # frequnencies @ bin centers , MHz
if ramod :
ras=hdulist[len(hdulist)-1].data #41 in 43
freqs= hdulist[len(hdulist)-2].data
#gets the index of the HDU in the fits file from the index on names
i=2*num
if num>7 :
i= 16 + 4*(num-8)
......@@ -58,11 +114,14 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
if (verb) :
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 hds ' + str(len(hdulist))
print ' number of hdus ' + str(len(hdulist))
nfre = np.size(freqs)
ntim = np.size(times)
......@@ -85,18 +144,19 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
if (verb) :
print 'Image shape :'
print np.shape(img)
# gets what is asked for
if mod :
if num>7 :
imgi= hdulist[ i+1].data
img = np.sqrt(img**2+imgi**2)
if len(cor) ==2 :
if len(cor) ==2 : # for correlation coeff
imca = hdulist[cor[0]].data
imcb = hdulist[4+cor[1]].data
img=img/np.sqrt(imca*imcb)
if arg :
if arg :
if num>7 :
imgi= hdulist[ i+1].data
img = np.arctan2(imgi,img) * 180. / np.pi
......@@ -126,7 +186,7 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
iimg.append( hls_to_rgb(aa[ii,jj], bb[ii,jj],cc[ii,jj] ))
img = np.reshape( iimg,(np.shape(aa)[0],np.shape(aa)[1],3))
if mask :
if mask : # reads mask file
if cmplx :
print "mask & cmplx not yet compatible "
exit
......@@ -137,13 +197,14 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
mskdat = np.fabs(hdmskl[num+1].data-1.)
img = img*mskdat
filename=filename+"_masked"
if not cmplx :
# set min and max of color scale
if not cmplx : # in "normal range"
vmin = np.asarray(img[bfnor[0]:bfnor[1],]).min()
vmax = np.asarray(img[bfnor[0]:bfnor[1],]).max()
else :
vmin=0
vmax=1
#or for extrernal values
if mvmin !=0. :
vmin = mvmin
......@@ -154,6 +215,7 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
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)
#m = ax1.imshow(img,aspect='auto' ,extent=(timin,timax,frmax,frmin), interpolation='none',vmin=vmin,vmax=vmax)
if (yrcm) :
......@@ -191,7 +253,7 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
ax3.xaxis.set_view_interval( timin,timax ,ignore=True)
ax3.set_xlabel("time bin",x=0.85)
if save:
if save: # saves png with some name
ext=""
if creal :
ext = ext+"_real"
......@@ -205,7 +267,7 @@ 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])
......@@ -255,13 +317,14 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
bnext.on_clicked(cm_next)
bprev.on_clicked(cm_prev)
#----------------- END CMAP SWITCH
#plt.tight_layout()
plt.subplots_adjust(bottom=0.15)
plt.show()
if raplot :
if raplot : # does plots of signal vs RA at some frequencies
fig,ax1=plt.subplots(figsize=winwid)
[ plt.plot(ras, img[b,],label="F=%6.2f MHz" %(freqs[b]) ) for b in bfreq ]
......@@ -285,12 +348,13 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
plt.savefig(filename+"_"+ names[num]+"_TFM_ra.png")
if timeplot :
if timeplot : # does plots of sig vs time at some frequencies
fig,ax1=plt.subplots(figsize=winwid)
[ plt.plot(times/3600., img[b,],label="F=%6.2f MHz" %(freqs[b]) ) for b in bfreq ]
plt.xlabel('Time (TU, hour)')
plt.ylabel("Signal ")
plt.ylim(np.asarray(img[[bfreq],]).min(),np.asarray(img[[bfreq],]).max())
#legend with thicker linewidth
leg = plt.legend()
for legobj in leg.legendHandles:
legobj.set_linewidth(2.0)
......@@ -310,6 +374,9 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
print times.size
if (verb) :
print np.shape(img)
#return values
if noplo :
return times/3600.,img,freqs,ras,[ img[b,] for b in bfreq],mskdat
elif cmplx:
......
......@@ -54,8 +54,8 @@ def scan_source(source="Moon",
times_period = start + delta_period
# Compute the Alt-Az frames during the survey
frame_period = AltAz(obstime=times_period, location=observer)
# Compute the Alt-Az frames during the survey pressure = 0 ==> no refraction
frame_period = AltAz(obstime=times_period, location=observer, pressure=0. )
# perform the survey
source = source.lower()
......
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