add a script to view tf maps

Python/mynormalize.py

+# 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)

Python/viewtfmap.py

+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')
+        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], axisbg = axcolor )
+        axmax = fig.add_axes([0.075, 0.02, 0.35, 0.03], axisbg = axcolor)
+#       axrange = fig.add_axes([0.6, 0.05, 0.35, 0.03], axisbg = axcolor)
+#       axcent = fig.add_axes([0.6, 0.02, 0.35, 0.03], axisbg = 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')])
+        
+        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)
+
+            change_cm(cm_index)
+            print cm_index
+            
+        def change_cm(index):
+            print "test "+str(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()
+
+        axnext = fig.add_axes([0.6, 0.04, 0.05, 0.05])
+        bnext = Button(axnext, 'Next cm')
+        axprev = fig.add_axes([0.7, 0.04, 0.05, 0.05])
+        bprev = Button(axprev, 'Prev. cm')
+        
+#       srng.on_changed(update_rng)
+#       scen.on_changed(update_rng)
+        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  
+
+
+ 
+