source: branches/hpc33/python/asapgrid.py@ 2443

Last change on this file since 2443 was 2397, 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...

Updated log messages and documentation.


File size: 13.7 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 """
11 The asapgrid class is defined to convolve data onto regular
12 spatial grid. Typical usage is as follows:
13
14 # create asapgrid instance with two input data
15 g = asapgrid( ['testimage1.asap','testimage2.asap'] )
16 # set IFNO if necessary
17 g.setIF( 0 )
18 # set POLNOs if necessary
19 g.setPolList( [0,1] )
20 # set SCANNOs if necessary
21 g.setScanList( [22,23,24] )
22 # define image with full specification
23 # you can skip some parameters (see help for defineImage)
24 g.defineImage( nx=12, ny=12, cellx='10arcsec', celly='10arcsec',
25 center='J2000 10h10m10s -5d05m05s' )
26 # set convolution function
27 g.setFunc( func='sf', width=3 )
28 # enable min/max clipping
29 g.enableClip()
30 # or, disable min/max clipping
31 #g.disableClip()
32 # actual gridding
33 g.grid()
34 # save result
35 g.save( outfile='grid.asap' )
36 # plot result
37 g.plot( plotchan=1246, plotpol=-1, plotgrid=True, plotobs=True )
38 """
39 def __init__( self, infile ):
40 """
41 Create asapgrid instance.
42
43 infile -- input data as a string or string list if you want
44 to grid more than one data at once.
45 """
46 self.outfile = None
47 self.ifno = None
48 self.gridder = stgrid()
49 self.setData( infile )
50
51 def setData( self, infile ):
52 """
53 Set data to be processed.
54
55 infile -- input data as a string or string list if you want
56 to grid more than one data at once.
57 """
58 if isinstance( infile, str ):
59 self.gridder._setin( infile )
60 else:
61 self.gridder._setfiles( infile )
62 self.infile = infile
63
64 def setIF( self, ifno ):
65 """
66 Set IFNO to be processed. Currently, asapgrid allows to process
67 only one IFNO for one gridding run even if the data contains
68 multiple IFs. If you didn't specify IFNO, default value, which
69 is IFNO in the first spectrum, will be processed.
70
71 ifno -- IFNO to be processed.
72 """
73 self.ifno = ifno
74 self.gridder._setif( self.ifno )
75
76 def setPolList( self, pollist ):
77 """
78 Set list of polarization components you want to process.
79 If not specified, all POLNOs will be processed.
80
81 pollist -- list of POLNOs.
82 """
83 self.gridder._setpollist( pollist )
84
85 def setScanList( self, scanlist ):
86 """
87 Set list of scans you want to process. If not specified, all
88 scans will be processed.
89
90 scanlist -- list of SCANNOs.
91 """
92 self.gridder._setscanlist( scanlist )
93
94 def defineImage( self, nx=-1, ny=-1, cellx='', celly='', center='' ):
95 """
96 Define spatial grid.
97
98 First two parameters, nx and ny, define number of pixels of
99 the grid. If which of those is not specified, it will be set
100 to the same value as the other. If none of them are specified,
101 it will be determined from map extent and cell size.
102
103 Next two parameters, cellx and celly, define size of pixel.
104 You should set those parameters as string, which is constructed
105 numerical value and unit, e.g. '0.5arcmin', or numerical value.
106 If those values are specified as numerical value, their units
107 will be assumed to 'arcmin'. If which of those is not specified,
108 it will be set to the same value as the other. If none of them
109 are specified, it will be determined from map extent and number
110 of pixels, or set to '1arcmin' if neither nx nor ny is set.
111
112 The last parameter, center, define the central coordinate of
113 the grid. You should specify its value as a string, like,
114
115 'J2000 05h08m50s -16d23m30s'
116
117 or
118
119 'J2000 05:08:50 -16.23.30'
120
121 You can omit equinox when you specify center coordinate. In that
122 case, J2000 is assumed. If center is not specified, it will be
123 determined from the observed positions of input data.
124
125 nx -- number of pixels along x (R.A.) direction.
126 ny -- number of pixels along y (Dec.) direction.
127 cellx -- size of pixel in x (R.A.) direction.
128 celly -- size of pixel in y (Dec.) direction.
129 center -- central position of the grid.
130 """
131 if not isinstance( cellx, str ):
132 cellx = '%sarcmin'%(cellx)
133 if not isinstance( celly, str ):
134 celly = '%sarcmin'%(celly)
135 self.gridder._defineimage( nx, ny, cellx, celly, center )
136
137 def setFunc( self, func='box', width=-1 ):
138 """
139 Set convolution function. Possible options are 'box' (Box-car,
140 default), 'sf' (prolate spheroidal), and 'gauss' (Gaussian).
141 Width of convolution function can be set using width parameter.
142 By default (-1), width is automatically set depending on each
143 convolution function. Default values for width are:
144
145 'box': 1 pixel
146 'sf': 3 pixels
147 'gauss': 3 pixels (width is used as HWHM)
148
149 func -- Function type ('box', 'sf', 'gauss').
150 width -- Width of convolution function. Default (-1) is to
151 choose pre-defined value for each convolution function.
152 """
153 self.gridder._setfunc( func, width )
154
155 def setWeight( self, weightType='uniform' ):
156 """
157 Set weight type. Possible options are 'uniform' (default),
158 'tint' (weight by integration time), 'tsys' (weight by
159 Tsys: 1/Tsys**2), and 'tintsys' (weight by integration time
160 as well as Tsys: tint/Tsys**2).
161
162 weightType -- weight type ('uniform', 'tint', 'tsys', 'tintsys')
163 """
164 self.gridder._setweight( weightType )
165
166 def enableClip( self ):
167 """
168 Enable min/max clipping.
169
170 By default, min/max clipping is disabled so that you should
171 call this method before actual gridding if you want to do
172 clipping.
173 """
174 self.gridder._enableclip()
175
176 def disableClip( self ):
177 """
178 Disable min/max clipping.
179 """
180 self.gridder._disableclip()
181
182 def grid( self ):
183 """
184 Actual gridding which will be done based on several user inputs.
185 """
186 self.gridder._grid()
187
188 def save( self, outfile='' ):
189 """
190 Save result. By default, output data name will be constructed
191 from first element of input data name list (e.g. 'input.asap.grid').
192
193 outfile -- output data name.
194 """
195 self.outfile = self.gridder._save( outfile )
196
197 def plot( self, plotchan=-1, plotpol=-1, plotobs=False, plotgrid=False ):
198 """
199 Plot gridded data.
200
201 plotchan -- Which channel you want to plot. Default (-1) is
202 to average all the channels.
203 plotpol -- Which polarization component you want to plot.
204 Default (-1) is to average all the polarization
205 components.
206 plotobs -- Also plot observed position if True. Default
207 is False. Setting True for large amount of spectra
208 might be time consuming.
209 plotgrid -- Also plot grid center if True. Default is False.
210 Setting True for large number of grids might be
211 time consuming.
212 """
213 import time
214 t0=time.time()
215 # to load scantable on disk
216 storg = rcParams['scantable.storage']
217 rcParams['scantable.storage'] = 'disk'
218 plotter = _SDGridPlotter( self.infile, self.outfile, self.ifno )
219 plotter.plot( chan=plotchan, pol=plotpol, plotobs=plotobs, plotgrid=plotgrid )
220 # back to original setup
221 rcParams['scantable.storage'] = storg
222 t1=time.time()
223 asaplog.push('plot: elapsed time %s sec'%(t1-t0))
224 asaplog.post('DEBUG','asapgrid.plot')
225
226class _SDGridPlotter:
227 def __init__( self, infile, outfile=None, ifno=-1 ):
228 if isinstance( infile, str ):
229 self.infile = [infile]
230 else:
231 self.infile = infile
232 self.outfile = outfile
233 if self.outfile is None:
234 self.outfile = self.infile[0].rstrip('/')+'.grid'
235 self.nx = -1
236 self.ny = -1
237 self.nchan = 0
238 self.npol = 0
239 self.pollist = []
240 self.cellx = 0.0
241 self.celly = 0.0
242 self.center = [0.0,0.0]
243 self.nonzero = [[0.0],[0.0]]
244 self.ifno = ifno
245 self.tablein = None
246 self.nrow = 0
247 self.blc = None
248 self.trc = None
249 self.get()
250
251 def get( self ):
252 s = scantable( self.outfile, average=False )
253 self.nchan = len(s._getspectrum(0))
254 nrow = s.nrow()
255 pols = numpy.ones( nrow, dtype=int )
256 for i in xrange(nrow):
257 pols[i] = s.getpol(i)
258 self.pollist, indices = numpy.unique( pols, return_inverse=True )
259 self.npol = len(self.pollist)
260 self.pollist = self.pollist[indices[:self.npol]]
261 #print 'pollist=',self.pollist
262 #print 'npol=',self.npol
263 #print 'nrow=',nrow
264
265 idx = 0
266 d0 = s.get_direction( 0 ).split()[-1]
267 while ( s.get_direction(self.npol*idx).split()[-1] == d0 ):
268 idx += 1
269
270 self.nx = idx
271 self.ny = nrow / (self.npol * idx )
272 #print 'nx,ny=',self.nx,self.ny
273
274 self.blc = s.get_directionval( 0 )
275 self.trc = s.get_directionval( nrow-self.npol )
276 #print self.blc
277 #print self.trc
278 incrx = s.get_directionval( self.npol )
279 incry = s.get_directionval( self.nx*self.npol )
280 self.cellx = abs( self.blc[0] - incrx[0] )
281 self.celly = abs( self.blc[1] - incry[1] )
282 #print 'cellx,celly=',self.cellx,self.celly
283
284 def plot( self, chan=-1, pol=-1, plotobs=False, plotgrid=False ):
285 if pol < 0:
286 opt = 'averaged over pol'
287 else:
288 opt = 'pol %s'%(pol)
289 if chan < 0:
290 opt += ', averaged over channel'
291 else:
292 opt += ', channel %s'%(chan)
293 data = self.getData( chan, pol )
294 title = 'Gridded Image (%s)'%(opt)
295 pl.figure(10)
296 pl.clf()
297 # plot grid position
298 if plotgrid:
299 x = numpy.arange(self.blc[0],self.trc[0]+0.5*self.cellx,self.cellx,dtype=float)
300 #print 'len(x)=',len(x)
301 #print 'x=',x
302 ybase = numpy.ones(self.nx,dtype=float)*self.blc[1]
303 #print 'len(ybase)=',len(ybase)
304 incr = self.celly
305 for iy in xrange(self.ny):
306 y = ybase + iy * incr
307 #print y
308 pl.plot(x,y,',',color='blue')
309 # plot observed position
310 if plotobs:
311 for i in xrange(len(self.infile)):
312 self.createTableIn( self.infile[i] )
313 irow = 0
314 while ( irow < self.nrow ):
315 chunk = self.getPointingChunk( irow )
316 #print chunk
317 pl.plot(chunk[0],chunk[1],',',color='green')
318 irow += chunk.shape[1]
319 #print irow
320 # show image
321 extent=[self.blc[0]-0.5*self.cellx,
322 self.trc[0]+0.5*self.cellx,
323 self.blc[1]-0.5*self.celly,
324 self.trc[1]+0.5*self.celly]
325 pl.imshow(data,extent=extent,origin='lower',interpolation='nearest')
326 pl.colorbar()
327 pl.xlabel('R.A. [rad]')
328 pl.ylabel('Dec. [rad]')
329 pl.title( title )
330
331 def createTableIn( self, tab ):
332 del self.tablein
333 self.tablein = scantable( tab, average=False )
334 if self.ifno < 0:
335 ifno = self.tablein.getif(0)
336 print 'ifno=',ifno
337 else:
338 ifno = self.ifno
339 sel = selector()
340 sel.set_ifs( ifno )
341 self.tablein.set_selection( sel )
342 self.nchan = len(self.tablein._getspectrum(0))
343 self.nrow = self.tablein.nrow()
344 del sel
345
346
347 def getPointingChunk( self, irow ):
348 numchunk = 1000
349 nrow = min( self.nrow-irow, numchunk )
350 #print 'nrow=',nrow
351 v = numpy.zeros( (2,nrow), dtype=float )
352 idx = 0
353 for i in xrange(irow,irow+nrow):
354 d = self.tablein.get_directionval( i )
355 v[0,idx] = d[0]
356 v[1,idx] = d[1]
357 idx += 1
358 return v
359
360 def getData( self, chan=-1, pol=-1 ):
361 if chan == -1:
362 spectra = self.__chanAverage()
363 else:
364 spectra = self.__chanIndex( chan )
365 data = spectra.reshape( (self.npol,self.ny,self.nx) )
366 if pol == -1:
367 retval = data.mean(axis=0)
368 else:
369 retval = data[pol]
370 return retval
371
372 def __chanAverage( self ):
373 s = scantable( self.outfile, average=False )
374 nrow = s.nrow()
375 spectra = numpy.zeros( (self.npol,nrow/self.npol), dtype=float )
376 irow = 0
377 sp = [0 for i in xrange(self.nchan)]
378 for i in xrange(nrow/self.npol):
379 for ip in xrange(self.npol):
380 sp = s._getspectrum( irow )
381 spectra[ip,i] = numpy.mean( sp )
382 irow += 1
383 return spectra
384
385 def __chanIndex( self, idx ):
386 s = scantable( self.outfile, average=False )
387 nrow = s.nrow()
388 spectra = numpy.zeros( (self.npol,nrow/self.npol), dtype=float )
389 irow = 0
390 sp = [0 for i in xrange(self.nchan)]
391 for i in xrange(nrow/self.npol):
392 for ip in xrange(self.npol):
393 sp = s._getspectrum( irow )
394 spectra[ip,i] = sp[idx]
395 irow += 1
396 return spectra
397
398
Note: See TracBrowser for help on using the repository browser.