source: trunk/python/asapgrid.py

Last change on this file was 2894, checked in by Takeshi Nakazato, 10 years ago

New Development: No

JIRA Issue: Yes CAS-6121

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

Fix for latitude inversion.
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.blc[0]-0.5*self.cellx,
489                self.trc[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.