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

Last change on this file since 2526 was 2459, checked in by Kana Sugimoto, 13 years ago

merged Takeshi's fix to asapgrid and STGrid class form trunk

File size: 14.3 KB
RevLine 
[2356]1import numpy
[2367]2from asap import rcParams
[2356]3from asap.scantable import scantable
[2367]4from asap.selector import selector
[2356]5from asap._asap import stgrid
6import pylab as pl
[2367]7from logging import asaplog
[2356]8
9class asapgrid:
[2391]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 )
[2397]28 # enable min/max clipping
29 g.enableClip()
30 # or, disable min/max clipping
31 #g.disableClip()
[2391]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 """
[2356]39 def __init__( self, infile ):
[2391]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 """
[2356]46 self.outfile = None
[2367]47 self.ifno = None
[2390]48 self.gridder = stgrid()
49 self.setData( infile )
[2356]50
51 def setData( self, infile ):
[2391]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 """
[2390]58 if isinstance( infile, str ):
59 self.gridder._setin( infile )
60 else:
61 self.gridder._setfiles( infile )
62 self.infile = infile
[2356]63
[2362]64 def setIF( self, ifno ):
[2391]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 """
[2367]73 self.ifno = ifno
74 self.gridder._setif( self.ifno )
[2362]75
[2360]76 def setPolList( self, pollist ):
[2391]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 """
[2360]83 self.gridder._setpollist( pollist )
84
[2364]85 def setScanList( self, scanlist ):
[2391]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 """
[2364]92 self.gridder._setscanlist( scanlist )
93
[2356]94 def defineImage( self, nx=-1, ny=-1, cellx='', celly='', center='' ):
[2391]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
[2459]107 will be assumed to 'arcsec'. If which of those is not specified,
[2391]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 ):
[2459]132 cellx = '%sarcsec'%(cellx)
[2391]133 if not isinstance( celly, str ):
[2459]134 celly = '%sarcsec'%(celly)
[2356]135 self.gridder._defineimage( nx, ny, cellx, celly, center )
136
[2364]137 def setFunc( self, func='box', width=-1 ):
[2391]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
[2459]147 'gauss': 1 pixel (width is used as HWHM)
[2391]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 """
[2364]153 self.gridder._setfunc( func, width )
[2356]154
[2361]155 def setWeight( self, weightType='uniform' ):
[2391]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 """
[2396]164 self.gridder._setweight( weightType )
[2361]165
[2396]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
[2356]182 def grid( self ):
[2391]183 """
184 Actual gridding which will be done based on several user inputs.
185 """
[2356]186 self.gridder._grid()
187
188 def save( self, outfile='' ):
[2391]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 """
[2356]195 self.outfile = self.gridder._save( outfile )
196
[2375]197 def plot( self, plotchan=-1, plotpol=-1, plotobs=False, plotgrid=False ):
[2391]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 """
[2367]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 )
[2375]219 plotter.plot( chan=plotchan, pol=plotpol, plotobs=plotobs, plotgrid=plotgrid )
[2367]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')
[2356]225
226class _SDGridPlotter:
[2373]227 def __init__( self, infile, outfile=None, ifno=-1 ):
[2390]228 if isinstance( infile, str ):
229 self.infile = [infile]
230 else:
231 self.infile = infile
[2356]232 self.outfile = outfile
233 if self.outfile is None:
[2390]234 self.outfile = self.infile[0].rstrip('/')+'.grid'
[2356]235 self.nx = -1
236 self.ny = -1
237 self.nchan = 0
[2360]238 self.npol = 0
239 self.pollist = []
[2356]240 self.cellx = 0.0
241 self.celly = 0.0
242 self.center = [0.0,0.0]
243 self.nonzero = [[0.0],[0.0]]
[2367]244 self.ifno = ifno
[2372]245 self.tablein = None
246 self.nrow = 0
247 self.blc = None
248 self.trc = None
[2356]249 self.get()
250
251 def get( self ):
252 s = scantable( self.outfile, average=False )
[2387]253 self.nchan = len(s._getspectrum(0))
[2356]254 nrow = s.nrow()
[2360]255 pols = numpy.ones( nrow, dtype=int )
[2356]256 for i in xrange(nrow):
[2360]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
[2356]264
265 idx = 0
[2372]266 d0 = s.get_direction( 0 ).split()[-1]
[2459]267 while ( s.get_direction(self.npol*idx) is not None \
268 and s.get_direction(self.npol*idx).split()[-1] == d0 ):
[2356]269 idx += 1
[2367]270
[2372]271 self.nx = idx
272 self.ny = nrow / (self.npol * idx )
[2360]273 #print 'nx,ny=',self.nx,self.ny
[2372]274
275 self.blc = s.get_directionval( 0 )
276 self.trc = s.get_directionval( nrow-self.npol )
277 #print self.blc
278 #print self.trc
[2459]279 if nrow > 1:
280 incrx = s.get_directionval( self.npol )
281 incry = s.get_directionval( self.nx*self.npol )
282 else:
283 incrx = [0.0,0.0]
284 incry = [0.0,0.0]
[2372]285 self.cellx = abs( self.blc[0] - incrx[0] )
286 self.celly = abs( self.blc[1] - incry[1] )
[2360]287 #print 'cellx,celly=',self.cellx,self.celly
[2356]288
[2375]289 def plot( self, chan=-1, pol=-1, plotobs=False, plotgrid=False ):
[2360]290 if pol < 0:
291 opt = 'averaged over pol'
[2356]292 else:
[2360]293 opt = 'pol %s'%(pol)
[2459]294 if type(chan) is list:
295 opt += ', averaged over channel %s-%s'%(chan[0],chan[1])
296 elif chan < 0:
[2360]297 opt += ', averaged over channel'
298 else:
299 opt += ', channel %s'%(chan)
[2459]300 data = self.getData( chan, pol )
301 data = numpy.fliplr( data )
[2360]302 title = 'Gridded Image (%s)'%(opt)
[2356]303 pl.figure(10)
304 pl.clf()
[2372]305 # plot grid position
[2375]306 if plotgrid:
307 x = numpy.arange(self.blc[0],self.trc[0]+0.5*self.cellx,self.cellx,dtype=float)
308 #print 'len(x)=',len(x)
309 #print 'x=',x
310 ybase = numpy.ones(self.nx,dtype=float)*self.blc[1]
311 #print 'len(ybase)=',len(ybase)
312 incr = self.celly
313 for iy in xrange(self.ny):
314 y = ybase + iy * incr
315 #print y
316 pl.plot(x,y,',',color='blue')
[2372]317 # plot observed position
[2375]318 if plotobs:
[2390]319 for i in xrange(len(self.infile)):
320 self.createTableIn( self.infile[i] )
321 irow = 0
322 while ( irow < self.nrow ):
323 chunk = self.getPointingChunk( irow )
324 #print chunk
325 pl.plot(chunk[0],chunk[1],',',color='green')
326 irow += chunk.shape[1]
327 #print irow
[2372]328 # show image
[2459]329 extent=[self.trc[0]+0.5*self.cellx,
330 self.blc[0]-0.5*self.cellx,
[2372]331 self.blc[1]-0.5*self.celly,
332 self.trc[1]+0.5*self.celly]
[2459]333 deccorr = 1.0/numpy.cos(0.5*(self.blc[1]+self.trc[1]))
[2356]334 pl.imshow(data,extent=extent,origin='lower',interpolation='nearest')
335 pl.colorbar()
336 pl.xlabel('R.A. [rad]')
337 pl.ylabel('Dec. [rad]')
[2459]338 ax = pl.axes()
339 ax.set_aspect(deccorr)
[2358]340 pl.title( title )
[2367]341
[2390]342 def createTableIn( self, tab ):
343 del self.tablein
344 self.tablein = scantable( tab, average=False )
[2387]345 if self.ifno < 0:
346 ifno = self.tablein.getif(0)
347 print 'ifno=',ifno
348 else:
349 ifno = self.ifno
350 sel = selector()
351 sel.set_ifs( ifno )
[2390]352 self.tablein.set_selection( sel )
[2387]353 self.nchan = len(self.tablein._getspectrum(0))
[2390]354 self.nrow = self.tablein.nrow()
[2387]355 del sel
356
357
[2372]358 def getPointingChunk( self, irow ):
359 numchunk = 1000
360 nrow = min( self.nrow-irow, numchunk )
361 #print 'nrow=',nrow
362 v = numpy.zeros( (2,nrow), dtype=float )
363 idx = 0
364 for i in xrange(irow,irow+nrow):
365 d = self.tablein.get_directionval( i )
366 v[0,idx] = d[0]
367 v[1,idx] = d[1]
368 idx += 1
369 return v
370
[2367]371 def getData( self, chan=-1, pol=-1 ):
[2459]372 if type(chan) == list:
373 spectra = self.__chanAverage(start=chan[0],end=chan[1])
374 elif chan == -1:
[2367]375 spectra = self.__chanAverage()
376 else:
377 spectra = self.__chanIndex( chan )
[2372]378 data = spectra.reshape( (self.npol,self.ny,self.nx) )
[2367]379 if pol == -1:
380 retval = data.mean(axis=0)
381 else:
382 retval = data[pol]
383 return retval
384
[2459]385 def __chanAverage( self, start=-1, end=-1 ):
[2367]386 s = scantable( self.outfile, average=False )
[2372]387 nrow = s.nrow()
[2367]388 spectra = numpy.zeros( (self.npol,nrow/self.npol), dtype=float )
389 irow = 0
390 sp = [0 for i in xrange(self.nchan)]
[2459]391 if start < 0:
392 start = 0
393 if end < 0:
394 end = self.nchan
[2367]395 for i in xrange(nrow/self.npol):
396 for ip in xrange(self.npol):
[2459]397 sp = s._getspectrum( irow )[start:end]
[2367]398 spectra[ip,i] = numpy.mean( sp )
399 irow += 1
[2459]400
[2367]401 return spectra
402
403 def __chanIndex( self, idx ):
404 s = scantable( self.outfile, average=False )
[2372]405 nrow = s.nrow()
[2367]406 spectra = numpy.zeros( (self.npol,nrow/self.npol), dtype=float )
407 irow = 0
408 sp = [0 for i in xrange(self.nchan)]
409 for i in xrange(nrow/self.npol):
410 for ip in xrange(self.npol):
411 sp = s._getspectrum( irow )
412 spectra[ip,i] = sp[idx]
413 irow += 1
414 return spectra
415
416
Note: See TracBrowser for help on using the repository browser.