source: trunk/python/asapgrid.py@ 2849

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

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: run sdgrid task with plot=True

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Fix for plotter failure.


File size: 20.7 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 [nx,ny] = self.gridder._get_resultant_map_size()
309 [cellx,celly] = self.gridder._get_resultant_cell_size()
310 plotter = _SDGridPlotter( self.infile, self.outfile, self.ifno,
311 nx=nx, ny=ny, cellx=cellx, celly=celly )
312 plotter.plot( chan=plotchan, pol=plotpol, plotobs=plotobs, plotgrid=plotgrid )
313 # back to original setup
314 rcParams['scantable.storage'] = storg
315 t1=time.time()
316 asaplog.push('plot: elapsed time %s sec'%(t1-t0))
317 asaplog.post('DEBUG','asapgrid.plot')
318
319class asapgrid2(asapgrid_base):
320 """
321 The asapgrid class is defined to convolve data onto regular
322 spatial grid. Typical usage is as follows:
323
324 # create asapgrid instance with input scantable
325 s = scantable( 'testimage1.asap', average=False )
326 g = asapgrid( s )
327 # set IFNO if necessary
328 g.setIF( 0 )
329 # set POLNOs if necessary
330 g.setPolList( [0,1] )
331 # set SCANNOs if necessary
332 g.setScanList( [22,23,24] )
333 # define image with full specification
334 # you can skip some parameters (see help for defineImage)
335 g.defineImage( nx=12, ny=12, cellx='10arcsec', celly='10arcsec',
336 center='J2000 10h10m10s -5d05m05s' )
337 # set convolution function
338 g.setFunc( func='sf', width=3 )
339 # enable min/max clipping
340 g.enableClip()
341 # or, disable min/max clipping
342 #g.disableClip()
343 # actual gridding
344 g.grid()
345 # get result as scantable
346 sg = g.getResult()
347 """
348 def __init__( self, scantab ):
349 """
350 Create asapgrid instance.
351
352 scantab -- input data as a scantable or a list of scantables
353 to grid more than one data at once.
354 """
355 super(asapgrid2,self).__init__()
356 self.gridder = stgrid2()
357 self.scantab = scantab
358 self.setData( scantab )
359
360 def setData( self, scantab ):
361 """
362 Set data to be processed.
363
364 scantab -- input data as a scantable or a list of scantables
365 to grid more than one data at once.
366 """
367 if isinstance( scantab, scantable ):
368 self.gridder._setin( scantab )
369 else:
370 self.gridder._setfiles( scantab )
371 self.scantab = scantab
372
373 def getResult( self ):
374 """
375 Return gridded data as a scantable.
376 """
377 tp = 0 if rcParams['scantable.storage']=='memory' else 1
378 return scantable( self.gridder._get( tp ), average=False )
379
380class _SDGridPlotter:
381 def __init__( self, infile, outfile=None, ifno=-1, nx=-1, ny=-1, cellx=0.0, celly=0.0 ):
382 if isinstance( infile, str ):
383 self.infile = [infile]
384 else:
385 self.infile = infile
386 self.outfile = outfile
387 if self.outfile is None:
388 self.outfile = self.infile[0].rstrip('/')+'.grid'
389 self.nx = nx
390 self.ny = ny if ny > 0 else nx
391 self.nchan = 0
392 self.npol = 0
393 self.pollist = []
394 self.cellx = cellx
395 self.celly = celly if celly > 0.0 else cellx
396 self.center = [0.0,0.0]
397 self.nonzero = [[0.0],[0.0]]
398 self.ifno = ifno
399 self.tablein = None
400 self.nrow = 0
401 self.blc = None
402 self.trc = None
403 self.get()
404
405 def get( self ):
406 s = scantable( self.outfile, average=False )
407 self.nchan = len(s._getspectrum(0))
408 nrow = s.nrow()
409 pols = numpy.ones( nrow, dtype=int )
410 for i in xrange(nrow):
411 pols[i] = s.getpol(i)
412 self.pollist, indices = numpy.unique( pols, return_inverse=True )
413 self.npol = len(self.pollist)
414 self.pollist = self.pollist[indices[:self.npol]]
415 #print 'pollist=',self.pollist
416 #print 'npol=',self.npol
417 #print 'nrow=',nrow
418
419 if self.nx <= 0 or self.ny <= 0:
420 idx = 1
421 d0 = s.get_direction( 0 ).split()[-2]
422 d = s.get_direction(self.npol*idx)
423 while( d is not None \
424 and d.split()[-2] != d0):
425 idx += 1
426 d = s.get_direction(self.npol*idx)
427
428 self.nx = idx
429 self.ny = nrow / (self.npol * idx )
430 #print 'nx,ny=',self.nx,self.ny
431
432 self.blc = s.get_directionval( 0 )
433 self.trc = s.get_directionval( nrow-self.npol )
434 #print self.blc
435 #print self.trc
436
437 if self.cellx <= 0.0 or self.celly <= 0.0:
438 if nrow > 1:
439 incrx = s.get_directionval( self.npol )
440 incry = s.get_directionval( self.nx*self.npol )
441 else:
442 incrx = [0.0,0.0]
443 incry = [0.0,0.0]
444 self.cellx = abs( self.blc[0] - incrx[0] )
445 self.celly = abs( self.blc[1] - incry[1] )
446 #print 'cellx,celly=',self.cellx,self.celly
447
448 def plot( self, chan=-1, pol=-1, plotobs=False, plotgrid=False ):
449 if pol < 0:
450 opt = 'averaged over pol'
451 else:
452 opt = 'pol %s'%(pol)
453 if type(chan) is list:
454 opt += ', averaged over channel %s-%s'%(chan[0],chan[1])
455 elif chan < 0:
456 opt += ', averaged over channel'
457 else:
458 opt += ', channel %s'%(chan)
459 data = self.getData( chan, pol )
460 #data = numpy.fliplr( data )
461 title = 'Gridded Image (%s)'%(opt)
462 pl.figure(10)
463 pl.clf()
464 # plot grid position
465 if plotgrid:
466 x = numpy.arange(self.blc[0],self.trc[0]+0.5*self.cellx,self.cellx,dtype=float)
467 #print 'len(x)=',len(x)
468 #print 'x=',x
469 ybase = numpy.ones(self.nx,dtype=float)*self.blc[1]
470 #print 'len(ybase)=',len(ybase)
471 incr = self.celly
472 for iy in xrange(self.ny):
473 y = ybase + iy * incr
474 #print y
475 pl.plot(x,y,',',color='blue')
476 # plot observed position
477 if plotobs:
478 for i in xrange(len(self.infile)):
479 self.createTableIn( self.infile[i] )
480 irow = 0
481 while ( irow < self.nrow ):
482 chunk = self.getPointingChunk( irow )
483 #print chunk
484 pl.plot(chunk[0],chunk[1],',',color='green')
485 irow += chunk.shape[1]
486 #print irow
487 # show image
488 extent=[self.trc[0]+0.5*self.cellx,
489 self.blc[0]-0.5*self.cellx,
490 self.blc[1]-0.5*self.celly,
491 self.trc[1]+0.5*self.celly]
492 deccorr = 1.0/numpy.cos(0.5*(self.blc[1]+self.trc[1]))
493 pl.imshow(data,extent=extent,origin='lower',interpolation='nearest')
494 pl.colorbar()
495 pl.xlabel('R.A. [rad]')
496 pl.ylabel('Dec. [rad]')
497 ax = pl.axes()
498 ax.set_aspect(deccorr)
499 pl.title( title )
500
501 def createTableIn( self, tab ):
502 del self.tablein
503 self.tablein = scantable( tab, average=False )
504 if self.ifno < 0:
505 ifno = self.tablein.getif(0)
506 #print 'ifno=',ifno
507 else:
508 ifno = self.ifno
509 sel = selector()
510 sel.set_ifs( ifno )
511 self.tablein.set_selection( sel )
512 self.nchan = len(self.tablein._getspectrum(0))
513 self.nrow = self.tablein.nrow()
514 del sel
515
516
517 def getPointingChunk( self, irow ):
518 numchunk = 1000
519 nrow = min( self.nrow-irow, numchunk )
520 #print 'nrow=',nrow
521 v = numpy.zeros( (2,nrow), dtype=float )
522 idx = 0
523 for i in xrange(irow,irow+nrow):
524 d = self.tablein.get_directionval( i )
525 v[0,idx] = d[0]
526 v[1,idx] = d[1]
527 idx += 1
528 return v
529
530 def getData( self, chan=-1, pol=-1 ):
531 if type(chan) == list:
532 spectra = self.__chanAverage(start=chan[0],end=chan[1])
533 elif chan == -1:
534 spectra = self.__chanAverage()
535 else:
536 spectra = self.__chanIndex( chan )
537 data = spectra.reshape( (self.npol,self.ny,self.nx) )
538 if pol == -1:
539 retval = data.mean(axis=0)
540 else:
541 retval = data[pol]
542 return retval
543
544 def __chanAverage( self, start=-1, end=-1 ):
545 s = scantable( self.outfile, average=False )
546 nrow = s.nrow()
547 spectra = numpy.zeros( (self.npol,nrow/self.npol), dtype=float )
548 irow = 0
549 sp = [0 for i in xrange(self.nchan)]
550 if start < 0:
551 start = 0
552 if end < 0:
553 end = self.nchan
554 for i in xrange(nrow/self.npol):
555 for ip in xrange(self.npol):
556 sp = s._getspectrum( irow )[start:end]
557 spectra[ip,i] = numpy.mean( sp )
558 irow += 1
559
560 return spectra
561
562 def __chanIndex( self, idx ):
563 s = scantable( self.outfile, average=False )
564 nrow = s.nrow()
565 spectra = numpy.zeros( (self.npol,nrow/self.npol), dtype=float )
566 irow = 0
567 sp = [0 for i in xrange(self.nchan)]
568 for i in xrange(nrow/self.npol):
569 for ip in xrange(self.npol):
570 sp = s._getspectrum( irow )
571 spectra[ip,i] = sp[idx]
572 irow += 1
573 return spectra
574
575
Note: See TracBrowser for help on using the repository browser.