source: trunk/python/asapgrid.py@ 2595

Last change on this file since 2595 was 2594, 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: sdgrid unit test with sd.rcParamsscantable.storage='memory'

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Bug fix on STGrid2.


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:
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 'arcsec'. 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 = '%sarcsec'%(cellx)
133 if not isinstance( celly, str ):
134 celly = '%sarcsec'%(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': 1 pixel (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 asapgrid2:
227 """
228 The asapgrid class is defined to convolve data onto regular
229 spatial grid. Typical usage is as follows:
230
231 # create asapgrid instance with input scantable
232 s = scantable( 'testimage1.asap', average=False )
233 g = asapgrid( s )
234 # set IFNO if necessary
235 g.setIF( 0 )
236 # set POLNOs if necessary
237 g.setPolList( [0,1] )
238 # set SCANNOs if necessary
239 g.setScanList( [22,23,24] )
240 # define image with full specification
241 # you can skip some parameters (see help for defineImage)
242 g.defineImage( nx=12, ny=12, cellx='10arcsec', celly='10arcsec',
243 center='J2000 10h10m10s -5d05m05s' )
244 # set convolution function
245 g.setFunc( func='sf', width=3 )
246 # enable min/max clipping
247 g.enableClip()
248 # or, disable min/max clipping
249 #g.disableClip()
250 # actual gridding
251 g.grid()
252 # get result as scantable
253 sg = g.getResult()
254 """
255 def __init__( self, scantab ):
256 """
257 Create asapgrid instance.
258
259 scantab -- input data as a scantable or a list of scantables
260 to grid more than one data at once.
261 """
262 self.outfile = None
263 self.ifno = None
264 self.gridder = stgrid2()
265 self.setData( scantab )
266
267 def setData( self, scantab ):
268 """
269 Set data to be processed.
270
271 scantab -- input data as a scantable or a list of scantables
272 to grid more than one data at once.
273 """
274 if isinstance( scantab, scantable ):
275 self.gridder._setin( scantab )
276 else:
277 self.gridder._setfiles( scantab )
278 self.scantab = scantab
279
280 def setIF( self, ifno ):
281 """
282 Set IFNO to be processed. Currently, asapgrid allows to process
283 only one IFNO for one gridding run even if the data contains
284 multiple IFs. If you didn't specify IFNO, default value, which
285 is IFNO in the first spectrum, will be processed.
286
287 ifno -- IFNO to be processed.
288 """
289 self.ifno = ifno
290 self.gridder._setif( self.ifno )
291
292 def setPolList( self, pollist ):
293 """
294 Set list of polarization components you want to process.
295 If not specified, all POLNOs will be processed.
296
297 pollist -- list of POLNOs.
298 """
299 self.gridder._setpollist( pollist )
300
301 def setScanList( self, scanlist ):
302 """
303 Set list of scans you want to process. If not specified, all
304 scans will be processed.
305
306 scanlist -- list of SCANNOs.
307 """
308 self.gridder._setscanlist( scanlist )
309
310 def defineImage( self, nx=-1, ny=-1, cellx='', celly='', center='' ):
311 """
312 Define spatial grid.
313
314 First two parameters, nx and ny, define number of pixels of
315 the grid. If which of those is not specified, it will be set
316 to the same value as the other. If none of them are specified,
317 it will be determined from map extent and cell size.
318
319 Next two parameters, cellx and celly, define size of pixel.
320 You should set those parameters as string, which is constructed
321 numerical value and unit, e.g. '0.5arcmin', or numerical value.
322 If those values are specified as numerical value, their units
323 will be assumed to 'arcsec'. If which of those is not specified,
324 it will be set to the same value as the other. If none of them
325 are specified, it will be determined from map extent and number
326 of pixels, or set to '1arcmin' if neither nx nor ny is set.
327
328 The last parameter, center, define the central coordinate of
329 the grid. You should specify its value as a string, like,
330
331 'J2000 05h08m50s -16d23m30s'
332
333 or
334
335 'J2000 05:08:50 -16.23.30'
336
337 You can omit equinox when you specify center coordinate. In that
338 case, J2000 is assumed. If center is not specified, it will be
339 determined from the observed positions of input data.
340
341 nx -- number of pixels along x (R.A.) direction.
342 ny -- number of pixels along y (Dec.) direction.
343 cellx -- size of pixel in x (R.A.) direction.
344 celly -- size of pixel in y (Dec.) direction.
345 center -- central position of the grid.
346 """
347 if not isinstance( cellx, str ):
348 cellx = '%sarcsec'%(cellx)
349 if not isinstance( celly, str ):
350 celly = '%sarcsec'%(celly)
351 self.gridder._defineimage( nx, ny, cellx, celly, center )
352
353 def setFunc( self, func='box', width=-1 ):
354 """
355 Set convolution function. Possible options are 'box' (Box-car,
356 default), 'sf' (prolate spheroidal), and 'gauss' (Gaussian).
357 Width of convolution function can be set using width parameter.
358 By default (-1), width is automatically set depending on each
359 convolution function. Default values for width are:
360
361 'box': 1 pixel
362 'sf': 3 pixels
363 'gauss': 1 pixel (width is used as HWHM)
364
365 func -- Function type ('box', 'sf', 'gauss').
366 width -- Width of convolution function. Default (-1) is to
367 choose pre-defined value for each convolution function.
368 """
369 self.gridder._setfunc( func, width )
370
371 def setWeight( self, weightType='uniform' ):
372 """
373 Set weight type. Possible options are 'uniform' (default),
374 'tint' (weight by integration time), 'tsys' (weight by
375 Tsys: 1/Tsys**2), and 'tintsys' (weight by integration time
376 as well as Tsys: tint/Tsys**2).
377
378 weightType -- weight type ('uniform', 'tint', 'tsys', 'tintsys')
379 """
380 self.gridder._setweight( weightType )
381
382 def enableClip( self ):
383 """
384 Enable min/max clipping.
385
386 By default, min/max clipping is disabled so that you should
387 call this method before actual gridding if you want to do
388 clipping.
389 """
390 self.gridder._enableclip()
391
392 def disableClip( self ):
393 """
394 Disable min/max clipping.
395 """
396 self.gridder._disableclip()
397
398 def grid( self ):
399 """
400 Actual gridding which will be done based on several user inputs.
401 """
402 self.gridder._grid()
403
404 def getResult( self ):
405 """
406 Return gridded data as a scantable.
407 """
408 tp = 0 if rcParams['scantable.storage']=='memory' else 1
409 return scantable( self.gridder._get( tp ), average=False )
410
411class _SDGridPlotter:
412 def __init__( self, infile, outfile=None, ifno=-1 ):
413 if isinstance( infile, str ):
414 self.infile = [infile]
415 else:
416 self.infile = infile
417 self.outfile = outfile
418 if self.outfile is None:
419 self.outfile = self.infile[0].rstrip('/')+'.grid'
420 self.nx = -1
421 self.ny = -1
422 self.nchan = 0
423 self.npol = 0
424 self.pollist = []
425 self.cellx = 0.0
426 self.celly = 0.0
427 self.center = [0.0,0.0]
428 self.nonzero = [[0.0],[0.0]]
429 self.ifno = ifno
430 self.tablein = None
431 self.nrow = 0
432 self.blc = None
433 self.trc = None
434 self.get()
435
436 def get( self ):
437 s = scantable( self.outfile, average=False )
438 self.nchan = len(s._getspectrum(0))
439 nrow = s.nrow()
440 pols = numpy.ones( nrow, dtype=int )
441 for i in xrange(nrow):
442 pols[i] = s.getpol(i)
443 self.pollist, indices = numpy.unique( pols, return_inverse=True )
444 self.npol = len(self.pollist)
445 self.pollist = self.pollist[indices[:self.npol]]
446 #print 'pollist=',self.pollist
447 #print 'npol=',self.npol
448 #print 'nrow=',nrow
449
450 idx = 0
451 d0 = s.get_direction( 0 ).split()[-1]
452 while ( s.get_direction(self.npol*idx) is not None \
453 and s.get_direction(self.npol*idx).split()[-1] == d0 ):
454 idx += 1
455
456 self.nx = idx
457 self.ny = nrow / (self.npol * idx )
458 #print 'nx,ny=',self.nx,self.ny
459
460 self.blc = s.get_directionval( 0 )
461 self.trc = s.get_directionval( nrow-self.npol )
462 #print self.blc
463 #print self.trc
464 if nrow > 1:
465 incrx = s.get_directionval( self.npol )
466 incry = s.get_directionval( self.nx*self.npol )
467 else:
468 incrx = [0.0,0.0]
469 incry = [0.0,0.0]
470 self.cellx = abs( self.blc[0] - incrx[0] )
471 self.celly = abs( self.blc[1] - incry[1] )
472 #print 'cellx,celly=',self.cellx,self.celly
473
474 def plot( self, chan=-1, pol=-1, plotobs=False, plotgrid=False ):
475 if pol < 0:
476 opt = 'averaged over pol'
477 else:
478 opt = 'pol %s'%(pol)
479 if type(chan) is list:
480 opt += ', averaged over channel %s-%s'%(chan[0],chan[1])
481 elif chan < 0:
482 opt += ', averaged over channel'
483 else:
484 opt += ', channel %s'%(chan)
485 data = self.getData( chan, pol )
486 data = numpy.fliplr( data )
487 title = 'Gridded Image (%s)'%(opt)
488 pl.figure(10)
489 pl.clf()
490 # plot grid position
491 if plotgrid:
492 x = numpy.arange(self.blc[0],self.trc[0]+0.5*self.cellx,self.cellx,dtype=float)
493 #print 'len(x)=',len(x)
494 #print 'x=',x
495 ybase = numpy.ones(self.nx,dtype=float)*self.blc[1]
496 #print 'len(ybase)=',len(ybase)
497 incr = self.celly
498 for iy in xrange(self.ny):
499 y = ybase + iy * incr
500 #print y
501 pl.plot(x,y,',',color='blue')
502 # plot observed position
503 if plotobs:
504 for i in xrange(len(self.infile)):
505 self.createTableIn( self.infile[i] )
506 irow = 0
507 while ( irow < self.nrow ):
508 chunk = self.getPointingChunk( irow )
509 #print chunk
510 pl.plot(chunk[0],chunk[1],',',color='green')
511 irow += chunk.shape[1]
512 #print irow
513 # show image
514 extent=[self.trc[0]+0.5*self.cellx,
515 self.blc[0]-0.5*self.cellx,
516 self.blc[1]-0.5*self.celly,
517 self.trc[1]+0.5*self.celly]
518 deccorr = 1.0/numpy.cos(0.5*(self.blc[1]+self.trc[1]))
519 pl.imshow(data,extent=extent,origin='lower',interpolation='nearest')
520 pl.colorbar()
521 pl.xlabel('R.A. [rad]')
522 pl.ylabel('Dec. [rad]')
523 ax = pl.axes()
524 ax.set_aspect(deccorr)
525 pl.title( title )
526
527 def createTableIn( self, tab ):
528 del self.tablein
529 self.tablein = scantable( tab, average=False )
530 if self.ifno < 0:
531 ifno = self.tablein.getif(0)
532 print 'ifno=',ifno
533 else:
534 ifno = self.ifno
535 sel = selector()
536 sel.set_ifs( ifno )
537 self.tablein.set_selection( sel )
538 self.nchan = len(self.tablein._getspectrum(0))
539 self.nrow = self.tablein.nrow()
540 del sel
541
542
543 def getPointingChunk( self, irow ):
544 numchunk = 1000
545 nrow = min( self.nrow-irow, numchunk )
546 #print 'nrow=',nrow
547 v = numpy.zeros( (2,nrow), dtype=float )
548 idx = 0
549 for i in xrange(irow,irow+nrow):
550 d = self.tablein.get_directionval( i )
551 v[0,idx] = d[0]
552 v[1,idx] = d[1]
553 idx += 1
554 return v
555
556 def getData( self, chan=-1, pol=-1 ):
557 if type(chan) == list:
558 spectra = self.__chanAverage(start=chan[0],end=chan[1])
559 elif chan == -1:
560 spectra = self.__chanAverage()
561 else:
562 spectra = self.__chanIndex( chan )
563 data = spectra.reshape( (self.npol,self.ny,self.nx) )
564 if pol == -1:
565 retval = data.mean(axis=0)
566 else:
567 retval = data[pol]
568 return retval
569
570 def __chanAverage( self, start=-1, end=-1 ):
571 s = scantable( self.outfile, average=False )
572 nrow = s.nrow()
573 spectra = numpy.zeros( (self.npol,nrow/self.npol), dtype=float )
574 irow = 0
575 sp = [0 for i in xrange(self.nchan)]
576 if start < 0:
577 start = 0
578 if end < 0:
579 end = self.nchan
580 for i in xrange(nrow/self.npol):
581 for ip in xrange(self.npol):
582 sp = s._getspectrum( irow )[start:end]
583 spectra[ip,i] = numpy.mean( sp )
584 irow += 1
585
586 return spectra
587
588 def __chanIndex( self, idx ):
589 s = scantable( self.outfile, average=False )
590 nrow = s.nrow()
591 spectra = numpy.zeros( (self.npol,nrow/self.npol), dtype=float )
592 irow = 0
593 sp = [0 for i in xrange(self.nchan)]
594 for i in xrange(nrow/self.npol):
595 for ip in xrange(self.npol):
596 sp = s._getspectrum( irow )
597 spectra[ip,i] = sp[idx]
598 irow += 1
599 return spectra
600
601
Note: See TracBrowser for help on using the repository browser.