Commit 4322001d authored by Jean-Eric Campagne's avatar Jean-Eric Campagne
Browse files

Merge branch 'master' of gitlab.in2p3.fr:baoradio/AnaPAON4

parents 386171d2 5b1c1e2d
# The Normalize class is largely based on code provided by Sarah Graves.
import numpy as np
import numpy.ma as ma
import matplotlib.cbook as cbook
from matplotlib.colors import Normalize
class MyNormalize(Normalize):
'''
A Normalize class for imshow that allows different stretching functions
for astronomical images.
'''
def __init__(self, stretch='linear', exponent=5, vmid=None, vmin=None,
vmax=None, clip=False):
'''
Initalize an APLpyNormalize instance.
Optional Keyword Arguments:
*vmin*: [ None | float ]
Minimum pixel value to use for the scaling.
*vmax*: [ None | float ]
Maximum pixel value to use for the scaling.
*stretch*: [ 'linear' | 'log' | 'sqrt' | 'arcsinh' | 'power' ]
The stretch function to use (default is 'linear').
*vmid*: [ None | float ]
Mid-pixel value used for the log and arcsinh stretches. If
set to None, a default value is picked.
*exponent*: [ float ]
if self.stretch is set to 'power', this is the exponent to use.
*clip*: [ True | False ]
If clip is True and the given value falls outside the range,
the returned value will be 0 or 1, whichever is closer.
'''
if vmax < vmin:
raise Exception("vmax should be larger than vmin")
# Call original initalization routine
Normalize.__init__(self, vmin=vmin, vmax=vmax, clip=clip)
# Save parameters
self.stretch = stretch
self.exponent = exponent
if stretch == 'power' and np.equal(self.exponent, None):
raise Exception("For stretch=='power', an exponent should be specified")
if np.equal(vmid, None):
if stretch == 'log':
if vmin > 0:
self.midpoint = vmax / vmin
else:
raise Exception("When using a log stretch, if vmin < 0, then vmid has to be specified")
elif stretch == 'arcsinh':
self.midpoint = -1. / 30.
else:
self.midpoint = None
else:
if stretch == 'log':
if vmin < vmid:
raise Exception("When using a log stretch, vmin should be larger than vmid")
self.midpoint = (vmax - vmid) / (vmin - vmid)
elif stretch == 'arcsinh':
self.midpoint = (vmid - vmin) / (vmax - vmin)
else:
self.midpoint = None
def __call__(self, value, clip=None):
#read in parameters
method = self.stretch
exponent = self.exponent
midpoint = self.midpoint
# ORIGINAL MATPLOTLIB CODE
if clip is None:
clip = self.clip
if cbook.iterable(value):
vtype = 'array'
val = ma.asarray(value).astype(np.float)
else:
vtype = 'scalar'
val = ma.array([value]).astype(np.float)
self.autoscale_None(val)
vmin, vmax = self.vmin, self.vmax
if vmin > vmax:
raise ValueError("minvalue must be less than or equal to maxvalue")
elif vmin == vmax:
return 0.0 * val
else:
if clip:
mask = ma.getmask(val)
val = ma.array(np.clip(val.filled(vmax), vmin, vmax),
mask=mask)
result = (val - vmin) * (1.0 / (vmax - vmin))
# CUSTOM APLPY CODE
# Keep track of negative values
negative = result < 0.
if self.stretch == 'linear':
pass
elif self.stretch == 'log':
result = ma.log10(result * (self.midpoint - 1.) + 1.) \
/ ma.log10(self.midpoint)
elif self.stretch == 'sqrt':
result = ma.sqrt(result)
elif self.stretch == 'arcsinh':
result = ma.arcsinh(result / self.midpoint) \
/ ma.arcsinh(1. / self.midpoint)
elif self.stretch == 'power':
result = ma.power(result, exponent)
else:
raise Exception("Unknown stretch in APLpyNormalize: %s" %
self.stretch)
# Now set previously negative values to 0, as these are
# different from true NaN values in the FITS image
result[negative] = -np.inf
if vtype == 'scalar':
result = result[0]
return result
def inverse(self, value):
# ORIGINAL MATPLOTLIB CODE
if not self.scaled():
raise ValueError("Not invertible until scaled")
vmin, vmax = self.vmin, self.vmax
# CUSTOM APLPY CODE
if cbook.iterable(value):
val = ma.asarray(value)
else:
val = value
if self.stretch == 'linear':
pass
elif self.stretch == 'log':
val = (ma.power(10., val * ma.log10(self.midpoint)) - 1.) / (self.midpoint - 1.)
elif self.stretch == 'sqrt':
val = val * val
elif self.stretch == 'arcsinh':
val = self.midpoint * \
ma.sinh(val * ma.arcsinh(1. / self.midpoint))
elif self.stretch == 'power':
val = ma.power(val, (1. / self.exponent))
else:
raise Exception("Unknown stretch in APLpyNormalize: %s" %
self.stretch)
return vmin + val * (vmax - vmin)
import math as math
import sys, platform, os
import numpy as np
from astropy.io import fits
from matplotlib import pyplot as plt
import itertools as itert
from skimage import data, img_as_float
from skimage import exposure
import mynormalize as mynorm
from matplotlib.widgets import Slider, Button, RadioButtons
import time
# Choices: ['auto', 'gtk', 'gtk3', 'inline', 'nbagg', 'notebook', 'osx', 'qt', 'qt4', 'qt5', 'tk', 'wx']
rfreq = [1370.,1380.,1421.,1440.]
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']
def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=False,mtitle="",clean=True,mod=False,arg=False,noplo=False,verb=False,yrcm=False,cor=[],creal=False,cimag=False):
print len(cor)
filename = folder
i=2*num
if num>7 :
i= 16 + 4*(num-8)
#print i
#print num
if clean:
plt.close("all")
hdulist = fits.open(folder+'.fits')
if (verb) :
print hdulist[i].header
freqs= hdulist[len(hdulist)-3].data
bfreq = [ (np.abs(freqs-rf)).argmin() for rf in rfreq ]
len(hdulist)
ras=hdulist[len(hdulist)-2].data #41 in 43
times=hdulist[len(hdulist)-1].data
nfre = np.size(freqs)
ntim = np.size(times)
timin = min(times)/3600.
timax = max(times)/3600.
frmin = min(freqs)
frmax=max(freqs)
if (not noplo):
fig,ax1=plt.subplots(figsize=(20,9))
img = hdulist[ i].data
if(cimag):
img = hdulist[ i+1].data
print np.shape(img)
if mod :
if num>7 :
imgi= hdulist[ i+1].data
img = sqrt(img**2+imgi**2)
if len(cor) ==2 :
imca = hdulist[cor[0]].data
imcb = hdulist[4+cor[1]].data
img=img/sqrt(imca*imcb)
if arg :
if num>7 :
imgi= hdulist[ i+1].data
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 "
exit
vmin = np.asarray(img[300:900,]).min()
if mvmin !=0. :
vmin = mvmin
vmax = np.asarray(img[300:900,]).min()*10.
if mvmax !=0. :
vmax = mvmax
if verb :
print vmin, np.asarray(img[300:900,]).min()
print vmax,np.asarray(img[300:900,]).max()
if (not noplo) :
im = ax1.imshow(img,aspect='auto' ,extent=(timin,timax,frmax,frmin), interpolation='none',vmin=vmin,vmax=vmax)
if (yrcm) :
im.set_cmap('YlOrBr_r')
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)
plt.show()
if mtitle =="":
plt.title(filename + ' ' + names[num],y=1.08)
else :
plt.title(filename + ' ' +mtitle,y=1.08)
plt.ylabel("Frequency (Mhz)")
plt.xlabel("Temps (h)")
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(len(freqs),0.,ignore=True)
ax2.set_ylabel("Freq. bin",rotation=270.,y=0.9)
ax3=ax1.twiny()
ax3.xaxis.set_view_interval(0.,len(times))
ax3.set_xlabel("time bin",x=0.85)
if save:
ext=""
if creal :
ext = ext+"_real"
plt.savefig(filename+"_"+ names[num]+"_rawTFM.png")
axcolor = 'lightgoldenrodyellow'
axmin = fig.add_axes([0.075, 0.05, 0.35, 0.03], facecolor = axcolor )
axmax = fig.add_axes([0.075, 0.02, 0.35, 0.03], facecolor = axcolor)
# axrange = fig.add_axes([0.6, 0.05, 0.35, 0.03], facecolor = axcolor)
# axcent = fig.add_axes([0.6, 0.02, 0.35, 0.03], facecolor = axcolor)
smin = Slider(axmin, 'Min', vmin,vmax, valinit=vmin)
smax = Slider(axmax, 'Max', vmin,vmax, valinit=vmax)
def update(val):
print "min=" + str(smin.val)
print "max=" + str(smax.val)
cbar.norm.vmin=smin.val
cbar.norm.vmax=smax.val
start = time.time()
cbar.draw_all()
#elp = time.time()
#print str(elp-start) + " seconds"
im.set_norm(cbar.norm)
#elp = time.time()
#print str(elp-start) + " seconds"
cbar.patch.figure.canvas.draw()
#elp = time.time()
#print "end->"+str(elp-start) + " seconds"
#cen.set_val( (smax.val + smin.val)/2.)
#rng.set_val( (smax.val - smin.val) )
#m.set_clim([smin.val,smax.val])
#ig.canvas.draw()
smin.on_changed(update)
smax.on_changed(update)
# class Index(object):
# ind = cm_cycle.index(cbar.get_cmap().name)
# cycle = sorted([i for i in dir(plt.cm) if hasattr(getattr(plt.cm,i),'N')])
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])
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)
cm_index = cm_index+1
if cm_index>= len(cm_cycle):
cm_index=0
change_cm(cm_index)
#print cm_index
def cm_prev(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)
cm_index = cm_index-1
if cm_index < 0 :
cm_index=len(cm_cycle)-1
change_cm(cm_index)
#print cm_index
def change_cm(index):
cm_cycle = sorted([i for i in dir(plt.cm) if hasattr(getattr(plt.cm,i),'N')])
cmap = cm_cycle[index]
#print cmap
cbar.set_cmap(cmap)
cbar.draw_all()
im.set_cmap(cmap)
cbar.patch.figure.canvas.draw()
bcur.label.set_text(cmap)
indprev=index-1
if indprev<0 :
indprev=len(cm_cycle)
indnext=index+1
if indnext>=len(cm_cycle):
indnext=0
bnext.label.set_text(cm_cycle[indnext])
bprev.label.set_text(cm_cycle[indprev])
bnext.on_clicked(cm_next)
bprev.on_clicked(cm_prev)
plt.tight_layout()
plt.subplots_adjust(bottom=0.15)
plt.show()
if raplot :
fig,ax1=plt.subplots(figsize=(20,8))
[ plt.plot(ras, img[b,],label="F=%6.2f MHz" %(freqs[b]) ) for b in bfreq ]
plt.xlabel("RA (hour)")
plt.ylabel("Signal ")
plt.ylim(np.asarray(img[[bfreq],]).min(),np.asarray(img[[bfreq],]).max())
leg = plt.legend()
for legobj in leg.legendHandles:
legobj.set_linewidth(2.0)
if mtitle =="":
plt.title(filename + ' ' + names[num],y=1.08)
else :
plt.title(filename + ' ' +mtitle,y=1.08)
if save:
plt.savefig(filename+"_"+ names[num]+"_rawTFM_freq.png")
if timeplot :
fig,ax1=plt.subplots(figsize=(20,8))
[ 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())
leg = plt.legend()
for legobj in leg.legendHandles:
legobj.set_linewidth(2.0)
if mtitle =="":
plt.title(filename + ' ' + names[num],y=1.08)
else :
plt.title(filename + ' ' +mtitle,y=1.08)
if save:
plt.savefig(filename+"_"+ names[num]+"_rawTFM_freq.png")
if (verb) :
print times.size
if (verb) :
print np.shape(img)
if noplo :
return times/3600.,img,freqs,ras,[ img[b,] for b in bfreq]
else :
return times/3600.,img,freqs,ras,[ img[b,] for b in bfreq],ax1,bnext,bprev#smin,smax
......@@ -37,3 +37,8 @@ norm g1H g2H g3H g4H g1V g2V g3V g4V
## @freqband CentralFreqMHz BandWidthMHz
@freqband 1390. 1
@freqband 1420.4 0.25
### Select VV and/or HV cross correlation processing
@dovv
@dohv
### Require computation of dispersion (sigma) TFM maps, in addition to mean (average) TFM maps
@dosigma
......@@ -10,13 +10,13 @@ MYINCLISTHERE = p4autils.h visip4reader.h visp4winreader.h p4gnugain.h p4gvcor.h
MYOLISTHERE = $(OBJ)/p4autils.o $(OBJ)/visip4reader.o $(OBJ)/visp4winreader.o $(OBJ)/p4gnugain.o $(OBJ)/p4gvcor.o $(OBJ)/p4phcor.o
# Define our target list
all : rdvisip4 visi2ntac visi2dtacx visi2tmfreq tstutls p4conv2fits msvis2dt visiavg filt_blind ckconvphases vis2ra
all : rdvisip4 visi2ntac visi2dtacx visi2tmfreq tstutls p4conv2fits msvis2dt visiavg filt_blind ckconvphases vis2ra p4vdblist
clean :
rm -f $(EXE)/rdvisip4 $(EXE)/visi2ntac $(EXE)/visi2dtacx $(EXE)/visi2tmfreq $(EXE)/p4conv2fits $(EXE)/msvis2dt $(EXE)/visiavg \
$(EXE)/filt_blind $(EXE)/ckconvphases $(EXE)/vis2ra
$(EXE)/filt_blind $(EXE)/ckconvphases $(EXE)/vis2ra $(EXE)/p4vdblist
rm -f $(OBJ)/rdvisip4.o $(OBJ)/visi2ntac.o $(OBJ)/visi2dtacx.o $(OBJ)/visi2tmfreq.o $(OBJ)/p4conv2fits.o $(OBJ)/msvis2dt.o $(OBJ)/visiavg.o \
$(OBJ)/filt_blind.o $(OBJ)/ckconvphases.o $(OBJ)/vis2ra.o
$(OBJ)/filt_blind.o $(OBJ)/ckconvphases.o $(OBJ)/vis2ra.o $(OBJ)/p4vdblist.o
rm -f $(MYOLISTHERE)
depend :
......@@ -136,6 +136,16 @@ $(OBJ)/rdvisip4.o : rdvisip4.cc $(MYINCLISTHERE)
######
## programme de tests de base pour diverses classes / fonctions
p4vdblist : $(EXE)/p4vdblist
echo '---p4vdblist made'
$(EXE)/p4vdblist : $(OBJ)/p4vdblist.o $(MYOLISTHERE)
$(CXXLINK) -o $(EXE)/p4vdblist $(OBJ)/p4vdblist.o $(MYOLISTHERE) $(SOPHYAEXTSLBLIST)
$(OBJ)/p4vdblist.o : tstutls.cc $(MYINCLISTHERE)
$(CXXCOMPILE) -o $(OBJ)/p4vdblist.o p4vdblist.cc
######
## programme pour lister la description des fichiers de la base de donnees de Visibilites
tstutls : $(EXE)/tstutls
echo '---tstutls made'
......
......@@ -194,6 +194,17 @@ std::vector<sa_size_t> P4AVisiNumEncoder::getAllVCrossCor()
}
return P4AVisiNumEncoder::Convert2VisiNumber(nms);
}
std::vector<sa_size_t> P4AVisiNumEncoder::getAllHVCrossCor()
{
const char * vcxnm[16]={"1H-1V","1H-2V","1H-3V","1H-4V","2H-1V","2H-2V","2H-3V","2H-4V","3H-1V","3H-2V","3H-3V","3H-4V","4H-1V","4H-2V","4H-3V","4H-4V"};
std::vector<std::string> nms;
for(size_t i=0; i<16; i++) {
nms.push_back(std::string(vcxnm[i]));
}
return P4AVisiNumEncoder::Convert2VisiNumber(nms);
}
//---------------------------------------------------------------------
//------------------------ classe P4AnaParams ------------------------
//---------------------------------------------------------------------
......@@ -203,7 +214,8 @@ P4AnaParams::P4AnaParams()
gain_gnu_file_(""), gain_var_file_(""), gain_blind_filt_file_(""), phase_corr_param_file_(""),
outfile_("p4a.ppf"), fitsoutfile_(""), prtlev_(0), prtmodulo_(1),
Nmean_(10), fgTFM_(false), TFMtimebin_(30), TFMfreqbin_(4), fgdtable_(false),
fgRA_(false), RAcenter_(12.), RAwidth_(24.), RAresol_mn_(1.)
fgRA_(false), RAcenter_(12.), RAwidth_(24.), RAresol_mn_(1.),
doVV_(false), doHV_(false), doSigma_(false)
{
}
......@@ -218,6 +230,7 @@ P4AnaParams::P4AnaParams(P4AnaParams const & a)
Nmean_(a.Nmean_), fgTFM_(a.fgTFM_), TFMtimebin_(a.TFMtimebin_), TFMfreqbin_(a.TFMfreqbin_),
fgdtable_(a.fgdtable_), fbands_(a.fbands_),
fgRA_(a.fgRA_), RAcenter_(a.RAcenter_), RAwidth_(a.RAwidth_), RAresol_mn_(a.RAresol_mn_),
doVV_(a.doVV_), doHV_(a.doHV_), doSigma_(a.doSigma_),
lastargs_(a.lastargs_)
{
}
......@@ -235,6 +248,7 @@ P4AnaParams& P4AnaParams::operator = (P4AnaParams const & a)
Nmean_=a.Nmean_; fgTFM_=a.fgTFM_; TFMfreqbin_=a.TFMfreqbin_; TFMtimebin_=a.TFMtimebin_;
fgdtable_=a.fgdtable_; fbands_=a.fbands_;
fgRA_=a.fgRA_; RAcenter_=a.RAcenter_; RAwidth_=a.RAwidth_; RAresol_mn_=a.RAresol_mn_;
doVV_=a.doVV_; doHV_=a.doHV_; doSigma_=a.doSigma_;
lastargs_=a.lastargs_;
return *this;
}
......@@ -275,6 +289,18 @@ int P4AnaParams::DecodeArgs(int narg, const char* arg[])
else if (fbo=="-reorderfreq") {
fgreorderfreq_=true; arg++; narg--; lastargs_.clear();
}
else if (fbo=="-dovv") {
doVV_=true; arg++; narg--; lastargs_.clear();
}
else if (fbo=="-dohv") {
doHV_=true; arg++; narg--; lastargs_.clear();
}
else if (fbo=="-dosigma") {
doSigma_=true; arg++; narg--; lastargs_.clear();
}
else if (fbo=="-dohv") {
doHV_=true; arg++; narg--; lastargs_.clear();
}
else if (fbo=="-gnu") {
if (narg<2) { cout << "P4AnaParams::DecodeArgs missing/bad argument, -h for help " << endl; return -1; }
gain_gnu_file_=arg[2]; arg+=2; narg-=2; lastargs_.clear();
......@@ -293,7 +319,7 @@ int P4AnaParams::DecodeArgs(int narg, const char* arg[])
}
else if (fbo=="-norm") {
if (narg<2) { cout << "P4AnaParams::DecodeArgs missing/bad argument, -h for help " << endl; return -1; }
normalisation_.resize((size_t)P4AVisiNumEncoder::getTotalNbVisibilities());
normalisation_.resize((size_t)P4AVisiNumEncoder::getTotalNbFeeds());
for(size_t jj=0; jj<normalisation_.size(); jj++) normalisation_[jj]=1.;
string sn=arg[2]; vector<string> vs;
SplitStringToVString(sn,vs,',');
......@@ -440,6 +466,15 @@ int P4AnaParams::ReadDatacardFile(const char * dcfilename)
prtlev_=dc.IParam("prtlev",0,1); prtmodulo_=dc.IParam("prtlev",1,1);
cnt++;
}
if (dc.HasKey("dovv")) { //@dovv : process also VV cross-correlations
doVV_=true; cnt++;
}
if (dc.HasKey("dohv")) { //@dohv : process also HV cross-correlations
doHV_=true; cnt++;
}
if (dc.HasKey("dosigma")) { //@dosigma : computes also dispersion (Sigma) TFM maps
doSigma_=true; cnt++;
}
for (size_t i=0; i<dc_fbands_.size(); i++) fbands_.push_back(dc_fbands_[i]);
cout << " P4AnaParams::ReadDatacardFile(FileName=" << dcfilename << " ) -> decoded " << cnt+dc_fbands_.size()<<" cards"<<endl;
dc_fbands_.clear();
......@@ -453,7 +488,9 @@ int P4AnaParams::UsageOptions()
<< " [-inrange Imin,Imax[,Istep] ] [-inall] [-inseltm YYYY-MM-DDThh:mm:ss duration_minutes] \n"
<< " [-gnu gain_filename] [-vgf filename] [-bff filename] [-phcor filename] [-norm fac1H,fac2H...,fac4V] \n"
<< " [-o outfilename] [-fo fits_filename] [-nmean N] [-reorderfreq] \n"
<< " [-tfm timebin,freqbin] [-ra RAcenterHours,RAwidthHours,RAresolMinutes] [-freq f0,df] [-prt level] \n"
<< " [-tfm timebin,freqbin] [-ra RAcenterHours,RAwidthHours,RAresolMinutes] [-freq f0,df] \n"
<< " [-dovv] [-dohv] [-dosigma] [-prt level] \n"
<< " -dcf DatacardFileName : read and decode parameters from datacard file \n"
<< " -in5/6 InputPath_BAO5/6 : input file path for BAO5 and BAO6 \n"
<< " -inrange Imin,Imax[,Istep] : define range of visibilite matrices to be processed (Imin<=I<=Imax with step Istep) \n"
......@@ -474,6 +511,9 @@ int P4AnaParams::UsageOptions()
<< " -ra RAcenterHours,RAwidthHours,RAresolMinutes: define RA (right ascension) map center,total width (hours), resolution (minutes) \n"
<< " -freq f0,df : definition of a frequency band centered on f0, with total width df (in MHz) \n"
<< " several frequency bands can be specified by mutiple -freq options \n"
<< " -dovv : processing also VV cross correlations \n"
<< " -dohv : processing also HV cross correlations \n"
<< " -dosigma : compute also TFM maps of dispersion (sigma) \n"