source: trunk/python/asapgrid.py @ 2391

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

New Development: No

JIRA Issue: Yes CAS-2816

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

Added documentation of toolkit level.


File size: 13.2 KB
Line 
1import numpy
2from asap import rcParams
3from asap.scantable import scantable
4from asap.selector import selector
5from asap._asap import stgrid
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       # actual gridding
29       g.grid()
30       # save result
31       g.save( outfile='grid.asap' )
32       # plot result
33       g.plot( plotchan=1246, plotpol=-1, plotgrid=True, plotobs=True )
34    """
35    def __init__( self, infile ):
36        """
37        Create asapgrid instance.
38
39        infile -- input data as a string or string list if you want
40                  to grid more than one data at once. 
41        """
42        self.outfile = None
43        self.ifno = None
44        self.gridder = stgrid()
45        self.setData( infile )
46
47    def setData( self, infile ):
48        """
49        Set data to be processed.
50
51        infile -- input data as a string or string list if you want
52                  to grid more than one data at once. 
53        """
54        if isinstance( infile, str ):
55            self.gridder._setin( infile )
56        else:
57            self.gridder._setfiles( infile )
58        self.infile = infile
59
60    def setIF( self, ifno ):
61        """
62        Set IFNO to be processed. Currently, asapgrid allows to process
63        only one IFNO for one gridding run even if the data contains
64        multiple IFs. If you didn't specify IFNO, default value, which
65        is IFNO in the first spectrum, will be processed.
66
67        ifno -- IFNO to be processed.
68        """
69        self.ifno = ifno
70        self.gridder._setif( self.ifno )
71
72    def setPolList( self, pollist ):
73        """
74        Set list of polarization components you want to process.
75        If not specified, all POLNOs will be processed.
76
77        pollist -- list of POLNOs.
78        """
79        self.gridder._setpollist( pollist )
80
81    def setScanList( self, scanlist ):
82        """
83        Set list of scans you want to process. If not specified, all
84        scans will be processed.
85
86        scanlist -- list of SCANNOs.
87        """
88        self.gridder._setscanlist( scanlist )
89
90    def defineImage( self, nx=-1, ny=-1, cellx='', celly='', center='' ):
91        """
92        Define spatial grid.
93
94        First two parameters, nx and ny, define number of pixels of
95        the grid. If which of those is not specified, it will be set
96        to the same value as the other. If none of them are specified,
97        it will be determined from map extent and cell size.
98
99        Next two parameters, cellx and celly, define size of pixel.
100        You should set those parameters as string, which is constructed
101        numerical value and unit, e.g. '0.5arcmin', or numerical value.
102        If those values are specified as numerical value, their units
103        will be assumed to 'arcmin'. If which of those is not specified,
104        it will be set to the same value as the other. If none of them
105        are specified, it will be determined from map extent and number
106        of pixels, or set to '1arcmin' if neither nx nor ny is set.
107
108        The last parameter, center, define the central coordinate of
109        the grid. You should specify its value as a string, like,
110
111           'J2000 05h08m50s -16d23m30s'
112
113        or
114
115           'J2000 05:08:50 -16.23.30'
116
117        You can omit equinox when you specify center coordinate. In that
118        case, J2000 is assumed. If center is not specified, it will be
119        determined from the observed positions of input data.
120
121        nx -- number of pixels along x (R.A.) direction.
122        ny -- number of pixels along y (Dec.) direction.
123        cellx -- size of pixel in x (R.A.) direction.
124        celly -- size of pixel in y (Dec.) direction.
125        center -- central position of the grid.
126        """
127        if not isinstance( cellx, str ):
128            cellx = '%sarcmin'%(cellx)
129        if not isinstance( celly, str ):
130            celly = '%sarcmin'%(celly)
131        self.gridder._defineimage( nx, ny, cellx, celly, center )
132
133    def setFunc( self, func='box', width=-1 ):
134        """
135        Set convolution function. Possible options are 'box' (Box-car,
136        default), 'sf' (prolate spheroidal), and 'gauss' (Gaussian).
137        Width of convolution function can be set using width parameter.
138        By default (-1), width is automatically set depending on each
139        convolution function. Default values for width are:
140
141           'box': 1 pixel
142           'sf': 3 pixels
143           'gauss': 3 pixels (width is used as HWHM)
144
145        func -- Function type ('box', 'sf', 'gauss').
146        width -- Width of convolution function. Default (-1) is to
147                 choose pre-defined value for each convolution function.
148        """
149        self.gridder._setfunc( func, width )
150
151    def setWeight( self, weightType='uniform' ):
152        """
153        Set weight type. Possible options are 'uniform' (default),
154        'tint' (weight by integration time), 'tsys' (weight by
155        Tsys: 1/Tsys**2), and 'tintsys' (weight by integration time
156        as well as Tsys: tint/Tsys**2).
157
158        weightType -- weight type ('uniform', 'tint', 'tsys', 'tintsys')
159        """
160        self.gridder._setweight( weightType )
161
162    def grid( self ):
163        """
164        Actual gridding which will be done based on several user inputs.
165        """
166        self.gridder._grid()
167
168    def save( self, outfile='' ):
169        """
170        Save result. By default, output data name will be constructed
171        from first element of input data name list (e.g. 'input.asap.grid').
172
173        outfile -- output data name.
174        """
175        self.outfile = self.gridder._save( outfile )
176
177    def plot( self, plotchan=-1, plotpol=-1, plotobs=False, plotgrid=False ):
178        """
179        Plot gridded data.
180
181        plotchan -- Which channel you want to plot. Default (-1) is
182                    to average all the channels.
183        plotpol -- Which polarization component you want to plot.
184                   Default (-1) is to average all the polarization
185                   components.
186        plotobs -- Also plot observed position if True. Default
187                   is False. Setting True for large amount of spectra
188                   might be time consuming.
189        plotgrid -- Also plot grid center if True. Default is False.
190                    Setting True for large number of grids might be
191                    time consuming.
192        """
193        import time
194        t0=time.time()
195        # to load scantable on disk
196        storg = rcParams['scantable.storage']
197        rcParams['scantable.storage'] = 'disk'
198        plotter = _SDGridPlotter( self.infile, self.outfile, self.ifno )
199        plotter.plot( chan=plotchan, pol=plotpol, plotobs=plotobs, plotgrid=plotgrid )
200        # back to original setup
201        rcParams['scantable.storage'] = storg
202        t1=time.time()
203        asaplog.push('plot: elapsed time %s sec'%(t1-t0))
204        asaplog.post('DEBUG','asapgrid.plot')
205       
206class _SDGridPlotter:
207    def __init__( self, infile, outfile=None, ifno=-1 ):
208        if isinstance( infile, str ):
209            self.infile = [infile]
210        else:
211            self.infile = infile
212        self.outfile = outfile
213        if self.outfile is None:
214            self.outfile = self.infile[0].rstrip('/')+'.grid'
215        self.nx = -1
216        self.ny = -1
217        self.nchan = 0
218        self.npol = 0
219        self.pollist = []
220        self.cellx = 0.0
221        self.celly = 0.0
222        self.center = [0.0,0.0]
223        self.nonzero = [[0.0],[0.0]]
224        self.ifno = ifno
225        self.tablein = None
226        self.nrow = 0
227        self.blc = None
228        self.trc = None
229        self.get()
230
231    def get( self ):
232        s = scantable( self.outfile, average=False )
233        self.nchan = len(s._getspectrum(0))
234        nrow = s.nrow()
235        pols = numpy.ones( nrow, dtype=int )
236        for i in xrange(nrow):
237            pols[i] = s.getpol(i)
238        self.pollist, indices = numpy.unique( pols, return_inverse=True )
239        self.npol = len(self.pollist)
240        self.pollist = self.pollist[indices[:self.npol]]
241        #print 'pollist=',self.pollist
242        #print 'npol=',self.npol
243        #print 'nrow=',nrow
244
245        idx = 0
246        d0 = s.get_direction( 0 ).split()[-1]
247        while ( s.get_direction(self.npol*idx).split()[-1] == d0 ): 
248            idx += 1
249       
250        self.nx = idx
251        self.ny = nrow / (self.npol * idx )
252        #print 'nx,ny=',self.nx,self.ny
253
254        self.blc = s.get_directionval( 0 )
255        self.trc = s.get_directionval( nrow-self.npol )
256        #print self.blc
257        #print self.trc
258        incrx = s.get_directionval( self.npol )
259        incry = s.get_directionval( self.nx*self.npol )
260        self.cellx = abs( self.blc[0] - incrx[0] )
261        self.celly = abs( self.blc[1] - incry[1] )
262        #print 'cellx,celly=',self.cellx,self.celly
263
264    def plot( self, chan=-1, pol=-1, plotobs=False, plotgrid=False ):
265        if pol < 0:
266            opt = 'averaged over pol'
267        else:
268            opt = 'pol %s'%(pol)
269        if chan < 0:
270            opt += ', averaged over channel'
271        else:
272            opt += ', channel %s'%(chan)
273        data = self.getData( chan, pol )
274        title = 'Gridded Image (%s)'%(opt)
275        pl.figure(10)
276        pl.clf()
277        # plot grid position
278        if plotgrid:
279            x = numpy.arange(self.blc[0],self.trc[0]+0.5*self.cellx,self.cellx,dtype=float)
280            #print 'len(x)=',len(x)
281            #print 'x=',x
282            ybase = numpy.ones(self.nx,dtype=float)*self.blc[1]
283            #print 'len(ybase)=',len(ybase)
284            incr = self.celly
285            for iy in xrange(self.ny):
286                y = ybase + iy * incr
287                #print y
288                pl.plot(x,y,',',color='blue')
289        # plot observed position
290        if plotobs:
291            for i in xrange(len(self.infile)):
292                self.createTableIn( self.infile[i] )
293                irow = 0
294                while ( irow < self.nrow ):
295                    chunk = self.getPointingChunk( irow )
296                    #print chunk
297                    pl.plot(chunk[0],chunk[1],',',color='green')
298                    irow += chunk.shape[1]
299                    #print irow
300        # show image
301        extent=[self.blc[0]-0.5*self.cellx,
302                self.trc[0]+0.5*self.cellx,
303                self.blc[1]-0.5*self.celly,
304                self.trc[1]+0.5*self.celly]
305        pl.imshow(data,extent=extent,origin='lower',interpolation='nearest')
306        pl.colorbar()
307        pl.xlabel('R.A. [rad]')
308        pl.ylabel('Dec. [rad]')
309        pl.title( title )
310
311    def createTableIn( self, tab ):
312        del self.tablein
313        self.tablein = scantable( tab, average=False )
314        if self.ifno < 0:
315            ifno = self.tablein.getif(0)
316            print 'ifno=',ifno
317        else:
318            ifno = self.ifno
319        sel = selector()
320        sel.set_ifs( ifno )
321        self.tablein.set_selection( sel )
322        self.nchan = len(self.tablein._getspectrum(0))
323        self.nrow = self.tablein.nrow()
324        del sel
325       
326
327    def getPointingChunk( self, irow ):
328        numchunk = 1000
329        nrow = min( self.nrow-irow, numchunk )
330        #print 'nrow=',nrow
331        v = numpy.zeros( (2,nrow), dtype=float )
332        idx = 0
333        for i in xrange(irow,irow+nrow):
334            d = self.tablein.get_directionval( i )
335            v[0,idx] = d[0]
336            v[1,idx] = d[1]
337            idx += 1
338        return v
339
340    def getData( self, chan=-1, pol=-1 ):
341        if chan == -1:
342            spectra = self.__chanAverage()
343        else:
344            spectra = self.__chanIndex( chan )
345        data = spectra.reshape( (self.npol,self.ny,self.nx) )
346        if pol == -1:
347            retval = data.mean(axis=0)
348        else:
349            retval = data[pol]
350        return retval
351
352    def __chanAverage( self ):
353        s = scantable( self.outfile, average=False )
354        nrow = s.nrow()
355        spectra = numpy.zeros( (self.npol,nrow/self.npol), dtype=float )
356        irow = 0
357        sp = [0 for i in xrange(self.nchan)]
358        for i in xrange(nrow/self.npol):
359            for ip in xrange(self.npol):
360                sp = s._getspectrum( irow )
361                spectra[ip,i] = numpy.mean( sp )
362                irow += 1
363        return spectra
364
365    def __chanIndex( self, idx ):
366        s = scantable( self.outfile, average=False )
367        nrow = s.nrow()
368        spectra = numpy.zeros( (self.npol,nrow/self.npol), dtype=float )
369        irow = 0
370        sp = [0 for i in xrange(self.nchan)]
371        for i in xrange(nrow/self.npol):
372            for ip in xrange(self.npol):
373                sp = s._getspectrum( irow )
374                spectra[ip,i] = sp[idx]
375                irow += 1
376        return spectra
377       
378           
Note: See TracBrowser for help on using the repository browser.