source: trunk/python/asapgrid.py@ 2383

Last change on this file since 2383 was 2375, checked in by Takeshi Nakazato, 13 years ago

New Development: No

JIRA Issue: Yes CAS-2816

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Several bug fixes and updates

  • fixed misunderstanding of array shape when getting data
  • load flag data as Int (originally uChar/uInt) at the beginning of gridding
  • do not plot observed position and grid position by default
  • use Array (instead of Vector/Matrix/Cube) as much as possible


File size: 7.2 KB
Line 
1import numpy
2from asap import rcParams
3from asap.scantable import scantable
4from asap.selector import selector
5from asap._asap import stgrid
6import pylab as pl
7from logging import asaplog
8
9class asapgrid:
10 def __init__( self, infile ):
11 self.infile = infile
12 self.outfile = None
13 self.gridder = stgrid( self.infile )
14 self.ifno = None
15
16 def setData( self, infile ):
17 self.gridder._setin( infile )
18
19 def setIF( self, ifno ):
20 self.ifno = ifno
21 self.gridder._setif( self.ifno )
22
23 def setPolList( self, pollist ):
24 self.gridder._setpollist( pollist )
25
26 def setScanList( self, scanlist ):
27 self.gridder._setscanlist( scanlist )
28
29 def defineImage( self, nx=-1, ny=-1, cellx='', celly='', center='' ):
30 self.gridder._defineimage( nx, ny, cellx, celly, center )
31
32 def setFunc( self, func='box', width=-1 ):
33 self.gridder._setfunc( func, width )
34
35 def setWeight( self, weightType='uniform' ):
36 self.gridder._setweight( weightType )
37
38 def grid( self ):
39 self.gridder._grid()
40
41 def save( self, outfile='' ):
42 self.outfile = self.gridder._save( outfile )
43
44 def plot( self, plotchan=-1, plotpol=-1, plotobs=False, plotgrid=False ):
45 import time
46 t0=time.time()
47 # to load scantable on disk
48 storg = rcParams['scantable.storage']
49 rcParams['scantable.storage'] = 'disk'
50 plotter = _SDGridPlotter( self.infile, self.outfile, self.ifno )
51 plotter.plot( chan=plotchan, pol=plotpol, plotobs=plotobs, plotgrid=plotgrid )
52 # back to original setup
53 rcParams['scantable.storage'] = storg
54 t1=time.time()
55 asaplog.push('plot: elapsed time %s sec'%(t1-t0))
56 asaplog.post('DEBUG','asapgrid.plot')
57
58class _SDGridPlotter:
59 def __init__( self, infile, outfile=None, ifno=-1 ):
60 self.infile = infile
61 self.outfile = outfile
62 if self.outfile is None:
63 self.outfile = self.infile.rstrip('/')+'.grid'
64 self.nx = -1
65 self.ny = -1
66 self.nchan = 0
67 self.npol = 0
68 self.pollist = []
69 self.cellx = 0.0
70 self.celly = 0.0
71 self.center = [0.0,0.0]
72 self.nonzero = [[0.0],[0.0]]
73 self.ifno = ifno
74 self.tablein = None
75 self.nrow = 0
76 self.blc = None
77 self.trc = None
78 self.get()
79
80 def get( self ):
81 self.tablein = scantable( self.infile, average=False )
82 if self.ifno < 0:
83 ifno = self.tablein.getif(0)
84 print 'ifno=',ifno
85 else:
86 ifno = self.ifno
87 sel = selector()
88 sel.set_ifs( ifno )
89 self.tablein.set_selection( sel )
90 self.nchan = len(self.tablein._getspectrum(0))
91 self.nrow = self.tablein.nrow()
92 del sel
93
94 s = scantable( self.outfile, average=False )
95 nrow = s.nrow()
96 pols = numpy.ones( nrow, dtype=int )
97 for i in xrange(nrow):
98 pols[i] = s.getpol(i)
99 self.pollist, indices = numpy.unique( pols, return_inverse=True )
100 self.npol = len(self.pollist)
101 self.pollist = self.pollist[indices[:self.npol]]
102 #print 'pollist=',self.pollist
103 #print 'npol=',self.npol
104 #print 'nrow=',nrow
105
106 idx = 0
107 d0 = s.get_direction( 0 ).split()[-1]
108 while ( s.get_direction(self.npol*idx).split()[-1] == d0 ):
109 idx += 1
110
111 self.nx = idx
112 self.ny = nrow / (self.npol * idx )
113 #print 'nx,ny=',self.nx,self.ny
114
115 self.blc = s.get_directionval( 0 )
116 self.trc = s.get_directionval( nrow-self.npol )
117 #print self.blc
118 #print self.trc
119 incrx = s.get_directionval( self.npol )
120 incry = s.get_directionval( self.nx*self.npol )
121 self.cellx = abs( self.blc[0] - incrx[0] )
122 self.celly = abs( self.blc[1] - incry[1] )
123 #print 'cellx,celly=',self.cellx,self.celly
124
125 def plot( self, chan=-1, pol=-1, plotobs=False, plotgrid=False ):
126 if pol < 0:
127 opt = 'averaged over pol'
128 else:
129 opt = 'pol %s'%(pol)
130 if chan < 0:
131 opt += ', averaged over channel'
132 else:
133 opt += ', channel %s'%(chan)
134 data = self.getData( chan, pol )
135 title = 'Gridded Image (%s)'%(opt)
136 pl.figure(10)
137 pl.clf()
138 # plot grid position
139 if plotgrid:
140 x = numpy.arange(self.blc[0],self.trc[0]+0.5*self.cellx,self.cellx,dtype=float)
141 #print 'len(x)=',len(x)
142 #print 'x=',x
143 ybase = numpy.ones(self.nx,dtype=float)*self.blc[1]
144 #print 'len(ybase)=',len(ybase)
145 incr = self.celly
146 for iy in xrange(self.ny):
147 y = ybase + iy * incr
148 #print y
149 pl.plot(x,y,',',color='blue')
150 # plot observed position
151 if plotobs:
152 irow = 0
153 while ( irow < self.nrow ):
154 chunk = self.getPointingChunk( irow )
155 #print chunk
156 pl.plot(chunk[0],chunk[1],',',color='green')
157 irow += chunk.shape[1]
158 #print irow
159 # show image
160 extent=[self.blc[0]-0.5*self.cellx,
161 self.trc[0]+0.5*self.cellx,
162 self.blc[1]-0.5*self.celly,
163 self.trc[1]+0.5*self.celly]
164 pl.imshow(data,extent=extent,origin='lower',interpolation='nearest')
165 pl.colorbar()
166 pl.xlabel('R.A. [rad]')
167 pl.ylabel('Dec. [rad]')
168 pl.title( title )
169
170 def getPointingChunk( self, irow ):
171 numchunk = 1000
172 nrow = min( self.nrow-irow, numchunk )
173 #print 'nrow=',nrow
174 v = numpy.zeros( (2,nrow), dtype=float )
175 idx = 0
176 for i in xrange(irow,irow+nrow):
177 d = self.tablein.get_directionval( i )
178 v[0,idx] = d[0]
179 v[1,idx] = d[1]
180 idx += 1
181 return v
182
183 def getData( self, chan=-1, pol=-1 ):
184 if chan == -1:
185 spectra = self.__chanAverage()
186 else:
187 spectra = self.__chanIndex( chan )
188 data = spectra.reshape( (self.npol,self.ny,self.nx) )
189 if pol == -1:
190 retval = data.mean(axis=0)
191 else:
192 retval = data[pol]
193 #retval[0][self.nx-1] = -1.0
194 return retval
195
196 def __chanAverage( self ):
197 s = scantable( self.outfile, average=False )
198 nrow = s.nrow()
199 spectra = numpy.zeros( (self.npol,nrow/self.npol), dtype=float )
200 irow = 0
201 sp = [0 for i in xrange(self.nchan)]
202 for i in xrange(nrow/self.npol):
203 for ip in xrange(self.npol):
204 sp = s._getspectrum( irow )
205 spectra[ip,i] = numpy.mean( sp )
206 irow += 1
207 return spectra
208
209 def __chanIndex( self, idx ):
210 s = scantable( self.outfile, average=False )
211 nrow = s.nrow()
212 spectra = numpy.zeros( (self.npol,nrow/self.npol), dtype=float )
213 irow = 0
214 sp = [0 for i in xrange(self.nchan)]
215 for i in xrange(nrow/self.npol):
216 for ip in xrange(self.npol):
217 sp = s._getspectrum( irow )
218 spectra[ip,i] = sp[idx]
219 irow += 1
220 return spectra
221
222
Note: See TracBrowser for help on using the repository browser.