source: trunk/python/asapgrid.py @ 2593

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

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: Yes/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...

New class asapgrid2 is defined. Difference from asapgrid class is
that input and output is a scantable instance, not a table name.


File size: 20.6 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        return scantable( self.gridder._get(), average=False )
409
410class _SDGridPlotter:
411    def __init__( self, infile, outfile=None, ifno=-1 ):
412        if isinstance( infile, str ):
413            self.infile = [infile]
414        else:
415            self.infile = infile
416        self.outfile = outfile
417        if self.outfile is None:
418            self.outfile = self.infile[0].rstrip('/')+'.grid'
419        self.nx = -1
420        self.ny = -1
421        self.nchan = 0
422        self.npol = 0
423        self.pollist = []
424        self.cellx = 0.0
425        self.celly = 0.0
426        self.center = [0.0,0.0]
427        self.nonzero = [[0.0],[0.0]]
428        self.ifno = ifno
429        self.tablein = None
430        self.nrow = 0
431        self.blc = None
432        self.trc = None
433        self.get()
434
435    def get( self ):
436        s = scantable( self.outfile, average=False )
437        self.nchan = len(s._getspectrum(0))
438        nrow = s.nrow()
439        pols = numpy.ones( nrow, dtype=int )
440        for i in xrange(nrow):
441            pols[i] = s.getpol(i)
442        self.pollist, indices = numpy.unique( pols, return_inverse=True )
443        self.npol = len(self.pollist)
444        self.pollist = self.pollist[indices[:self.npol]]
445        #print 'pollist=',self.pollist
446        #print 'npol=',self.npol
447        #print 'nrow=',nrow
448
449        idx = 0
450        d0 = s.get_direction( 0 ).split()[-1]
451        while ( s.get_direction(self.npol*idx) is not None \
452                and s.get_direction(self.npol*idx).split()[-1] == d0 ):
453            idx += 1
454       
455        self.nx = idx
456        self.ny = nrow / (self.npol * idx )
457        #print 'nx,ny=',self.nx,self.ny
458
459        self.blc = s.get_directionval( 0 )
460        self.trc = s.get_directionval( nrow-self.npol )
461        #print self.blc
462        #print self.trc
463        if nrow > 1:
464            incrx = s.get_directionval( self.npol )
465            incry = s.get_directionval( self.nx*self.npol )
466        else:
467            incrx = [0.0,0.0]
468            incry = [0.0,0.0]
469        self.cellx = abs( self.blc[0] - incrx[0] )
470        self.celly = abs( self.blc[1] - incry[1] )
471        #print 'cellx,celly=',self.cellx,self.celly
472
473    def plot( self, chan=-1, pol=-1, plotobs=False, plotgrid=False ):
474        if pol < 0:
475            opt = 'averaged over pol'
476        else:
477            opt = 'pol %s'%(pol)
478        if type(chan) is list:
479            opt += ', averaged over channel %s-%s'%(chan[0],chan[1])
480        elif chan < 0:
481            opt += ', averaged over channel'
482        else:
483            opt += ', channel %s'%(chan)
484        data = self.getData( chan, pol )
485        data = numpy.fliplr( data )
486        title = 'Gridded Image (%s)'%(opt)
487        pl.figure(10)
488        pl.clf()
489        # plot grid position
490        if plotgrid:
491            x = numpy.arange(self.blc[0],self.trc[0]+0.5*self.cellx,self.cellx,dtype=float)
492            #print 'len(x)=',len(x)
493            #print 'x=',x
494            ybase = numpy.ones(self.nx,dtype=float)*self.blc[1]
495            #print 'len(ybase)=',len(ybase)
496            incr = self.celly
497            for iy in xrange(self.ny):
498                y = ybase + iy * incr
499                #print y
500                pl.plot(x,y,',',color='blue')
501        # plot observed position
502        if plotobs:
503            for i in xrange(len(self.infile)):
504                self.createTableIn( self.infile[i] )
505                irow = 0
506                while ( irow < self.nrow ):
507                    chunk = self.getPointingChunk( irow )
508                    #print chunk
509                    pl.plot(chunk[0],chunk[1],',',color='green')
510                    irow += chunk.shape[1]
511                    #print irow
512        # show image
513        extent=[self.trc[0]+0.5*self.cellx,
514                self.blc[0]-0.5*self.cellx,
515                self.blc[1]-0.5*self.celly,
516                self.trc[1]+0.5*self.celly]
517        deccorr = 1.0/numpy.cos(0.5*(self.blc[1]+self.trc[1]))
518        pl.imshow(data,extent=extent,origin='lower',interpolation='nearest')
519        pl.colorbar()
520        pl.xlabel('R.A. [rad]')
521        pl.ylabel('Dec. [rad]')
522        ax = pl.axes()
523        ax.set_aspect(deccorr)
524        pl.title( title )
525
526    def createTableIn( self, tab ):
527        del self.tablein
528        self.tablein = scantable( tab, average=False )
529        if self.ifno < 0:
530            ifno = self.tablein.getif(0)
531            print 'ifno=',ifno
532        else:
533            ifno = self.ifno
534        sel = selector()
535        sel.set_ifs( ifno )
536        self.tablein.set_selection( sel )
537        self.nchan = len(self.tablein._getspectrum(0))
538        self.nrow = self.tablein.nrow()
539        del sel
540       
541
542    def getPointingChunk( self, irow ):
543        numchunk = 1000
544        nrow = min( self.nrow-irow, numchunk )
545        #print 'nrow=',nrow
546        v = numpy.zeros( (2,nrow), dtype=float )
547        idx = 0
548        for i in xrange(irow,irow+nrow):
549            d = self.tablein.get_directionval( i )
550            v[0,idx] = d[0]
551            v[1,idx] = d[1]
552            idx += 1
553        return v
554
555    def getData( self, chan=-1, pol=-1 ):
556        if type(chan) == list:
557            spectra = self.__chanAverage(start=chan[0],end=chan[1])
558        elif chan == -1:
559            spectra = self.__chanAverage()
560        else:
561            spectra = self.__chanIndex( chan )
562        data = spectra.reshape( (self.npol,self.ny,self.nx) )
563        if pol == -1:
564            retval = data.mean(axis=0)
565        else:
566            retval = data[pol]
567        return retval
568
569    def __chanAverage( self, start=-1, end=-1 ):
570        s = scantable( self.outfile, average=False )
571        nrow = s.nrow()
572        spectra = numpy.zeros( (self.npol,nrow/self.npol), dtype=float )
573        irow = 0
574        sp = [0 for i in xrange(self.nchan)]
575        if start < 0:
576            start = 0
577        if end < 0:
578            end = self.nchan
579        for i in xrange(nrow/self.npol):
580            for ip in xrange(self.npol):
581                sp = s._getspectrum( irow )[start:end]
582                spectra[ip,i] = numpy.mean( sp )
583                irow += 1
584           
585        return spectra
586
587    def __chanIndex( self, idx ):
588        s = scantable( self.outfile, average=False )
589        nrow = s.nrow()
590        spectra = numpy.zeros( (self.npol,nrow/self.npol), dtype=float )
591        irow = 0
592        sp = [0 for i in xrange(self.nchan)]
593        for i in xrange(nrow/self.npol):
594            for ip in xrange(self.npol):
595                sp = s._getspectrum( irow )
596                spectra[ip,i] = sp[idx]
597                irow += 1
598        return spectra
599       
600           
Note: See TracBrowser for help on using the repository browser.