source: trunk/python/asapgrid.py@ 2421

Last change on this file since 2421 was 2421, checked in by Takeshi Nakazato, 13 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...

Support to plot only 1 pixel result.

Default unit for cell size is changed to 'arcsec' (was 'arcmin').

File size: 14.3 KB
RevLine 
[2356]1import numpy
[2367]2from asap import rcParams
[2356]3from asap.scantable import scantable
[2367]4from asap.selector import selector
[2356]5from asap._asap import stgrid
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
[2364]137 def setFunc( self, func='box', width=-1 ):
[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
147 'gauss': 3 pixels (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 """
[2364]153 self.gridder._setfunc( func, width )
[2356]154
[2361]155 def setWeight( self, weightType='uniform' ):
[2391]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 """
[2396]164 self.gridder._setweight( weightType )
[2361]165
[2396]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
[2356]182 def grid( self ):
[2391]183 """
184 Actual gridding which will be done based on several user inputs.
185 """
[2356]186 self.gridder._grid()
187
188 def save( self, outfile='' ):
[2391]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 """
[2356]195 self.outfile = self.gridder._save( outfile )
196
[2375]197 def plot( self, plotchan=-1, plotpol=-1, plotobs=False, plotgrid=False ):
[2391]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 """
[2367]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 )
[2375]219 plotter.plot( chan=plotchan, pol=plotpol, plotobs=plotobs, plotgrid=plotgrid )
[2367]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')
[2356]225
226class _SDGridPlotter:
[2373]227 def __init__( self, infile, outfile=None, ifno=-1 ):
[2390]228 if isinstance( infile, str ):
229 self.infile = [infile]
230 else:
231 self.infile = infile
[2356]232 self.outfile = outfile
233 if self.outfile is None:
[2390]234 self.outfile = self.infile[0].rstrip('/')+'.grid'
[2356]235 self.nx = -1
236 self.ny = -1
237 self.nchan = 0
[2360]238 self.npol = 0
239 self.pollist = []
[2356]240 self.cellx = 0.0
241 self.celly = 0.0
242 self.center = [0.0,0.0]
243 self.nonzero = [[0.0],[0.0]]
[2367]244 self.ifno = ifno
[2372]245 self.tablein = None
246 self.nrow = 0
247 self.blc = None
248 self.trc = None
[2356]249 self.get()
250
251 def get( self ):
252 s = scantable( self.outfile, average=False )
[2387]253 self.nchan = len(s._getspectrum(0))
[2356]254 nrow = s.nrow()
[2360]255 pols = numpy.ones( nrow, dtype=int )
[2356]256 for i in xrange(nrow):
[2360]257 pols[i] = s.getpol(i)
258 self.pollist, indices = numpy.unique( pols, return_inverse=True )
259 self.npol = len(self.pollist)
260 self.pollist = self.pollist[indices[:self.npol]]
261 #print 'pollist=',self.pollist
262 #print 'npol=',self.npol
263 #print 'nrow=',nrow
[2356]264
265 idx = 0
[2372]266 d0 = s.get_direction( 0 ).split()[-1]
[2421]267 while ( s.get_direction(self.npol*idx) is not None \
268 and s.get_direction(self.npol*idx).split()[-1] == d0 ):
[2356]269 idx += 1
[2367]270
[2372]271 self.nx = idx
272 self.ny = nrow / (self.npol * idx )
[2360]273 #print 'nx,ny=',self.nx,self.ny
[2372]274
275 self.blc = s.get_directionval( 0 )
276 self.trc = s.get_directionval( nrow-self.npol )
277 #print self.blc
278 #print self.trc
[2421]279 if nrow > 1:
280 incrx = s.get_directionval( self.npol )
281 incry = s.get_directionval( self.nx*self.npol )
282 else:
283 incrx = [0.0,0.0]
284 incry = [0.0,0.0]
[2372]285 self.cellx = abs( self.blc[0] - incrx[0] )
286 self.celly = abs( self.blc[1] - incry[1] )
[2360]287 #print 'cellx,celly=',self.cellx,self.celly
[2356]288
[2375]289 def plot( self, chan=-1, pol=-1, plotobs=False, plotgrid=False ):
[2360]290 if pol < 0:
291 opt = 'averaged over pol'
[2356]292 else:
[2360]293 opt = 'pol %s'%(pol)
[2419]294 if type(chan) is list:
295 opt += ', averaged over channel %s-%s'%(chan[0],chan[1])
296 elif chan < 0:
[2360]297 opt += ', averaged over channel'
298 else:
299 opt += ', channel %s'%(chan)
[2367]300 data = self.getData( chan, pol )
[2360]301 title = 'Gridded Image (%s)'%(opt)
[2356]302 pl.figure(10)
303 pl.clf()
[2372]304 # plot grid position
[2375]305 if plotgrid:
306 x = numpy.arange(self.blc[0],self.trc[0]+0.5*self.cellx,self.cellx,dtype=float)
307 #print 'len(x)=',len(x)
308 #print 'x=',x
309 ybase = numpy.ones(self.nx,dtype=float)*self.blc[1]
310 #print 'len(ybase)=',len(ybase)
311 incr = self.celly
312 for iy in xrange(self.ny):
313 y = ybase + iy * incr
314 #print y
315 pl.plot(x,y,',',color='blue')
[2372]316 # plot observed position
[2375]317 if plotobs:
[2390]318 for i in xrange(len(self.infile)):
319 self.createTableIn( self.infile[i] )
320 irow = 0
321 while ( irow < self.nrow ):
322 chunk = self.getPointingChunk( irow )
323 #print chunk
324 pl.plot(chunk[0],chunk[1],',',color='green')
325 irow += chunk.shape[1]
326 #print irow
[2372]327 # show image
328 extent=[self.blc[0]-0.5*self.cellx,
329 self.trc[0]+0.5*self.cellx,
330 self.blc[1]-0.5*self.celly,
331 self.trc[1]+0.5*self.celly]
[2420]332 deccorr = 1.0/numpy.cos(0.5*(self.blc[1]+self.trc[1]))
[2356]333 pl.imshow(data,extent=extent,origin='lower',interpolation='nearest')
334 pl.colorbar()
335 pl.xlabel('R.A. [rad]')
336 pl.ylabel('Dec. [rad]')
[2420]337 ax = pl.axes()
338 ax.set_aspect(deccorr)
[2358]339 pl.title( title )
[2367]340
[2390]341 def createTableIn( self, tab ):
342 del self.tablein
343 self.tablein = scantable( tab, average=False )
[2387]344 if self.ifno < 0:
345 ifno = self.tablein.getif(0)
346 print 'ifno=',ifno
347 else:
348 ifno = self.ifno
349 sel = selector()
350 sel.set_ifs( ifno )
[2390]351 self.tablein.set_selection( sel )
[2387]352 self.nchan = len(self.tablein._getspectrum(0))
[2390]353 self.nrow = self.tablein.nrow()
[2387]354 del sel
355
356
[2372]357 def getPointingChunk( self, irow ):
358 numchunk = 1000
359 nrow = min( self.nrow-irow, numchunk )
360 #print 'nrow=',nrow
361 v = numpy.zeros( (2,nrow), dtype=float )
362 idx = 0
363 for i in xrange(irow,irow+nrow):
364 d = self.tablein.get_directionval( i )
365 v[0,idx] = d[0]
366 v[1,idx] = d[1]
367 idx += 1
368 return v
369
[2367]370 def getData( self, chan=-1, pol=-1 ):
[2419]371 if type(chan) == list:
372 spectra = self.__chanAverage(start=chan[0],end=chan[1])
373 elif chan == -1:
[2367]374 spectra = self.__chanAverage()
375 else:
376 spectra = self.__chanIndex( chan )
[2372]377 data = spectra.reshape( (self.npol,self.ny,self.nx) )
[2367]378 if pol == -1:
379 retval = data.mean(axis=0)
380 else:
381 retval = data[pol]
382 return retval
383
[2419]384 def __chanAverage( self, start=-1, end=-1 ):
[2367]385 s = scantable( self.outfile, average=False )
[2372]386 nrow = s.nrow()
[2367]387 spectra = numpy.zeros( (self.npol,nrow/self.npol), dtype=float )
388 irow = 0
389 sp = [0 for i in xrange(self.nchan)]
[2419]390 if start < 0:
391 start = 0
392 if end < 0:
393 end = self.nchan
[2367]394 for i in xrange(nrow/self.npol):
395 for ip in xrange(self.npol):
[2419]396 sp = s._getspectrum( irow )[start:end]
[2367]397 spectra[ip,i] = numpy.mean( sp )
398 irow += 1
[2419]399
[2367]400 return spectra
401
402 def __chanIndex( self, idx ):
403 s = scantable( self.outfile, average=False )
[2372]404 nrow = s.nrow()
[2367]405 spectra = numpy.zeros( (self.npol,nrow/self.npol), dtype=float )
406 irow = 0
407 sp = [0 for i in xrange(self.nchan)]
408 for i in xrange(nrow/self.npol):
409 for ip in xrange(self.npol):
410 sp = s._getspectrum( irow )
411 spectra[ip,i] = sp[idx]
412 irow += 1
413 return spectra
414
415
Note: See TracBrowser for help on using the repository browser.