source: trunk/python/asapgrid.py@ 2676

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

Implemented gauss and gjinc gridding. Added some arguments to
pass necessary parameters for those grid functions.

Also, defined helper function that plots current grid function
against radius in pixel.


File size: 22.4 KB
RevLine 
[2356]1import numpy
[2367]2from asap import rcParams
[2356]3from asap.scantable import scantable
[2367]4from asap.selector import selector
[2593]5from asap._asap import stgrid, stgrid2
[2356]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
[2421]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 ):
[2421]132 cellx = '%sarcsec'%(cellx)
[2391]133 if not isinstance( celly, str ):
[2421]134 celly = '%sarcsec'%(celly)
[2356]135 self.gridder._defineimage( nx, ny, cellx, celly, center )
136
[2671]137 def setFunc( self, func='box', width=-1, gwidth="", jwidth="" ):
[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
[2450]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 """
[2671]153 fname = func.upper()
154 if fname == 'GAUSS' or fname == 'GJINC':
155 gw = str(gwidth)
156 jw = str(jwidth)
157 w = str(width)
158 if w[0] == '-': w = ''
159 #self.gridder._setfunc(fname, -1, w, gw, jw)
160 self.gridder._setfunc(fname,convsupport=-1,gwidth=gw,jwidth=jw,truncate=w)
161 else:
162 self.gridder._setfunc( func, convsupport=width )
[2356]163
[2361]164 def setWeight( self, weightType='uniform' ):
[2391]165 """
166 Set weight type. Possible options are 'uniform' (default),
167 'tint' (weight by integration time), 'tsys' (weight by
168 Tsys: 1/Tsys**2), and 'tintsys' (weight by integration time
169 as well as Tsys: tint/Tsys**2).
170
171 weightType -- weight type ('uniform', 'tint', 'tsys', 'tintsys')
172 """
[2396]173 self.gridder._setweight( weightType )
[2361]174
[2396]175 def enableClip( self ):
176 """
177 Enable min/max clipping.
178
179 By default, min/max clipping is disabled so that you should
180 call this method before actual gridding if you want to do
181 clipping.
182 """
183 self.gridder._enableclip()
184
185 def disableClip( self ):
186 """
187 Disable min/max clipping.
188 """
189 self.gridder._disableclip()
190
[2356]191 def grid( self ):
[2391]192 """
193 Actual gridding which will be done based on several user inputs.
194 """
[2356]195 self.gridder._grid()
196
197 def save( self, outfile='' ):
[2391]198 """
199 Save result. By default, output data name will be constructed
200 from first element of input data name list (e.g. 'input.asap.grid').
201
202 outfile -- output data name.
203 """
[2356]204 self.outfile = self.gridder._save( outfile )
205
[2375]206 def plot( self, plotchan=-1, plotpol=-1, plotobs=False, plotgrid=False ):
[2391]207 """
208 Plot gridded data.
209
210 plotchan -- Which channel you want to plot. Default (-1) is
211 to average all the channels.
212 plotpol -- Which polarization component you want to plot.
213 Default (-1) is to average all the polarization
214 components.
215 plotobs -- Also plot observed position if True. Default
216 is False. Setting True for large amount of spectra
217 might be time consuming.
218 plotgrid -- Also plot grid center if True. Default is False.
219 Setting True for large number of grids might be
220 time consuming.
221 """
[2367]222 import time
223 t0=time.time()
224 # to load scantable on disk
225 storg = rcParams['scantable.storage']
226 rcParams['scantable.storage'] = 'disk'
227 plotter = _SDGridPlotter( self.infile, self.outfile, self.ifno )
[2375]228 plotter.plot( chan=plotchan, pol=plotpol, plotobs=plotobs, plotgrid=plotgrid )
[2367]229 # back to original setup
230 rcParams['scantable.storage'] = storg
231 t1=time.time()
232 asaplog.push('plot: elapsed time %s sec'%(t1-t0))
233 asaplog.post('DEBUG','asapgrid.plot')
[2671]234
235 def plotFunc(self, clear=True):
236 """
237 Support function to see the shape of current grid function.
238
239 clear -- clear panel if True. Default is True.
240 """
241 pl.figure(11)
242 if clear:
243 pl.clf()
244 f = self.gridder._getfunc()
245 convsampling = 100
246 a = numpy.arange(0,len(f)/convsampling,1./convsampling,dtype=float)
247 pl.plot(a,f,'.-')
248 pl.xlabel('pixel')
249 pl.ylabel('convFunc')
[2356]250
[2593]251class asapgrid2:
252 """
253 The asapgrid class is defined to convolve data onto regular
254 spatial grid. Typical usage is as follows:
255
256 # create asapgrid instance with input scantable
257 s = scantable( 'testimage1.asap', average=False )
258 g = asapgrid( s )
259 # set IFNO if necessary
260 g.setIF( 0 )
261 # set POLNOs if necessary
262 g.setPolList( [0,1] )
263 # set SCANNOs if necessary
264 g.setScanList( [22,23,24] )
265 # define image with full specification
266 # you can skip some parameters (see help for defineImage)
267 g.defineImage( nx=12, ny=12, cellx='10arcsec', celly='10arcsec',
268 center='J2000 10h10m10s -5d05m05s' )
269 # set convolution function
270 g.setFunc( func='sf', width=3 )
271 # enable min/max clipping
272 g.enableClip()
273 # or, disable min/max clipping
274 #g.disableClip()
275 # actual gridding
276 g.grid()
277 # get result as scantable
278 sg = g.getResult()
279 """
280 def __init__( self, scantab ):
281 """
282 Create asapgrid instance.
283
284 scantab -- input data as a scantable or a list of scantables
285 to grid more than one data at once.
286 """
287 self.outfile = None
288 self.ifno = None
289 self.gridder = stgrid2()
290 self.setData( scantab )
291
292 def setData( self, scantab ):
293 """
294 Set data to be processed.
295
296 scantab -- input data as a scantable or a list of scantables
297 to grid more than one data at once.
298 """
299 if isinstance( scantab, scantable ):
300 self.gridder._setin( scantab )
301 else:
302 self.gridder._setfiles( scantab )
303 self.scantab = scantab
304
305 def setIF( self, ifno ):
306 """
307 Set IFNO to be processed. Currently, asapgrid allows to process
308 only one IFNO for one gridding run even if the data contains
309 multiple IFs. If you didn't specify IFNO, default value, which
310 is IFNO in the first spectrum, will be processed.
311
312 ifno -- IFNO to be processed.
313 """
314 self.ifno = ifno
315 self.gridder._setif( self.ifno )
316
317 def setPolList( self, pollist ):
318 """
319 Set list of polarization components you want to process.
320 If not specified, all POLNOs will be processed.
321
322 pollist -- list of POLNOs.
323 """
324 self.gridder._setpollist( pollist )
325
326 def setScanList( self, scanlist ):
327 """
328 Set list of scans you want to process. If not specified, all
329 scans will be processed.
330
331 scanlist -- list of SCANNOs.
332 """
333 self.gridder._setscanlist( scanlist )
334
335 def defineImage( self, nx=-1, ny=-1, cellx='', celly='', center='' ):
336 """
337 Define spatial grid.
338
339 First two parameters, nx and ny, define number of pixels of
340 the grid. If which of those is not specified, it will be set
341 to the same value as the other. If none of them are specified,
342 it will be determined from map extent and cell size.
343
344 Next two parameters, cellx and celly, define size of pixel.
345 You should set those parameters as string, which is constructed
346 numerical value and unit, e.g. '0.5arcmin', or numerical value.
347 If those values are specified as numerical value, their units
348 will be assumed to 'arcsec'. If which of those is not specified,
349 it will be set to the same value as the other. If none of them
350 are specified, it will be determined from map extent and number
351 of pixels, or set to '1arcmin' if neither nx nor ny is set.
352
353 The last parameter, center, define the central coordinate of
354 the grid. You should specify its value as a string, like,
355
356 'J2000 05h08m50s -16d23m30s'
357
358 or
359
360 'J2000 05:08:50 -16.23.30'
361
362 You can omit equinox when you specify center coordinate. In that
363 case, J2000 is assumed. If center is not specified, it will be
364 determined from the observed positions of input data.
365
366 nx -- number of pixels along x (R.A.) direction.
367 ny -- number of pixels along y (Dec.) direction.
368 cellx -- size of pixel in x (R.A.) direction.
369 celly -- size of pixel in y (Dec.) direction.
370 center -- central position of the grid.
371 """
372 if not isinstance( cellx, str ):
373 cellx = '%sarcsec'%(cellx)
374 if not isinstance( celly, str ):
375 celly = '%sarcsec'%(celly)
376 self.gridder._defineimage( nx, ny, cellx, celly, center )
377
378 def setFunc( self, func='box', width=-1 ):
379 """
380 Set convolution function. Possible options are 'box' (Box-car,
381 default), 'sf' (prolate spheroidal), and 'gauss' (Gaussian).
382 Width of convolution function can be set using width parameter.
383 By default (-1), width is automatically set depending on each
384 convolution function. Default values for width are:
385
386 'box': 1 pixel
387 'sf': 3 pixels
388 'gauss': 1 pixel (width is used as HWHM)
389
390 func -- Function type ('box', 'sf', 'gauss').
391 width -- Width of convolution function. Default (-1) is to
392 choose pre-defined value for each convolution function.
393 """
[2671]394 fname = func.upper()
395 if fname == 'GAUSS' or fname == 'GJINC':
396 gw = str(gwidth)
397 jw = str(jwidth)
398 w = str(width)
399 if w[0] == '-': w = ''
400 #self.gridder._setfunc(fname, -1, w, gw, jw)
401 self.gridder._setfunc(fname,convsupport=-1,gwidth=gw,jwidth=jw,truncate=w)
402 else:
403 self.gridder._setfunc( func, convsupport=width )
[2593]404
405 def setWeight( self, weightType='uniform' ):
406 """
407 Set weight type. Possible options are 'uniform' (default),
408 'tint' (weight by integration time), 'tsys' (weight by
409 Tsys: 1/Tsys**2), and 'tintsys' (weight by integration time
410 as well as Tsys: tint/Tsys**2).
411
412 weightType -- weight type ('uniform', 'tint', 'tsys', 'tintsys')
413 """
414 self.gridder._setweight( weightType )
415
416 def enableClip( self ):
417 """
418 Enable min/max clipping.
419
420 By default, min/max clipping is disabled so that you should
421 call this method before actual gridding if you want to do
422 clipping.
423 """
424 self.gridder._enableclip()
425
426 def disableClip( self ):
427 """
428 Disable min/max clipping.
429 """
430 self.gridder._disableclip()
431
432 def grid( self ):
433 """
434 Actual gridding which will be done based on several user inputs.
435 """
436 self.gridder._grid()
437
438 def getResult( self ):
439 """
440 Return gridded data as a scantable.
441 """
[2594]442 tp = 0 if rcParams['scantable.storage']=='memory' else 1
443 return scantable( self.gridder._get( tp ), average=False )
[2671]444
445 def plotFunc(self, clear=True):
446 """
447 Support function to see the shape of current grid function.
[2593]448
[2671]449 clear -- clear panel if True. Default is True.
450 """
451 pl.figure(11)
452 if clear:
453 pl.clf()
454 f = self.gridder._getfunc()
455 convsampling = 100
456 a = numpy.arange(0,len(f)/convsampling,1./convsampling,dtype=float)
457 pl.plot(a,f,'.-')
458 pl.xlabel('pixel')
459 pl.ylabel('convFunc')
460
[2356]461class _SDGridPlotter:
[2373]462 def __init__( self, infile, outfile=None, ifno=-1 ):
[2390]463 if isinstance( infile, str ):
464 self.infile = [infile]
465 else:
466 self.infile = infile
[2356]467 self.outfile = outfile
468 if self.outfile is None:
[2390]469 self.outfile = self.infile[0].rstrip('/')+'.grid'
[2356]470 self.nx = -1
471 self.ny = -1
472 self.nchan = 0
[2360]473 self.npol = 0
474 self.pollist = []
[2356]475 self.cellx = 0.0
476 self.celly = 0.0
477 self.center = [0.0,0.0]
478 self.nonzero = [[0.0],[0.0]]
[2367]479 self.ifno = ifno
[2372]480 self.tablein = None
481 self.nrow = 0
482 self.blc = None
483 self.trc = None
[2356]484 self.get()
485
486 def get( self ):
487 s = scantable( self.outfile, average=False )
[2387]488 self.nchan = len(s._getspectrum(0))
[2356]489 nrow = s.nrow()
[2360]490 pols = numpy.ones( nrow, dtype=int )
[2356]491 for i in xrange(nrow):
[2360]492 pols[i] = s.getpol(i)
493 self.pollist, indices = numpy.unique( pols, return_inverse=True )
494 self.npol = len(self.pollist)
495 self.pollist = self.pollist[indices[:self.npol]]
496 #print 'pollist=',self.pollist
497 #print 'npol=',self.npol
498 #print 'nrow=',nrow
[2356]499
[2669]500 idx = 1
501 d0 = s.get_direction( 0 ).split()[-2]
502 d = s.get_direction(self.npol*idx)
503 while( d is not None \
504 and d.split()[-2] != d0):
[2356]505 idx += 1
[2669]506 d = s.get_direction(self.npol*idx)
[2367]507
[2372]508 self.nx = idx
509 self.ny = nrow / (self.npol * idx )
[2360]510 #print 'nx,ny=',self.nx,self.ny
[2372]511
512 self.blc = s.get_directionval( 0 )
513 self.trc = s.get_directionval( nrow-self.npol )
514 #print self.blc
515 #print self.trc
[2421]516 if nrow > 1:
517 incrx = s.get_directionval( self.npol )
518 incry = s.get_directionval( self.nx*self.npol )
519 else:
520 incrx = [0.0,0.0]
521 incry = [0.0,0.0]
[2372]522 self.cellx = abs( self.blc[0] - incrx[0] )
523 self.celly = abs( self.blc[1] - incry[1] )
[2360]524 #print 'cellx,celly=',self.cellx,self.celly
[2356]525
[2375]526 def plot( self, chan=-1, pol=-1, plotobs=False, plotgrid=False ):
[2360]527 if pol < 0:
528 opt = 'averaged over pol'
[2356]529 else:
[2360]530 opt = 'pol %s'%(pol)
[2419]531 if type(chan) is list:
532 opt += ', averaged over channel %s-%s'%(chan[0],chan[1])
533 elif chan < 0:
[2360]534 opt += ', averaged over channel'
535 else:
536 opt += ', channel %s'%(chan)
[2423]537 data = self.getData( chan, pol )
[2669]538 #data = numpy.fliplr( data )
[2360]539 title = 'Gridded Image (%s)'%(opt)
[2356]540 pl.figure(10)
541 pl.clf()
[2372]542 # plot grid position
[2375]543 if plotgrid:
544 x = numpy.arange(self.blc[0],self.trc[0]+0.5*self.cellx,self.cellx,dtype=float)
545 #print 'len(x)=',len(x)
546 #print 'x=',x
547 ybase = numpy.ones(self.nx,dtype=float)*self.blc[1]
548 #print 'len(ybase)=',len(ybase)
549 incr = self.celly
550 for iy in xrange(self.ny):
551 y = ybase + iy * incr
552 #print y
553 pl.plot(x,y,',',color='blue')
[2372]554 # plot observed position
[2375]555 if plotobs:
[2390]556 for i in xrange(len(self.infile)):
557 self.createTableIn( self.infile[i] )
558 irow = 0
559 while ( irow < self.nrow ):
560 chunk = self.getPointingChunk( irow )
561 #print chunk
562 pl.plot(chunk[0],chunk[1],',',color='green')
563 irow += chunk.shape[1]
564 #print irow
[2372]565 # show image
[2423]566 extent=[self.trc[0]+0.5*self.cellx,
567 self.blc[0]-0.5*self.cellx,
[2372]568 self.blc[1]-0.5*self.celly,
569 self.trc[1]+0.5*self.celly]
[2420]570 deccorr = 1.0/numpy.cos(0.5*(self.blc[1]+self.trc[1]))
[2356]571 pl.imshow(data,extent=extent,origin='lower',interpolation='nearest')
572 pl.colorbar()
573 pl.xlabel('R.A. [rad]')
574 pl.ylabel('Dec. [rad]')
[2420]575 ax = pl.axes()
576 ax.set_aspect(deccorr)
[2358]577 pl.title( title )
[2367]578
[2390]579 def createTableIn( self, tab ):
580 del self.tablein
581 self.tablein = scantable( tab, average=False )
[2387]582 if self.ifno < 0:
583 ifno = self.tablein.getif(0)
[2669]584 #print 'ifno=',ifno
[2387]585 else:
586 ifno = self.ifno
587 sel = selector()
588 sel.set_ifs( ifno )
[2390]589 self.tablein.set_selection( sel )
[2387]590 self.nchan = len(self.tablein._getspectrum(0))
[2390]591 self.nrow = self.tablein.nrow()
[2387]592 del sel
593
594
[2372]595 def getPointingChunk( self, irow ):
596 numchunk = 1000
597 nrow = min( self.nrow-irow, numchunk )
598 #print 'nrow=',nrow
599 v = numpy.zeros( (2,nrow), dtype=float )
600 idx = 0
601 for i in xrange(irow,irow+nrow):
602 d = self.tablein.get_directionval( i )
603 v[0,idx] = d[0]
604 v[1,idx] = d[1]
605 idx += 1
606 return v
607
[2367]608 def getData( self, chan=-1, pol=-1 ):
[2419]609 if type(chan) == list:
610 spectra = self.__chanAverage(start=chan[0],end=chan[1])
611 elif chan == -1:
[2367]612 spectra = self.__chanAverage()
613 else:
614 spectra = self.__chanIndex( chan )
[2372]615 data = spectra.reshape( (self.npol,self.ny,self.nx) )
[2367]616 if pol == -1:
617 retval = data.mean(axis=0)
618 else:
619 retval = data[pol]
620 return retval
621
[2419]622 def __chanAverage( self, start=-1, end=-1 ):
[2367]623 s = scantable( self.outfile, average=False )
[2372]624 nrow = s.nrow()
[2367]625 spectra = numpy.zeros( (self.npol,nrow/self.npol), dtype=float )
626 irow = 0
627 sp = [0 for i in xrange(self.nchan)]
[2419]628 if start < 0:
629 start = 0
630 if end < 0:
631 end = self.nchan
[2367]632 for i in xrange(nrow/self.npol):
633 for ip in xrange(self.npol):
[2419]634 sp = s._getspectrum( irow )[start:end]
[2367]635 spectra[ip,i] = numpy.mean( sp )
636 irow += 1
[2419]637
[2367]638 return spectra
639
640 def __chanIndex( self, idx ):
641 s = scantable( self.outfile, average=False )
[2372]642 nrow = s.nrow()
[2367]643 spectra = numpy.zeros( (self.npol,nrow/self.npol), dtype=float )
644 irow = 0
645 sp = [0 for i in xrange(self.nchan)]
646 for i in xrange(nrow/self.npol):
647 for ip in xrange(self.npol):
648 sp = s._getspectrum( irow )
649 spectra[ip,i] = sp[idx]
650 irow += 1
651 return spectra
652
653
Note: See TracBrowser for help on using the repository browser.