source: trunk/python/asapgrid.py@ 2682

Last change on this file since 2682 was 2680, checked in by Takeshi Nakazato, 12 years ago

New Development: No

JIRA Issue: Yes CAS-4429

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...

Refactored asapgrid module by defining asapgrid_base class.


File size: 20.3 KB
Line 
1import numpy
2from asap import rcParams
3from asap.scantable import scantable
4from asap.selector import selector
5from asap._asap import stgrid, stgrid2
6import pylab as pl
7from logging import asaplog
8
9class asapgrid_base(object):
10 def __init__( self ):
11 self.outfile = None
12 self.ifno = None
13 self.gridder = None
14 self.infile = None
15 self.scantab = None
16
17 def setData( self, infile ):
18 raise NotImplementedError('setData is not implemented')
19
20 def setIF( self, ifno ):
21 """
22 Set IFNO to be processed. Currently, asapgrid allows to process
23 only one IFNO for one gridding run even if the data contains
24 multiple IFs. If you didn't specify IFNO, default value, which
25 is IFNO in the first spectrum, will be processed.
26
27 ifno -- IFNO to be processed.
28 """
29 self.ifno = ifno
30 self.gridder._setif( self.ifno )
31
32 def setPolList( self, pollist ):
33 """
34 Set list of polarization components you want to process.
35 If not specified, all POLNOs will be processed.
36
37 pollist -- list of POLNOs.
38 """
39 self.gridder._setpollist( pollist )
40
41 def setScanList( self, scanlist ):
42 """
43 Set list of scans you want to process. If not specified, all
44 scans will be processed.
45
46 scanlist -- list of SCANNOs.
47 """
48 self.gridder._setscanlist( scanlist )
49
50 def defineImage( self, nx=-1, ny=-1, cellx='', celly='', center='' ):
51 """
52 Define spatial grid.
53
54 First two parameters, nx and ny, define number of pixels of
55 the grid. If which of those is not specified, it will be set
56 to the same value as the other. If none of them are specified,
57 it will be determined from map extent and cell size.
58
59 Next two parameters, cellx and celly, define size of pixel.
60 You should set those parameters as string, which is constructed
61 numerical value and unit, e.g. '0.5arcmin', or numerical value.
62 If those values are specified as numerical value, their units
63 will be assumed to 'arcsec'. If which of those is not specified,
64 it will be set to the same value as the other. If none of them
65 are specified, it will be determined from map extent and number
66 of pixels, or set to '1arcmin' if neither nx nor ny is set.
67
68 The last parameter, center, define the central coordinate of
69 the grid. You should specify its value as a string, like,
70
71 'J2000 05h08m50s -16d23m30s'
72
73 or
74
75 'J2000 05:08:50 -16.23.30'
76
77 You can omit equinox when you specify center coordinate. In that
78 case, J2000 is assumed. If center is not specified, it will be
79 determined from the observed positions of input data.
80
81 nx -- number of pixels along x (R.A.) direction.
82 ny -- number of pixels along y (Dec.) direction.
83 cellx -- size of pixel in x (R.A.) direction.
84 celly -- size of pixel in y (Dec.) direction.
85 center -- central position of the grid.
86 """
87 if not isinstance( cellx, str ):
88 cellx = '%sarcsec'%(cellx)
89 if not isinstance( celly, str ):
90 celly = '%sarcsec'%(celly)
91 self.gridder._defineimage( nx, ny, cellx, celly, center )
92
93 def setFunc( self, func='box', convsupport=-1, truncate="-1", gwidth="-1", jwidth="-1" ):
94 """
95 Set convolution function. Possible options are 'box' (Box-car,
96 default), 'sf' (prolate spheroidal), 'gauss' (Gaussian), and
97 'gjinc' (Gaussian * Jinc).
98 Width of convolution function can be set using several parameters.
99 For 'box' and 'sf', we have one parameter, convsupport, that
100 specifies a cut-off radius of the convlolution function. By default
101 (-1), convsupport is automatically set depending on each convolution
102 function. Default values for convsupport are:
103
104 'box': 1 pixel
105 'sf': 3 pixels
106
107 For 'gauss', we have two parameters for convolution function,
108 truncate and gwidth. The truncate is similar to convsupport
109 except that truncate allows to specify its value as float or
110 string consisting of numeric and unit (e.g. '10arcsec' or
111 '3pixel'). Available units are angular units ('arcsec', 'arcmin',
112 'deg', etc.) and 'pixel'. Default unit is 'pixel' so that if
113 you specify numerical value or string without unit to gwidth,
114 the value will be interpreted as 'pixel'. gwidth is an HWHM of
115 gaussian. It also allows string value. Interpretation of the
116 value for gwidth is same as truncate. Default value for 'gauss'
117 is
118
119 gwidth: '-1' ---> sqrt(log(2.0)) pixel
120 truncate: '-1' ---> 3*gwidth pixel
121
122 For 'gjinc', there is an additional parameter jwidth that
123 specifies a width of the jinc function whose functional form is
124
125 jinc(x) = J_1(pi*x/jwidth) / (pi*x/jwidth)
126
127 Default values for 'gjinc' is
128
129 gwidth: '-1' ---> 2.52*sqrt(log(2.0)) pixel
130 jwidth: '-1' ---> 1.55
131 truncate: '-1' ---> automatically truncate at first null
132
133 Default values for gwidth and jwidth are taken from Mangum et al.
134 (2007).
135
136 func -- Function type ('box', 'sf', 'gauss', 'gjinc').
137 convsupport -- Width of convolution function. Default (-1) is
138 to choose pre-defined value for each convolution
139 function. Effective only for 'box' and 'sf'.
140 truncate -- Truncation radius of the convolution function.
141 Acceptable value is an integer or a float in units of
142 pixel, or a string consisting of numeric plus unit.
143 Default unit for the string is 'pixel'. Default (-1)
144 is to choose pre-defined value for each convolution
145 function. Effective only for 'gauss' and 'gjinc'.
146 gwidth -- The HWHM of the gaussian. Acceptable value is an integer
147 or a float in units of pixel, or a string consisting of
148 numeric plus unit. Default unit for the string is 'pixel'.
149 Default (-1) is to choose pre-defined value for each
150 convolution function. Effective only for 'gauss' and
151 'gjinc'.
152 jwidth -- The width of the jinc function. Acceptable value is an
153 integer or a float in units of pixel, or a string
154 consisting of numeric plus unit. Default unit for the
155 string is 'pixel'. Default (-1) is to choose pre-defined
156 value for each convolution function. Effective only for
157 'gjinc'.
158 """
159 self.gridder._setfunc(func,
160 convsupport=convsupport,
161 truncate=truncate,
162 gwidth=gwidth,
163 jwidth=jwidth)
164
165 def setWeight( self, weightType='uniform' ):
166 """
167 Set weight type. Possible options are 'uniform' (default),
168 'tint' (weight by integration time), 'tsys' (weight by
169 Tsys: 1/Tsys**2), and 'tintsys' (weight by integration time
170 as well as Tsys: tint/Tsys**2).
171
172 weightType -- weight type ('uniform', 'tint', 'tsys', 'tintsys')
173 """
174 self.gridder._setweight( weightType )
175
176 def enableClip( self ):
177 """
178 Enable min/max clipping.
179
180 By default, min/max clipping is disabled so that you should
181 call this method before actual gridding if you want to do
182 clipping.
183 """
184 self.gridder._enableclip()
185
186 def disableClip( self ):
187 """
188 Disable min/max clipping.
189 """
190 self.gridder._disableclip()
191
192 def grid( self ):
193 """
194 Actual gridding which will be done based on several user inputs.
195 """
196 self.gridder._grid()
197
198 def plotFunc(self, clear=True):
199 """
200 Support function to see the shape of current grid function.
201
202 clear -- clear panel if True. Default is True.
203 """
204 pl.figure(11)
205 if clear:
206 pl.clf()
207 f = self.gridder._getfunc()
208 convsampling = 100
209 a = numpy.arange(0,len(f)/convsampling,1./convsampling,dtype=float)
210 pl.plot(a,f,'.-')
211 pl.xlabel('pixel')
212 pl.ylabel('convFunc')
213
214 def save( self, outfile='' ):
215 raise NotImplementedError('save is not implemented')
216
217 def plot( self, plotchan=-1, plotpol=-1, plotobs=False, plotgrid=False ):
218 raise NotImplementedError('plot is not implemented')
219
220 def getResult( self ):
221 raise NotImplementedError('getResult is not implemented')
222
223class asapgrid(asapgrid_base):
224 """
225 The asapgrid class is defined to convolve data onto regular
226 spatial grid. Typical usage is as follows:
227
228 # create asapgrid instance with two input data
229 g = asapgrid( ['testimage1.asap','testimage2.asap'] )
230 # set IFNO if necessary
231 g.setIF( 0 )
232 # set POLNOs if necessary
233 g.setPolList( [0,1] )
234 # set SCANNOs if necessary
235 g.setScanList( [22,23,24] )
236 # define image with full specification
237 # you can skip some parameters (see help for defineImage)
238 g.defineImage( nx=12, ny=12, cellx='10arcsec', celly='10arcsec',
239 center='J2000 10h10m10s -5d05m05s' )
240 # set convolution function
241 g.setFunc( func='sf', convsupport=3 )
242 # enable min/max clipping
243 g.enableClip()
244 # or, disable min/max clipping
245 #g.disableClip()
246 # actual gridding
247 g.grid()
248 # save result
249 g.save( outfile='grid.asap' )
250 # plot result
251 g.plot( plotchan=1246, plotpol=-1, plotgrid=True, plotobs=True )
252 """
253 def __init__( self, infile ):
254 """
255 Create asapgrid instance.
256
257 infile -- input data as a string or string list if you want
258 to grid more than one data at once.
259 """
260 super(asapgrid,self).__init__()
261 self.gridder = stgrid()
262 self.infile=infile
263 self.setData(infile)
264
265 def setData( self, infile ):
266 """
267 Set data to be processed.
268
269 infile -- input data as a string or string list if you want
270 to grid more than one data at once.
271 """
272 if isinstance( infile, str ):
273 self.gridder._setin( infile )
274 else:
275 self.gridder._setfiles( infile )
276 self.infile = infile
277
278 def save( self, outfile='' ):
279 """
280 Save result. By default, output data name will be constructed
281 from first element of input data name list (e.g. 'input.asap.grid').
282
283 outfile -- output data name.
284 """
285 self.outfile = self.gridder._save( outfile )
286
287 def plot( self, plotchan=-1, plotpol=-1, plotobs=False, plotgrid=False ):
288 """
289 Plot gridded data.
290
291 plotchan -- Which channel you want to plot. Default (-1) is
292 to average all the channels.
293 plotpol -- Which polarization component you want to plot.
294 Default (-1) is to average all the polarization
295 components.
296 plotobs -- Also plot observed position if True. Default
297 is False. Setting True for large amount of spectra
298 might be time consuming.
299 plotgrid -- Also plot grid center if True. Default is False.
300 Setting True for large number of grids might be
301 time consuming.
302 """
303 import time
304 t0=time.time()
305 # to load scantable on disk
306 storg = rcParams['scantable.storage']
307 rcParams['scantable.storage'] = 'disk'
308 plotter = _SDGridPlotter( self.infile, self.outfile, self.ifno )
309 plotter.plot( chan=plotchan, pol=plotpol, plotobs=plotobs, plotgrid=plotgrid )
310 # back to original setup
311 rcParams['scantable.storage'] = storg
312 t1=time.time()
313 asaplog.push('plot: elapsed time %s sec'%(t1-t0))
314 asaplog.post('DEBUG','asapgrid.plot')
315
316class asapgrid2(asapgrid_base):
317 """
318 The asapgrid class is defined to convolve data onto regular
319 spatial grid. Typical usage is as follows:
320
321 # create asapgrid instance with input scantable
322 s = scantable( 'testimage1.asap', average=False )
323 g = asapgrid( s )
324 # set IFNO if necessary
325 g.setIF( 0 )
326 # set POLNOs if necessary
327 g.setPolList( [0,1] )
328 # set SCANNOs if necessary
329 g.setScanList( [22,23,24] )
330 # define image with full specification
331 # you can skip some parameters (see help for defineImage)
332 g.defineImage( nx=12, ny=12, cellx='10arcsec', celly='10arcsec',
333 center='J2000 10h10m10s -5d05m05s' )
334 # set convolution function
335 g.setFunc( func='sf', width=3 )
336 # enable min/max clipping
337 g.enableClip()
338 # or, disable min/max clipping
339 #g.disableClip()
340 # actual gridding
341 g.grid()
342 # get result as scantable
343 sg = g.getResult()
344 """
345 def __init__( self, scantab ):
346 """
347 Create asapgrid instance.
348
349 scantab -- input data as a scantable or a list of scantables
350 to grid more than one data at once.
351 """
352 super(asapgrid2,self).__init__()
353 self.gridder = stgrid2()
354 self.scantab = scantab
355 self.setData( scantab )
356
357 def setData( self, scantab ):
358 """
359 Set data to be processed.
360
361 scantab -- input data as a scantable or a list of scantables
362 to grid more than one data at once.
363 """
364 if isinstance( scantab, scantable ):
365 self.gridder._setin( scantab )
366 else:
367 self.gridder._setfiles( scantab )
368 self.scantab = scantab
369
370 def getResult( self ):
371 """
372 Return gridded data as a scantable.
373 """
374 tp = 0 if rcParams['scantable.storage']=='memory' else 1
375 return scantable( self.gridder._get( tp ), average=False )
376
377class _SDGridPlotter:
378 def __init__( self, infile, outfile=None, ifno=-1 ):
379 if isinstance( infile, str ):
380 self.infile = [infile]
381 else:
382 self.infile = infile
383 self.outfile = outfile
384 if self.outfile is None:
385 self.outfile = self.infile[0].rstrip('/')+'.grid'
386 self.nx = -1
387 self.ny = -1
388 self.nchan = 0
389 self.npol = 0
390 self.pollist = []
391 self.cellx = 0.0
392 self.celly = 0.0
393 self.center = [0.0,0.0]
394 self.nonzero = [[0.0],[0.0]]
395 self.ifno = ifno
396 self.tablein = None
397 self.nrow = 0
398 self.blc = None
399 self.trc = None
400 self.get()
401
402 def get( self ):
403 s = scantable( self.outfile, average=False )
404 self.nchan = len(s._getspectrum(0))
405 nrow = s.nrow()
406 pols = numpy.ones( nrow, dtype=int )
407 for i in xrange(nrow):
408 pols[i] = s.getpol(i)
409 self.pollist, indices = numpy.unique( pols, return_inverse=True )
410 self.npol = len(self.pollist)
411 self.pollist = self.pollist[indices[:self.npol]]
412 #print 'pollist=',self.pollist
413 #print 'npol=',self.npol
414 #print 'nrow=',nrow
415
416 idx = 1
417 d0 = s.get_direction( 0 ).split()[-2]
418 d = s.get_direction(self.npol*idx)
419 while( d is not None \
420 and d.split()[-2] != d0):
421 idx += 1
422 d = s.get_direction(self.npol*idx)
423
424 self.nx = idx
425 self.ny = nrow / (self.npol * idx )
426 #print 'nx,ny=',self.nx,self.ny
427
428 self.blc = s.get_directionval( 0 )
429 self.trc = s.get_directionval( nrow-self.npol )
430 #print self.blc
431 #print self.trc
432 if nrow > 1:
433 incrx = s.get_directionval( self.npol )
434 incry = s.get_directionval( self.nx*self.npol )
435 else:
436 incrx = [0.0,0.0]
437 incry = [0.0,0.0]
438 self.cellx = abs( self.blc[0] - incrx[0] )
439 self.celly = abs( self.blc[1] - incry[1] )
440 #print 'cellx,celly=',self.cellx,self.celly
441
442 def plot( self, chan=-1, pol=-1, plotobs=False, plotgrid=False ):
443 if pol < 0:
444 opt = 'averaged over pol'
445 else:
446 opt = 'pol %s'%(pol)
447 if type(chan) is list:
448 opt += ', averaged over channel %s-%s'%(chan[0],chan[1])
449 elif chan < 0:
450 opt += ', averaged over channel'
451 else:
452 opt += ', channel %s'%(chan)
453 data = self.getData( chan, pol )
454 #data = numpy.fliplr( data )
455 title = 'Gridded Image (%s)'%(opt)
456 pl.figure(10)
457 pl.clf()
458 # plot grid position
459 if plotgrid:
460 x = numpy.arange(self.blc[0],self.trc[0]+0.5*self.cellx,self.cellx,dtype=float)
461 #print 'len(x)=',len(x)
462 #print 'x=',x
463 ybase = numpy.ones(self.nx,dtype=float)*self.blc[1]
464 #print 'len(ybase)=',len(ybase)
465 incr = self.celly
466 for iy in xrange(self.ny):
467 y = ybase + iy * incr
468 #print y
469 pl.plot(x,y,',',color='blue')
470 # plot observed position
471 if plotobs:
472 for i in xrange(len(self.infile)):
473 self.createTableIn( self.infile[i] )
474 irow = 0
475 while ( irow < self.nrow ):
476 chunk = self.getPointingChunk( irow )
477 #print chunk
478 pl.plot(chunk[0],chunk[1],',',color='green')
479 irow += chunk.shape[1]
480 #print irow
481 # show image
482 extent=[self.trc[0]+0.5*self.cellx,
483 self.blc[0]-0.5*self.cellx,
484 self.blc[1]-0.5*self.celly,
485 self.trc[1]+0.5*self.celly]
486 deccorr = 1.0/numpy.cos(0.5*(self.blc[1]+self.trc[1]))
487 pl.imshow(data,extent=extent,origin='lower',interpolation='nearest')
488 pl.colorbar()
489 pl.xlabel('R.A. [rad]')
490 pl.ylabel('Dec. [rad]')
491 ax = pl.axes()
492 ax.set_aspect(deccorr)
493 pl.title( title )
494
495 def createTableIn( self, tab ):
496 del self.tablein
497 self.tablein = scantable( tab, average=False )
498 if self.ifno < 0:
499 ifno = self.tablein.getif(0)
500 #print 'ifno=',ifno
501 else:
502 ifno = self.ifno
503 sel = selector()
504 sel.set_ifs( ifno )
505 self.tablein.set_selection( sel )
506 self.nchan = len(self.tablein._getspectrum(0))
507 self.nrow = self.tablein.nrow()
508 del sel
509
510
511 def getPointingChunk( self, irow ):
512 numchunk = 1000
513 nrow = min( self.nrow-irow, numchunk )
514 #print 'nrow=',nrow
515 v = numpy.zeros( (2,nrow), dtype=float )
516 idx = 0
517 for i in xrange(irow,irow+nrow):
518 d = self.tablein.get_directionval( i )
519 v[0,idx] = d[0]
520 v[1,idx] = d[1]
521 idx += 1
522 return v
523
524 def getData( self, chan=-1, pol=-1 ):
525 if type(chan) == list:
526 spectra = self.__chanAverage(start=chan[0],end=chan[1])
527 elif chan == -1:
528 spectra = self.__chanAverage()
529 else:
530 spectra = self.__chanIndex( chan )
531 data = spectra.reshape( (self.npol,self.ny,self.nx) )
532 if pol == -1:
533 retval = data.mean(axis=0)
534 else:
535 retval = data[pol]
536 return retval
537
538 def __chanAverage( self, start=-1, end=-1 ):
539 s = scantable( self.outfile, average=False )
540 nrow = s.nrow()
541 spectra = numpy.zeros( (self.npol,nrow/self.npol), dtype=float )
542 irow = 0
543 sp = [0 for i in xrange(self.nchan)]
544 if start < 0:
545 start = 0
546 if end < 0:
547 end = self.nchan
548 for i in xrange(nrow/self.npol):
549 for ip in xrange(self.npol):
550 sp = s._getspectrum( irow )[start:end]
551 spectra[ip,i] = numpy.mean( sp )
552 irow += 1
553
554 return spectra
555
556 def __chanIndex( self, idx ):
557 s = scantable( self.outfile, average=False )
558 nrow = s.nrow()
559 spectra = numpy.zeros( (self.npol,nrow/self.npol), dtype=float )
560 irow = 0
561 sp = [0 for i in xrange(self.nchan)]
562 for i in xrange(nrow/self.npol):
563 for ip in xrange(self.npol):
564 sp = s._getspectrum( irow )
565 spectra[ip,i] = sp[idx]
566 irow += 1
567 return spectra
568
569
Note: See TracBrowser for help on using the repository browser.