| 1 | import numpy
 | 
|---|
| 2 | from asap.scantable import scantable
 | 
|---|
| 3 | from asap._asap import stgrid
 | 
|---|
| 4 | import pylab as pl
 | 
|---|
| 5 | 
 | 
|---|
| 6 | class asapgrid:
 | 
|---|
| 7 |     def __init__( self, infile ):
 | 
|---|
| 8 |         self.infile = infile
 | 
|---|
| 9 |         self.outfile = None
 | 
|---|
| 10 |         self.gridder = stgrid( self.infile )
 | 
|---|
| 11 | 
 | 
|---|
| 12 |     def setData( self, infile ):
 | 
|---|
| 13 |         self.gridder._setin( infile )
 | 
|---|
| 14 | 
 | 
|---|
| 15 |     def setPolList( self, pollist ):
 | 
|---|
| 16 |         self.gridder._setpollist( pollist )
 | 
|---|
| 17 | 
 | 
|---|
| 18 |     def defineImage( self, nx=-1, ny=-1, cellx='', celly='', center='' ):
 | 
|---|
| 19 |         self.gridder._defineimage( nx, ny, cellx, celly, center )
 | 
|---|
| 20 | 
 | 
|---|
| 21 |     def setOption( self, convType='box', convSupport=-1 ):
 | 
|---|
| 22 |         self.gridder._setoption( convType, convSupport )
 | 
|---|
| 23 | 
 | 
|---|
| 24 |     def grid( self ):
 | 
|---|
| 25 |         self.gridder._grid()
 | 
|---|
| 26 | 
 | 
|---|
| 27 |     def save( self, outfile='' ):
 | 
|---|
| 28 |         self.outfile = self.gridder._save( outfile ) 
 | 
|---|
| 29 | 
 | 
|---|
| 30 |     def plot( self, plotchan=-1, plotpol=-1 ):
 | 
|---|
| 31 |         plotter = _SDGridPlotter( self.infile, self.outfile )
 | 
|---|
| 32 |         plotter.plot( chan=plotchan, pol=plotpol )
 | 
|---|
| 33 |         
 | 
|---|
| 34 | class _SDGridPlotter:
 | 
|---|
| 35 |     def __init__( self, infile, outfile=None ):
 | 
|---|
| 36 |         self.infile = infile
 | 
|---|
| 37 |         self.outfile = outfile
 | 
|---|
| 38 |         if self.outfile is None:
 | 
|---|
| 39 |             self.outfile = self.infile.rstrip('/')+'.grid'
 | 
|---|
| 40 |         self.grid = None
 | 
|---|
| 41 |         self.pointing = None
 | 
|---|
| 42 |         self.data = None
 | 
|---|
| 43 |         self.nx = -1
 | 
|---|
| 44 |         self.ny = -1
 | 
|---|
| 45 |         self.nchan = 0
 | 
|---|
| 46 |         self.npol = 0
 | 
|---|
| 47 |         self.pollist = []
 | 
|---|
| 48 |         self.cellx = 0.0
 | 
|---|
| 49 |         self.celly = 0.0
 | 
|---|
| 50 |         self.center = [0.0,0.0]
 | 
|---|
| 51 |         self.nonzero = [[0.0],[0.0]]
 | 
|---|
| 52 |         self.get()
 | 
|---|
| 53 | 
 | 
|---|
| 54 |     def get( self ):
 | 
|---|
| 55 |         s = scantable( self.infile, average=False )
 | 
|---|
| 56 |         self.pointing = numpy.array( s.get_directionval() ).transpose()
 | 
|---|
| 57 |         spectra = []
 | 
|---|
| 58 |         for i in xrange(s.nrow()):
 | 
|---|
| 59 |             spectra.append( s._getspectrum( i ) )
 | 
|---|
| 60 |         spectra = numpy.array( spectra ).transpose()
 | 
|---|
| 61 |         self.nchan = spectra.shape[0]
 | 
|---|
| 62 |         del s
 | 
|---|
| 63 | 
 | 
|---|
| 64 |         s = scantable( self.outfile, average=False )
 | 
|---|
| 65 |         nrow = s.nrow()
 | 
|---|
| 66 |         pols = numpy.ones( nrow, dtype=int )
 | 
|---|
| 67 |         for i in xrange(nrow):
 | 
|---|
| 68 |             pols[i] = s.getpol(i)
 | 
|---|
| 69 |         self.pollist, indices = numpy.unique( pols, return_inverse=True )
 | 
|---|
| 70 |         self.npol = len(self.pollist)
 | 
|---|
| 71 |         self.pollist = self.pollist[indices[:self.npol]]
 | 
|---|
| 72 |         #print 'pollist=',self.pollist
 | 
|---|
| 73 |         #print 'npol=',self.npol
 | 
|---|
| 74 |         #print 'nrow=',nrow
 | 
|---|
| 75 |         dirstring = numpy.array(s.get_direction()).take(range(0,nrow,self.npol))
 | 
|---|
| 76 |         self.grid = numpy.array( s.get_directionval() ).take(range(0,nrow,self.npol),axis=0).transpose()
 | 
|---|
| 77 |         spectra = numpy.zeros( (self.npol,self.nchan,nrow/self.npol), dtype=float )
 | 
|---|
| 78 |         irow = 0 
 | 
|---|
| 79 |         for i in xrange(nrow/self.npol):
 | 
|---|
| 80 |             for ip in xrange(self.npol):
 | 
|---|
| 81 |                 spectra[ip,:,i] = s._getspectrum( irow )
 | 
|---|
| 82 |                 irow += 1
 | 
|---|
| 83 | 
 | 
|---|
| 84 |         idx = 0
 | 
|---|
| 85 |         d0 = dirstring[0].split()[-1]
 | 
|---|
| 86 |         while ( dirstring[idx].split()[-1] == d0 ):  
 | 
|---|
| 87 |             idx += 1
 | 
|---|
| 88 |         self.ny = idx
 | 
|---|
| 89 |         self.nx = nrow / (self.npol * idx )
 | 
|---|
| 90 |         #print 'nx,ny=',self.nx,self.ny
 | 
|---|
| 91 |         
 | 
|---|
| 92 |         self.cellx = abs( self.grid[0][0] - self.grid[0][1] )
 | 
|---|
| 93 |         self.celly = abs( self.grid[1][0] - self.grid[1][self.ny] )
 | 
|---|
| 94 |         #print 'cellx,celly=',self.cellx,self.celly
 | 
|---|
| 95 | 
 | 
|---|
| 96 |         self.data = spectra.reshape( (self.npol,self.nchan,self.nx,self.ny) )
 | 
|---|
| 97 | 
 | 
|---|
| 98 |     def plot( self, chan=-1, pol=-1 ):
 | 
|---|
| 99 |         if pol < 0:
 | 
|---|
| 100 |             data = self.data.mean(axis=0)
 | 
|---|
| 101 |             opt = 'averaged over pol'
 | 
|---|
| 102 |         else:
 | 
|---|
| 103 |             idx = self.pollist.tolist().index( pol )
 | 
|---|
| 104 |             #print 'idx=',idx
 | 
|---|
| 105 |             data = self.data[idx]
 | 
|---|
| 106 |             opt = 'pol %s'%(pol)
 | 
|---|
| 107 |         if chan < 0:
 | 
|---|
| 108 |             data = data.mean(axis=0)
 | 
|---|
| 109 |             opt += ', averaged over channel'
 | 
|---|
| 110 |         else:
 | 
|---|
| 111 |             data = data[chan]
 | 
|---|
| 112 |             opt += ', channel %s'%(chan)
 | 
|---|
| 113 |         title = 'Gridded Image (%s)'%(opt)
 | 
|---|
| 114 |         pl.figure(10)
 | 
|---|
| 115 |         pl.clf()
 | 
|---|
| 116 |         pl.plot(self.grid[0],self.grid[1],'.',color='blue')
 | 
|---|
| 117 |         pl.plot(self.pointing[0],self.pointing[1],'.',color='red')
 | 
|---|
| 118 |         extent=[self.grid[0].min()-0.5*self.cellx,
 | 
|---|
| 119 |                 self.grid[0].max()+0.5*self.cellx,
 | 
|---|
| 120 |                 self.grid[1].min()-0.5*self.celly,
 | 
|---|
| 121 |                 self.grid[1].max()+0.5*self.celly]
 | 
|---|
| 122 |         pl.imshow(data,extent=extent,origin='lower',interpolation='nearest')
 | 
|---|
| 123 |         pl.colorbar()
 | 
|---|
| 124 |         pl.xlabel('R.A. [rad]')
 | 
|---|
| 125 |         pl.ylabel('Dec. [rad]')
 | 
|---|
| 126 |         pl.title( title )
 | 
|---|