Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/simutil.py: 4%
2174 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-10-31 19:53 +0000
1# geodesy and pointing and other helper functions that are useful
2# to be available outside of the simdata task
3# geodesy from NGS: http://www.ngs.noaa.gov/TOOLS/program_descriptions.html
4import os
5import shutil
6import pylab as pl
7import glob
8import re
9from scipy.stats import scoreatpercentile
10import scipy.special as spspec
11import scipy.signal as spsig
12import scipy.interpolate as spintrp
13from collections import OrderedDict
15from casatools import table, image, imagepol, regionmanager, calibrater, measures, quanta, coordsys, componentlist, simulator, synthesisutils, ctsys
16from casatasks import casalog, tclean
17from casatasks.private.cleanhelper import cleanhelper
18tb = table( )
19ia = image( )
20po = imagepol( )
21rg = regionmanager( )
22cb = calibrater( )
23me = measures( )
24qa = quanta( )
25cs = coordsys( )
26cl = componentlist( )
27sm = simulator( )
28_su = synthesisutils( )
30# functions defined outside of the simutil class
31def is_array_type(value):
32 array_type = [list, tuple, pl.ndarray]
33 if type(value) in array_type:
34 return True
35 else:
36 return False
38# the method which returns a string of task function call with parameter
39def get_taskstr(taskname, params):
40 """
41 Returns a task call string.
42 Arguments
43 taskname : task name string
44 params : a dictionary of parameter (key) and value pairs.
45 Example
46 get_taskstr('mytask', {'vis': 'foo.ms', 'factor': 2.0})
47 returns a string, 'mytask(vis='foo.ms', factor=2.0)'
48 """
49 out = ("%s(" % taskname)
50 sep = ", "
51 for key, val in params.items():
52 out += (key + "=" + __get_str(val) + sep)
54 return ( out.rstrip(sep) + ")" )
56def __get_str(paramval):
57 if type(paramval) == str:
58 return ("'%s'" % paramval)
59 # else
60 return str(paramval)
63class compositenumber:
64 def __init__(self, maxval=100):
65 self.generate(maxval)
66 def generate(self,maxval):
67 n2 = int(log(float(maxval))/log(2.0) + 1) +1
68 n3 = int(log(float(maxval))/log(3.0) + 1) +1
69 n5 = int(log(float(maxval))/log(5.0) + 1) +1
70 itsnumbers=pl.zeros(n2*n3*n5)
71 n = 0
72 for i2 in range(n2):
73 for i3 in range(n3):
74 for i5 in range(n5):
75 composite=( 2.**i2 * 3.**i3 * 5.**i5 )
76 itsnumbers[n] = composite
77 #casalog.post... i2,i3,i5,composite
78 n=n+1
79 itsnumbers.sort()
80 maxi=0
81 while maxi<(n2*n3*n5) and itsnumbers[maxi]<=maxval: maxi=maxi+1
82 self.itsnumbers=pl.int64(itsnumbers[0:maxi])
83 def list(self):
84 casalog.post(self.itsnumbers)
85 def nextlarger(self,x):
86 if x>max(self.itsnumbers): self.generate(2*x)
87 xi=0
88 n=self.itsnumbers.__len__()
89 while xi<n and self.itsnumbers[xi]<x: xi=xi+1
90 return self.itsnumbers[xi]
96class simutil:
97 """
98 simutil contains methods to facilitate simulation.
99 To use these, create a simutil instance e.g., in CASA 6:
100 CASA> from casatasks.private import simutil
101 CASA> u=simutil.simutil()
102 CASA> x,y,z,d,padnames,antnames,telescope,posobs = u.readantenna("myconfig.cfg")
103 (in CASA5: )
104 CASA> from simutil import simutil
105 CASA> u=simutil.simutil()
106 """
107 def __init__(self, direction="",
108 centerfreq=qa.quantity("245GHz"),
109 bandwidth=qa.quantity("1GHz"),
110 totaltime=qa.quantity("0h"),
111 verbose=False):
112 self.direction=direction
113 self.verbose=verbose
114 self.centerfreq=centerfreq
115 self.bandwidth=bandwidth
116 self.totaltime=totaltime
117 self.pmulti=0 # rows, cols, currsubplot
120 def newfig(self,multi=0,show=True): # new graphics window/file
121 if show:
122 pl.ion() # creates a fig if ness
123 else:
124 pl.ioff()
125 pl.clf()
127 if multi!=0:
128 if type(multi)!=type([]):
129 self.msg("internal error setting multi-panel figure with multi="+str(multi),priority="warn")
130 if len(multi)!=3:
131 self.msg("internal error setting multi-panel figure with multi="+str(multi),priority="warn")
132 self.pmulti=multi
133 pl.subplot(multi[0],multi[1],multi[2])
134 pl.subplots_adjust(left=0.05,right=0.98,bottom=0.09,top=0.95,hspace=0.2,wspace=0.2)
137 ###########################################################
139 def nextfig(self): # advance subwindow
140 ax=pl.gca()
141 l=ax.get_xticklabels()
142 pl.setp(l,fontsize="x-small")
143 l=ax.get_yticklabels()
144 pl.setp(l,fontsize="x-small")
145 if self.pmulti!=0:
146 self.pmulti[2] += 1
147 multi=self.pmulti
148 if multi[2] <= multi[0]*multi[1]:
149 pl.subplot(multi[0],multi[1],multi[2])
150 # consider pl.draw() here - may be slow
152 ###########################################################
154 def endfig(self,show=True,filename=""): # set margins to smaller, save to file if required
155 ax=pl.gca()
156 l=ax.get_xticklabels()
157 pl.setp(l,fontsize="x-small")
158 l=ax.get_yticklabels()
159 pl.setp(l,fontsize="x-small")
160 if show:
161 pl.draw()
162 if len(filename)>0:
163 pl.savefig(filename)
164 pl.ioff()
168 ###########################################################
170 def msg(self, s, origin=None, priority=None):
171 # everything goes to logger with priority=priority
172 # priority error: raise an exception,
173 # casa 5.0.0: default is now no term unless toterm=True
174 # priority warn: change color to magenta, send to terminal
175 # priority info: change color to green, send to terminal
176 # priority none: not normally to terminal unless verbose=True
178 # ansi color codes:
179 # Foreground colors
180 # 30 Black
181 # 31 Red
182 # 32 Green
183 # 33 Yellow
184 # 34 Blue
185 # 35 Magenta
186 # 36 Cyan
187 # 37 White
189 clr=""
190 if self.verbose:
191 toterm=True
192 else:
193 toterm=False
195 if priority==None:
196 priority="INFO"
197 else:
198 priority=priority.upper()
199 #toterm=True
200 if priority=="INFO":
201 clr="\x1b[32m"
202 elif priority.count("WARN")>0:
203 clr="\x1b[35m"
204 #toterm=True
205 if self.verbose: priority="INFO" # otherwise casalog will spew to term also.
206 elif priority=="ERROR":
207 clr="\x1b[31m"
208 toterm=False # casalog spews severe to term already
209 else:
210 if not (priority=="DEBUG" or priority[:-1]=="DEBUG"):
211 priority="INFO"
212 bw="\x1b[0m"
214 if toterm:
215 if self.isreport():
216 if origin:
217 self.report.write("["+origin+"] "+s+"\n")
218 else:
219 self.report.write(s+"\n")
221 if s.count("ERROR"):
222 foo=s.split("ERROR")
223 s=foo[0]+"\x1b[31mERROR\x1b[0m"+foo[1]
224 if s.count("WARNING"):
225 foo=s.split("WARNING")
226 s=foo[0]+"\x1b[35mWARNING\x1b[0m"+foo[1]
228 if origin:
229 casalog.post(clr+"["+origin+"] "+bw+s)
230 else:
231 casalog.post(s)
234 if priority=="ERROR":
235 raise Exception(s)
236 else:
237 if origin==None:
238 origin="simutil"
239 casalog.post(s,priority=priority,origin=origin)
242 ###########################################################
244 def isreport(self):
245 # is there an open report file?
246 try:
247 if self.report.name == self.reportfile:
248 return True
249 else:
250 return False
251 except:
252 return False
255 def openreport(self):
256 try:
257 if os.path.exists(self.reportfile):
258 self.report=open(self.reportfile,"a")
259 #self.msg("Report file "+self.reportfile+"already exists - delete or change reportfile",priority="ERROR",origin="simutil")
260 else:
261 self.report=open(self.reportfile,"w")
262 except:
263 self.msg("Can't open reportfile because it's not defined",priority="ERROR",origin="simutil")
266 def closereport(self):
267 self.report.close()
273 ###########################################################
275 def isquantity(self,s,halt=True):
276 if type(s)!=type([]):
277 t=[s]
278 else:
279 t=s
280 for t0 in t:
281 if not (len(t0)>0 and qa.isquantity(t0)):
282 if halt:
283 self.msg("can't interpret '"+str(t0)+"' as a CASA quantity",priority="error")
284 return False
285 return True
287 ###########################################################
289 def isdirection(self,s,halt=True):
290 if type(s)==type([]):
291 t=s[0]
292 else:
293 t=s
294 try:
295 x=self.direction_splitter(s)
296 y=me.direction(x[0],x[1],x[2])
297 except:
298 if halt:
299 self.msg("can't interpret '"+str(s)+"' as a direction",priority="error")
300 return False
301 if not me.measure(y,y['refer']):
302 if halt:
303 self.msg("can't interpret '"+str(s)+"' as a direction",priority="error")
304 return False
305 return True
307 ###########################################################
309 def ismstp(self,s,halt=False):
310 try:
311 istp = False
312 # check if the ms is tp data or not.
313 tb.open(s+'/ANTENNA')
314 antname = tb.getcol('NAME')
315 tb.close()
316 if antname[0].find('TP') > -1:
317 istp = True
318 elif len(antname) == 1:
319 istp = True
320 else:
321 # need complete testing of UVW
322 tb.open(s)
323 uvw = tb.getcol("UVW")
324 tb.close()
325 if uvw.all() == 0.:
326 istp = True
327 except:
328 if halt:
329 self.msg("can't understand the file '"+str(s)+"'",priority="error")
330 return False
331 if not istp:
332 if halt:
333 self.msg("input file '"+str(s)+"' is not a totalpower ms",priority="error")
334 return False
335 return True
339 ###########################################################
340 # plot an image (optionally), and calculate its statistics
342 def statim(self,image,plot=True,incell=None,disprange=None,bar=True,showstats=True):
343 pix=self.cellsize(image) # cell positive by convention
344 pixarea=abs(qa.convert(pix[0],'arcsec')['value']*
345 qa.convert(pix[1],'arcsec')['value'])
346 ia.open(image)
347 imunit=ia.brightnessunit()
348 if imunit == 'Jy/beam':
349 bm=ia.restoringbeam()
350 if len(bm)>0:
351 toJyarcsec=1./pixarea
352 else:
353 toJyarcsec=1.
354 toJypix=toJyarcsec*pixarea
355 elif imunit == 'Jy/pixel':
356 toJyarcsec=1./pixarea
357 toJypix=1.
358 else:
359 self.msg("%s: unknown units" % image,origin="statim")
360 toJyarcsec=1.
361 toJypix=1.
362 stats=ia.statistics(robust=True,verbose=False,list=False)
363 im_min=stats['min']*toJypix
364 plarr=pl.zeros(1)
365 badim=False
366 if type(im_min)==type([]) or type(im_min)==type(plarr):
367 if len(im_min)<1:
368 badim=True
369 im_min=0.
370 im_max=stats['max']*toJypix
371 if type(im_max)==type([]) or type(im_max)==type(plarr):
372 if len(im_max)<1:
373 badim=True
374 im_max=1.
375 imsize=ia.shape()[0:2]
376 reg1=rg.box([0,0],[imsize[0]*.25,imsize[1]*.25])
377 stats=ia.statistics(region=reg1,verbose=False,list=False)
378 im_rms=stats['rms']*toJypix
379 if type(im_rms)==type([]) or type(im_rms)==type(plarr):
380 if len(im_rms)==0:
381 badim=True
382 im_rms=0.
383 data_array=ia.getchunk([-1,-1,1,1],[-1,-1,1,1],[1],[],True,True,False)
384 data_array=pl.array(data_array)
385 tdata_array=pl.transpose(data_array)
386 ttrans_array=tdata_array.tolist()
387 ttrans_array.reverse()
389 # get and apply mask
390 mask_array=ia.getchunk([-1,-1,1,1],[-1,-1,1,1],[1],[],True,True,True)
391 mask_array=pl.array(mask_array)
392 tmask_array=pl.transpose(mask_array)
393 tmtrans_array=tmask_array.tolist()
394 tmtrans_array.reverse()
395 ttrans_array=pl.array(ttrans_array)
396 z=pl.where(pl.array(tmtrans_array)==False)
397 ttrans_array[z[0],z[1]]=0.
399 if (plot):
400 pixsize=[qa.convert(pix[0],'arcsec')['value'],qa.convert(pix[1],'arcsec')['value']]
401 xextent=imsize[0]*abs(pixsize[0])*0.5
402 yextent=imsize[1]*abs(pixsize[1])*0.5
403 if self.verbose:
404 self.msg("plotting %fx%f\" im with %fx%f\" pix" %
405 (xextent,yextent,pixsize[0],pixsize[1]),origin="statim")
406 xextent=[xextent,-xextent]
407 yextent=[-yextent,yextent]
408 # remove top .5% of pixels:
409 nbin=200
410 if badim:
411 highvalue=im_max
412 lowvalue=im_min
413 else:
414 #imhist=ia.histograms(cumu=True,nbins=nbin,list=False)#['histout']
415 imhist=ia.histograms(cumu=True,nbins=nbin)#['histout']
416 ii=0
417 lowcounts=imhist['counts'][ii]
418 while imhist['counts'][ii]<0.005*lowcounts and ii<nbin:
419 ii=ii+1
420 lowvalue=imhist['values'][ii]
421 ii=nbin-1
422 highcounts=imhist['counts'][ii]
423 while imhist['counts'][ii]>0.995*highcounts and ii>0 and 0.995*highcounts>lowcounts:
424 ii=ii-1
425 highvalue=imhist['values'][ii]
426 if disprange != None:
427 if type(disprange)==type([]):
428 if len(disprange)>0:
429 highvalue=disprange[-1]
430 if len(disprange)>1:
431 lowvalue=disprange[0]
432 if len(disprange)>2:
433 throw("internal error disprange="+str(disprange)+" has too many elements")
434 else: # if passed an empty list [], return low.high
435 disprange.append(lowvalue)
436 disprange.append(highvalue)
437 else:
438 highvalue=disprange # assume if scalar passed its the max
440 if plot:
441 img=pl.imshow(ttrans_array,interpolation='bilinear',cmap=pl.cm.jet,extent=xextent+yextent,vmax=highvalue,vmin=lowvalue)
442 ax=pl.gca()
443 #l=ax.get_xticklabels()
444 #pl.setp(l,fontsize="x-small")
445 #l=ax.get_yticklabels()
446 #pl.setp(l,fontsize="x-small")
447 foo=image.split("/")
448 if len(foo)==1:
449 imagestrip=image
450 else:
451 imagestrip=foo[1]
452 pl.title(imagestrip,fontsize="x-small")
453 if showstats:
454 pl.text(0.05,0.95,"min=%7.1e\nmax=%7.1e\nRMS=%7.1e\n%s" % (im_min,im_max,im_rms,imunit),transform = ax.transAxes,bbox=dict(facecolor='white', alpha=0.7),size="x-small",verticalalignment="top")
455 if bar:
456 cb=pl.colorbar(pad=0)
457 cl = pl.getp(cb.ax,'yticklabels')
458 pl.setp(cl,fontsize='x-small')
459 ia.done()
460 return im_min,im_max,im_rms,imunit
468 ###########################################################
470 def calc_pointings2(self, spacing, size, maptype="hex", direction=None, relmargin=0.5, beam=0.):
471 """
472 If direction is a list, simply returns direction and the number of
473 pointings in it.
475 Otherwise, returns a hexagonally packed list of pointings separated by
476 spacing and fitting inside an area specified by direction and mapsize,
477 as well as the number of pointings. The hexagonal packing starts with a
478 horizontal row centered on direction, and the other rows alternate
479 being horizontally offset by a half spacing.
480 """
481 # make size 2-dimensional and ensure it is quantity
482 if type(size) != type([]):
483 size=[size,size]
484 if len(size) <2:
485 size=[size[0],size[0]]
486 self.isquantity(size)
488 # parse and check direction
489 if direction==None:
490 # if no direction is specified, use the object's direction
491 direction=self.direction
492 else:
493 # if one is specified, use it to set the object's direction
494 self.direction=direction
495 self.isdirection(direction)
497 # direction is always a list of strings (defined by .xml)
498 if type(direction)==type([]):
499 if len(direction) > 1:
500 if self.verbose: self.msg("you are inputing the precise pointings in 'direction' - if you want to calculate a mosaic, give a single direction string and set maptype",priority="warn")
501 return direction
502 else: direction=direction[0]
505 # haveing eliminated other options, we need to calculate:
506 epoch, centx, centy = self.direction_splitter()
508 pointings = []
510 shorttype=str.upper(maptype[0:3])
511# if not shorttype=="HEX":
512# self.msg("can't calculate map of maptype "+maptype,priority="error")
513 if shorttype == "HEX":
514 # this is hexagonal grid - Kana will add other types here
515 self.isquantity(spacing)
516 spacing = qa.quantity(spacing)
517 yspacing = qa.mul(0.866025404, spacing)
519 xsize=qa.quantity(size[0])
520 ysize=qa.quantity(size[1])
522 nrows = 1+ int(pl.floor(qa.convert(qa.div(ysize, yspacing), '')['value']
523 - 2.309401077 * relmargin))
525 availcols = 1 + qa.convert(qa.div(xsize, spacing),
526 '')['value'] - 2.0 * relmargin
527 ncols = int(pl.floor(availcols))
529 # By making the even rows shifted spacing/2 ahead, and possibly shorter,
530 # the top and bottom rows (nrows odd), are guaranteed to be short.
531 if availcols - ncols >= 0.5 and nrows>1: # O O O
532 evencols = ncols # O O O
533 ncolstomin = 0.5 * (ncols - 0.5)
534 else:
535 evencols = ncols - 1 # O O
536 ncolstomin = 0.5 * (ncols - 1) # O O O
537 pointings = []
539 # Start from the top because in the Southern hemisphere it sets first.
540 y = qa.add(centy, qa.mul(0.5 * (nrows - 1), yspacing))
541 for row in range(0, nrows): # xrange stops early.
542 xspacing = qa.mul(1.0 / pl.cos(qa.convert(y, 'rad')['value']),spacing)
543 ystr = qa.formxxx(y, format='dms',prec=5)
545 if row % 2: # Odd
546 xmin = qa.sub(centx, qa.mul(ncolstomin, xspacing))
547 stopcolp1 = ncols
548 else: # Even (including 0)
549 xmin = qa.sub(centx, qa.mul(ncolstomin - 0.5,
550 xspacing))
551 stopcolp1 = evencols
552 for col in range(0, stopcolp1): # xrange stops early.
553 x = qa.formxxx(qa.add(xmin, qa.mul(col, xspacing)),
554 format='hms',prec=5)
555 pointings.append("%s%s %s" % (epoch, x, ystr))
556 y = qa.sub(y, yspacing)
557 elif shorttype =="SQU":
558 # lattice gridding
559 self.isquantity(spacing)
560 spacing = qa.quantity(spacing)
561 yspacing = spacing
563 xsize=qa.quantity(size[0])
564 ysize=qa.quantity(size[1])
566 nrows = 1+ int(pl.floor(qa.convert(qa.div(ysize, yspacing), '')['value']
567 - 2.0 * relmargin))
569 availcols = 1 + qa.convert(qa.div(xsize, spacing), '')['value'] \
570 - 2.0 * relmargin
571 ncols = int(pl.floor(availcols))
574 ncolstomin = 0.5 * (ncols - 1) # O O O O
575 pointings = [] # O O O O
577 # Start from the top because in the Southern hemisphere it sets first.
578 y = qa.add(centy, qa.mul(0.5 * (nrows - 1), yspacing))
579 for row in range(0, nrows): # xrange stops early.
580 xspacing = qa.mul(1.0 / pl.cos(qa.convert(y, 'rad')['value']),spacing)
581 ystr = qa.formxxx(y, format='dms',prec=5)
583 xmin = qa.sub(centx, qa.mul(ncolstomin, xspacing))
584 stopcolp1 = ncols
586 for col in range(0, stopcolp1): # xrange stops early.
587 x = qa.formxxx(qa.add(xmin, qa.mul(col, xspacing)),
588 format='hms',prec=5)
589 pointings.append("%s%s %s" % (epoch, x, ystr))
590 y = qa.sub(y, yspacing)
591 if shorttype == "ALM":
592 # OT algorithm
593 self.isquantity(spacing)
594 spacing = qa.quantity(spacing)
595 xsize = qa.quantity(size[0])
596 ysize = qa.quantity(size[1])
598 spacing_asec=qa.convert(spacing,'arcsec')['value']
599 xsize_asec=qa.convert(xsize,'arcsec')['value']
600 ysize_asec=qa.convert(ysize,'arcsec')['value']
601 angle = 0. # deg
603 if str.upper(maptype[0:8]) == 'ALMA2012':
604 x,y = self.getTrianglePoints(xsize_asec, ysize_asec, angle, spacing_asec)
605 else:
606 if beam<=0:
607 beam=spacing_asec*pbcoeff*pl.sqrt(3) # ASSUMES ALMA default and arcsec
608 x,y = self.getTriangularTiling(xsize_asec, ysize_asec, angle, spacing_asec, beam)
610 pointings = []
611 nx=len(x)
612 for i in range(nx):
613 # Start from the top because in the Southern hemisphere it sets first.
614 y1=qa.add(centy, str(y[nx-i-1])+"arcsec")
615 ycos=pl.cos(qa.convert(y1,"rad")['value'])
616 ystr = qa.formxxx(y1, format='dms',prec=5)
617 xstr = qa.formxxx(qa.add(centx, str(x[nx-i-1]/ycos)+"arcsec"), format='hms',prec=5)
618 pointings.append("%s%s %s" % (epoch, xstr, ystr))
621 # if could not fit any pointings, then return single pointing
622 if(len(pointings)==0):
623 pointings.append(direction)
625 self.msg("using %i generated pointing(s)" % len(pointings),origin='calc_pointings')
626 self.pointings=pointings
627 return pointings
638 ###########################################################
640 def read_pointings(self, filename):
641 """
642 read pointing list from file containing epoch, ra, dec,
643 and scan time (optional,in sec).
644 Parameter:
645 filename: (str) the name of input file
647 The input file (ASCII) should contain at least 3 fields separated
648 by a space which specify positions with epoch, ra and dec (in dms
649 or hms).
650 The optional field, time, shoud be a list of decimal numbers
651 which specifies integration time at each position in second.
652 The lines which start with '#' is ignored and can be used
653 as comment lines.
655 Example of an input file:
656 #Epoch RA DEC TIME(optional)
657 J2000 23h59m28.10 -019d52m12.35 10.0
658 J2000 23h59m32.35 -019d52m12.35 10.0
659 J2000 23h59m36.61 -019d52m12.35 60.0
661 """
662 f=open(filename)
663 line= ' '
664 time=[]
665 pointings=[]
667 # add option of different epoch in a header line like read_antenna?
669 while (len(line)>0):
670 try:
671 line=f.readline()
672 if (not line.startswith('#')) and (not line.startswith("RA")) and (not line.startswith("--")):
673 ### it could be the new OT format
674 if (line.find('SEXAGESIMAL')>0):
675 splitline = line.split(',')
676 if len(splitline)>4:
677 time.append(float(splitline[4]))
678 else:
679 time.append(0.)
680 xstr = qa.formxxx(qa.quantity(splitline[0]), format='hms',prec=5)
681 de0=splitline[1]
682 # casa insists that strings with : are RA, so...
683 if de0.count(":")>0:
684 ystr = qa.formxxx(qa.div(qa.quantity(de0),15), format='dms',prec=5)
685 else:
686 ystr = qa.formxxx(qa.quantity(de0), format='dms',prec=5)
687 # ASSUME ICRS
688 pointings.append("ICRS %s %s" % (xstr,ystr))
689 ### ignoring line that has less than 3 elements
690 elif(len(line.split()) >2):
691 splitline=line.split()
692 epoch=splitline[0]
693 ra0=splitline[1]
694 de0=splitline[2]
695 if len(splitline)>3:
696 time.append(float(splitline[3]))
697 else:
698 time.append(0.)
699 xstr = qa.formxxx(qa.quantity(ra0), format='hms',prec=5)
700 ystr = qa.formxxx(qa.quantity(de0), format='dms',prec=5)
701 pointings.append("%s %s %s" % (epoch,xstr,ystr))
702 except:
703 break
704 f.close()
706 # need an error check here if zero valid pointings were read
707 if len(pointings) < 1:
708 s="No valid lines found in pointing file"
709 self.msg(s,priority="error")
710 self.msg("read in %i pointing(s) from file" % len(pointings),origin="read_pointings")
711 self.pointings=pointings
712 #self.direction=pointings
714 return len(pointings), pointings, time
722 ###########################################################
724 def write_pointings(self, filename,pointings,time=1.):
725 """
726 write pointing list to file containing epoch, ra, dec,
727 and scan time (optional,in sec).
729 Example of an output file:
730 #Epoch RA DEC TIME(optional)
731 J2000 23h59m28.10 -019d52m12.35 10.0
732 J2000 23h59m32.35 -019d52m12.35 10.0
733 J2000 23h59m36.61 -019d52m12.35 60.0
735 """
736 f=open(filename,'w')
737 f.write('#Epoch RA DEC TIME[sec]\n')
738 if type(pointings)!=type([]):
739 pointings=[pointings]
740 npos=len(pointings)
741 if type(time)!=type([]):
742 time=[time]
743 if len(time)==1:
744 time=list(time[0] for x in range(npos))
746 for i in range(npos):
747 #epoch,ra,dec=direction_splitter(pointings[i])
748 #xstr = qa.formxxx(qa.quantity(ra), format='hms')
749 #ystr = qa.formxxx(qa.quantity(dec), format='dms')
750 #line = "%s %s %s" % (epoch,xstr,ystr)
751 #self.isdirection(line) # extra check
752 #f.write(line+" "+str(time[i])+"\n")
753 f.write(pointings[i]+" "+str(time[i])+"\n")
755 f.close()
756 return
760 ###########################################################
762 def average_direction(self, directions=None):
763 # RI TODO make deal with list of measures as well as list of strings
764 """
765 Returns the average of directions as a string, and relative offsets
766 """
767 if directions==None:
768 directions=self.direction
769 epoch0, x, y = self.direction_splitter(directions[0])
770 i = 1
771 avgx = 0.0
772 avgy = 0.0
773 for drn in directions:
774 epoch, x, y = self.direction_splitter(drn)
775 # in principle direction_splitter returns directions in degrees,
776 # but can we be sure?
777 x=qa.convert(x,'deg')
778 y=qa.convert(y,'deg')
779 x = x['value']
780 y = y['value']
781 if epoch != epoch0: # Paranoia
782 casalog.post("[simutil] WARN: precession not handled by average_direction()",
783 'WARN')
784 x = self.wrapang(x, avgx, 360.0)
785 avgx += (x - avgx) / i
786 avgy += (y - avgy) / i
787 i += 1
788 offsets=pl.zeros([2,i-1])
789 i=0
790 cosdec=pl.cos(avgy*pl.pi/180.)
791 for drn in directions:
792 epoch, x, y = self.direction_splitter(drn)
793 x=qa.convert(x,'deg')
794 y=qa.convert(y,'deg')
795 x = x['value']
796 y = y['value']
797 x = self.wrapang(x, avgx, 360.0)
798 offsets[:,i]=[(x-avgx)*cosdec,y-avgy] # apply cosdec to make offsets on sky
799 i+=1
800 avgx = qa.toangle('%17.12fdeg' % avgx)
801 avgy = qa.toangle('%17.12fdeg' % avgy)
802 avgx = qa.formxxx(avgx, format='hms',prec=5)
803 avgy = qa.formxxx(avgy, format='dms',prec=5)
804 return "%s%s %s" % (epoch0, avgx, avgy), offsets
808 ###########################################################
810 def median_direction(self, directions=None):
811 # RI TODO make deal with list of measures as well as list of strings
812 """
813 Returns the median of directions as a string, and relative offsets
814 """
815 if directions==None:
816 directions=self.direction
817 epoch0, x, y = self.direction_splitter(directions[0])
818 i = 1
819 avgx = 0.0
820 avgy = 0.0
821 xx=[]
822 yy=[]
823 for drn in directions:
824 epoch, x, y = self.direction_splitter(drn)
825 # in principle direction_splitter returns directions in degrees,
826 # but can we be sure?
827 x=qa.convert(x,'deg')
828 y=qa.convert(y,'deg')
829 x = x['value']
830 y = y['value']
831 if epoch != epoch0: # Paranoia
832 casalog.post("[simutil] WARN: precession not handled by average_direction()",
833 'WARN')
834 x = self.wrapang(x, avgx, 360.0)
835 xx.append(x)
836 yy.append(y)
837 i += 1
838 avgx = pl.median(xx)
839 avgy = pl.median(yy)
840 offsets=pl.zeros([2,i-1])
841 i=0
842 cosdec=pl.cos(avgy*pl.pi/180.)
843 for drn in directions:
844 epoch, x, y = self.direction_splitter(drn)
845 x=qa.convert(x,'deg')
846 y=qa.convert(y,'deg')
847 x = x['value']
848 y = y['value']
849 x = self.wrapang(x, avgx, 360.0)
850 offsets[:,i]=[(x-avgx)*cosdec,y-avgy] # apply cosdec to make offsets on sky
851 i+=1
852 avgx = qa.toangle('%17.12fdeg' % avgx)
853 avgy = qa.toangle('%17.12fdeg' % avgy)
854 avgx = qa.formxxx(avgx, format='hms',prec=5)
855 avgy = qa.formxxx(avgy, format='dms',prec=5)
856 return "%s%s %s" % (epoch0, avgx, avgy), offsets
860 ###########################################################
862 def direction_splitter(self, direction=None):
863 """
864 Given a direction, return its epoch, x, and y parts. Epoch will be ''
865 if absent, or '%s ' % epoch if present. x and y will be angle qa's in
866 degrees.
867 """
868 if direction == None:
869 direction=self.direction
870 if type(direction) == type([]):
871 direction=self.average_direction(direction)[0]
872 dirl = direction.split()
873 if len(dirl) == 3:
874 epoch = dirl[0] + ' '
875 else:
876 epoch = ''
877 # x, y = map(qa.toangle, dirl[-2:])
878 x=qa.toangle(dirl[1])
879 # qa is stupid when it comes to dms vs hms, and assumes hms with colons and dms for periods.
880 decstr=dirl[2]
881 # parse with regex to get three numbers and reinstall them as dms
882 q=re.compile('([+-]?\d+).(\d+).(\d+\.?\d*)')
883 qq=q.match(decstr)
884 z=qq.groups()
885 decstr=z[0]+'d'
886 if len(z)>1:
887 decstr=decstr+z[1]+'m'
888 if len(z)>2:
889 decstr=decstr+z[2]+'s'
890 y=qa.toangle(decstr)
892 return epoch, qa.convert(x, 'deg'), qa.convert(y, 'deg')
895 ###########################################################
897 def dir_s2m(self, direction=None):
898 """
899 Given a direction as a string 'refcode lon lat', return it as qa measure.
900 """
901 if direction == None:
902 direction=self.direction
903 if type(direction) == type([]):
904 direction=self.average_direction(direction)[0]
905 dirl = direction.split()
906 if len(dirl) == 3:
907 refcode = dirl[0] + ' '
908 else:
909 refcode = 'J2000'
910 if self.verbose: self.msg("assuming J2000 for "+direction,origin="simutil.s2m")
911 x, y = map(qa.quantity, dirl[-2:])
912 if x['unit'] == '': x['unit']='deg'
913 if y['unit'] == '': y['unit']='deg'
914 return me.direction(refcode,qa.toangle(x),qa.toangle(y))
917 ###########################################################
919 def dir_m2s(self, dir):
920 """
921 Given a direction as a measure, return it as astring 'refcode lon lat'.
922 """
923 if dir['type'] != 'direction':
924 self.msg("converting direction measure",priority="error",origin="simutil.m2s")
925 return False
926 ystr = qa.formxxx(dir['m1'], format='dms',prec=5)
927 xstr = qa.formxxx(dir['m0'], format='hms',prec=5)
928 return "%s %s %s" % (dir['refer'], xstr, ystr)
930 ###########################################################
932 def wrapang(self, ang, target, period = 360.0):
933 """
934 Returns ang wrapped so that it is within +-period/2 of target.
935 """
936 dang = ang - target
937 period = pl.absolute(period)
938 halfperiod = 0.5 * period
939 if pl.absolute(dang) > halfperiod:
940 nwraps = pl.floor(0.5 + float(dang) / period)
941 ang -= nwraps * period
942 return ang
953 ###########################################################
954 #========================== tsys ==========================
956 def noisetemp(self, telescope=None, freq=None,
957 diam=None, epsilon=None):
959 """
960 Noise temperature and efficiencies for several telescopes:
961 ALMA, ACA, EVLA, VLA, and SMA
962 Input: telescope name, frequency as a quantity string "300GHz",
963 dish diameter (optional - knows diameters for arrays above)
964 epsilon = rms surface accuracy in microns (also optional -
965 this method contains the spec values for each telescope)
966 Output: eta_p phase efficieny (from Ruze formula),
967 eta_s spill (main beam) efficiency,
968 eta_b geometrical blockage efficiency,
969 eta_t taper efficiency,
970 eta_q correlator efficiency including quantization,
971 t_rx reciever temperature.
972 antenna efficiency = eta_p * eta_s * eta_b * eta_t
973 Notes: VLA correlator efficiency includes waveguide loss
974 EVLA correlator efficiency is probably optimistic at 0.88
975 """
977 if telescope==None: telescope=self.telescopename
978 # Noisetemp only knows about 6 observatories,
979 # none of which have measure.observatory dicts that
980 # contain lowercase characters so str.upper is harmless here
981 telescope=str.upper(telescope)
983 obs =['ALMASD','ALMA','ACA','EVLA','VLA','SMA']
984 d =[ 12. ,12. ,7. ,25. ,25. , 6. ]
985 ds =[ 0.75 ,0.75 ,0.75 ,0.364 ,0.364,0.35] # subreflector size for ACA?
986 eps =[ 25. ,25. ,20. ,300 ,300 ,15. ] # antenna surface accuracy
988 cq =[ 0.845, 0.845, 0.88, 0.79, 0.86, 0.88] # correlator eff
989 # SMA is from Wright http://astro.berkeley.edu/~wright/obsrms.py
990 # ALMA includes quantization eff of 0.96
991 # VLA includes additional waveguide loss from correlator loss of 0.809
992 # EVLA is probably optimistic
994 # things hardcoded in ALMA etimecalculator
995 # t_ground=270.
996 # t_cmb=2.73
997 # eta_q*eta_corr = 0.88*.961
998 # eta_ap = 0.72*eta_ruze
1000 if obs.count(telescope)>0:
1001 iobs=obs.index(telescope)
1002 else:
1003 if self.verbose: self.msg("I don't know much about "+telescope+" so I'm going to use ALMA specs")
1004 iobs=1 # ALMA is the default ;)
1006 if diam==None: diam=d[iobs]
1007 diam_subreflector=ds[iobs]
1008 if self.verbose: self.msg("subreflector diameter="+str(diam_subreflector),origin="noisetemp")
1010 # blockage efficiency.
1011 eta_b = 1.-(diam_subreflector/diam)**2
1013 # spillover efficiency.
1014 eta_s = 0.95 # these are ALMA values
1015 # taper efficiency.
1016 #eta_t = 0.86 # these are ALMA values
1017 eta_t = 0.819 # 20100914 OT value
1018 eta_t = 0.72
1020 # Ruze phase efficiency.
1021 if epsilon==None: epsilon = eps[iobs] # microns RMS
1022 if freq==None:
1023 freq_ghz=qa.convert(qa.quantity(self.centerfreq),'GHz')
1024 bw_ghz=qa.convert(qa.quantity(self.bandwidth),'GHz')
1025 else:
1026 freq_ghz=qa.convert(qa.quantity(freq),'GHz')
1027 eta_p = pl.exp(-(4.0*3.1415926535*epsilon*freq_ghz.get("value")/2.99792458e5)**2)
1028 if self.verbose: self.msg("ruze phase efficiency for surface accuracy of "+str(epsilon)+"um = " + str(eta_p) + " at "+str(freq),origin="noisetemp")
1030 # antenna efficiency
1031 # eta_a = eta_p*eta_s*eta_b*eta_t
1033 # correlator quantization efficiency.
1034 eta_q = cq[iobs]
1036 # Receiver radiation temperature in K.
1037 if telescope=='ALMA' or telescope=='ACA' or telescope=='ALMASD':
1038 # ALMA-40.00.00.00-001-A-SPE.pdf
1039 # http://www.eso.org/sci/facilities/alma/system/frontend/
1041 # lower limits
1042 #f0=[ 31, 67, 84, 125, 162, 211, 275, 385, 602, 787, 950]
1043 # go to higher band in gaps
1044 #f0=[ 31, 45, 84, 116, 162, 211, 275, 373, 500, 720, 950]
1045 # CASR-669
1046 f0=[ 35, 50, 84, 116, 162, 211, 275, 373, 500, 720, 950]
1047 # 80% spec
1048 #t0=[ 17, 30, 37, 51, 65, 83, 147, 196, 175, 230]
1049 # cycle 1 OT values 7/12
1050 #t0=[ 17, 30, 45, 51, 65, 55, 75, 196, 100, 230]
1051 # CAS-8802 B5 = 163-211GHz:
1052 #t0=[ 17, 30, 45, 51, 55, 55, 75, 196, 100, 230]
1053 # CAS-12387 match Cycle 7 OT and ASC
1054 #t0=[ 25, 30, 40, 42, 50, 50, 72, 135, 105, 230]
1055 # CASR-669
1056 t0=[ 28, 30, 40, 42, 50, 50, 72, 135, 105, 230]
1058 #flim=[31.3,950]
1059 # CASR-669
1060 flim=[35.0,950]
1061 if self.verbose: self.msg("using ALMA/ACA Rx specs",origin="noisetemp")
1062 else:
1063 if telescope=='EVLA':
1064 # 201009114 from rick perley:
1065 # f0=[1.5,3,6,10,15,23,33,45]
1066 t0=[10.,15,12,15,10,12,15,28]
1067 # limits
1068 f0=[1,2,4,8,12,18,26.5,40,50]
1070 flim=[0.8,50]
1071 if self.verbose: self.msg("using EVLA Rx specs",origin="noisetemp")
1072 else:
1073 if telescope=='VLA':
1074 # http://www.vla.nrao.edu/genpub/overview/
1075 # f0=[0.0735,0.32, 1.5, 4.75, 8.4, 14.9, 23, 45 ]
1076 # t0=[5000, 165, 56, 44, 34, 110, 110, 110]
1077 # exclude P band for now
1078 # f0=[0.32, 1.5, 4.75, 8.4, 14.9, 23, 45 ]
1079 # limits
1080 # http://www.vla.nrao.edu/genpub/overview/
1081 f0=[0.30,0.34,1.73,5,8.8,15.4,24,50]
1082 t0=[165, 56, 44, 34, 110, 110, 110]
1083 flim=[0.305,50]
1084 if self.verbose: self.msg("using old VLA Rx specs",origin="noisetemp")
1085 else:
1086 if telescope=='SMA':
1087 # f0=[212.,310.,383.,660.]
1088 # limits
1089 f0=[180,250,320,420,720]
1090 t0=[67, 116, 134, 500]
1091 flim=[180.,720]
1092 else:
1093 self.msg("I don't know about the "+
1094 telescope+" receivers, using 200K",
1095 priority="warn",origin="noisetemp")
1096 f0=[10,900]
1097 t0=[200,200]
1098 flim=[0,5000]
1100 obsfreq=freq_ghz.get("value")
1101 # z=pl.where(abs(obsfreq-pl.array(f0)) == min(abs(obsfreq-pl.array(f0))))
1102 # t_rx=t0[z[0]]
1104 if obsfreq<flim[0]:
1105 t_rx=t0[0]
1106 self.msg("observing freqency is lower than expected for "+
1107 telescope,priority="warn",origin="noise")
1108 self.msg("proceeding with extrapolated receiver temp="+
1109 str(t_rx),priority="warn",origin="noise")
1110 else:
1111 z=0
1112 while(f0[z]<obsfreq and z<len(t0)):
1113 z+=1
1114 t_rx=t0[z-1]
1116 if obsfreq>flim[1]:
1117 self.msg("observing freqency is higher than expected for "+
1118 telescope,priority="warn",origin="noise")
1119 self.msg("proceeding with extrapolated receiver temp="+
1120 str(t_rx),priority="warn",origin="noise")
1121 if obsfreq<=flim[1] and obsfreq>=flim[0]:
1122 self.msg("interpolated receiver temp="+str(t_rx),origin="noise")
1124 return eta_p, eta_s, eta_b, eta_t, eta_q, t_rx
1132 def sensitivity(self, freq, bandwidth, etime, elevation, pwv=None,
1133 telescope=None, diam=None, nant=None,
1134 antennalist=None,
1135 doimnoise=None,
1136 integration=None,debug=None,
1137 method="tsys-atm",tau0=None,t_sky=None):
1139 if qa.convert(elevation,"deg")['value']<3:
1140 self.msg("Elevation < ALMA limit of 3 deg",priority="error")
1141 return False
1143 tmpname="tmp"+str(os.getpid())
1144 i=0
1145 while i<10 and len(glob.glob(tmpname+"*"))>0:
1146 tmpname="tmp"+str(os.getpid())+str(i)
1147 i=i+1
1148 if i>=9:
1149 xx=glob.glob(tmpname+"*")
1150 for k in range(len(xx)):
1151 if os.path.isdir(xx[k]):
1152 ctsys.removetable(xx[k])
1153 else:
1154 os.remove(xx[k])
1156 msfile=tmpname+".ms"
1157 sm.open(msfile)
1159 rxtype=0 # 2SB
1160 if antennalist==None:
1161 if telescope==None:
1162 self.msg("Telescope name has not been set.",priority="error")
1163 return False
1164 if diam==None:
1165 self.msg("Antenna diameter has not been set.",priority="error")
1166 return False
1167 if nant==None:
1168 self.msg("Number of antennas has not been set.",priority="error")
1169 return False
1171 known, found_match, found_partial_match = False, False, False
1172 # check if telescope is known to measures tool
1173 # ensure case insensitivity - CAS-12753
1174 obslist_lower = [obs.lower() for obs in me.obslist()]
1175 if telescope.lower() in obslist_lower:
1176 found_match = True
1177 else:
1178 # e.g., aca.tp.cfg, where a substring of a known obsname (ALMASD) is known
1179 for obsname in obslist_lower:
1180 if obsname in telescope.lower():
1181 found_partial_match = True
1183 if found_match == True or found_partial_match == True:
1184 known = True
1186 if known == True:
1187 posobs = me.measure(me.observatory(telescope), 'WGS84')
1188 else:
1189 self.msg("Unknown telescope and no antenna list.",
1190 priority="error")
1191 return False
1193 obslat=qa.convert(posobs['m1'],'deg')['value']
1194 obslon=qa.convert(posobs['m0'],'deg')['value']
1195 obsalt=qa.convert(posobs['m2'],'m')['value']
1196 stnx,stny,stnz = self.locxyz2itrf(obslat,obslon,obsalt,0,0,0)
1197 antnames="A00"
1198 # use a single dish of diameter diam - must be set if antennalist
1199 # is not set.
1200 stnd=[diam]
1202 else:
1203 if str.upper(antennalist[0:4])=="ALMA":
1204 tail=antennalist[5:]
1205 if self.isquantity(tail,halt=False):
1206 resl=qa.convert(tail,"arcsec")['value']
1207 repodir=os.getenv("CASAPATH").split(' ')[0]+"/data/alma/simmos/"
1208 if os.path.exists(repodir):
1209 confnum=(2.867-pl.log10(resl*1000*qa.convert(freq,"GHz")['value']/672.))/0.0721
1210 confnum=max(1,min(28,confnum))
1211 conf=str(int(round(confnum)))
1212 if len(conf)<2: conf='0'+conf
1213 antennalist=repodir+"alma.out"+conf+".cfg"
1214 self.msg("converted resolution to antennalist "+antennalist)
1216 repodir = os.getenv("CASAPATH").split(' ')[0] + "/data/alma/simmos/"
1217 if not os.path.exists(antennalist) and \
1218 os.path.exists(repodir+antennalist):
1219 antennalist = repodir + antennalist
1220 if os.path.exists(antennalist):
1221 stnx, stny, stnz, stnd, padnames, antnames, telescope, posobs = self.readantenna(antennalist)
1222 else:
1223 self.msg("antennalist "+antennalist+" not found",priority="error")
1224 return False
1226 nant=len(stnx)
1227 # diam is only used as a test below, not quantitatively
1228 diam = pl.average(stnd)
1231 if (telescope==None or diam==None):
1232 self.msg("Telescope name or antenna diameter have not been set.",priority="error")
1233 return False
1235 # copied from task_simdata:
1237 self.setcfg(sm, telescope, stnx, stny, stnz, stnd, padnames, posobs)
1239 model_nchan=1
1240 # RI TODO isquantity checks
1241 model_width=qa.quantity(bandwidth) # note: ATM uses band center
1243 # start is center of first channel. for nch=1, that equals center
1244 model_start=qa.quantity(freq)
1246 stokes, feeds = self.polsettings(telescope)
1247 sm.setspwindow(spwname="band1", freq=qa.tos(model_start),
1248 deltafreq=qa.tos(model_width),
1249 freqresolution=qa.tos(model_width),
1250 nchannels=model_nchan, stokes=stokes)
1251 sm.setfeed(mode=feeds, pol=[''])
1253 sm.setlimits(shadowlimit=0.01, elevationlimit='3deg')
1254 sm.setauto(0.0)
1256 obslat=qa.convert(posobs['m1'],'deg')
1257 dec=qa.add(obslat, qa.add(qa.quantity("90deg"),qa.mul(elevation,-1)))
1259 sm.setfield(sourcename="src1",
1260 sourcedirection="J2000 00:00:00.00 "+qa.angle(dec)[0],
1261 calcode="OBJ", distance='0m')
1262 reftime = me.epoch('TAI', "2012/01/01/00:00:00")
1263 if integration==None:
1264 integration=qa.mul(etime,0.01)
1265 self.msg("observing for "+qa.tos(etime)+" with integration="+qa.tos(integration))
1266 sm.settimes(integrationtime=integration, usehourangle=True,
1267 referencetime=reftime)
1269 sm.observe(sourcename="src1", spwname="band1",
1270 starttime=qa.quantity(0, "s"),
1271 stoptime=qa.quantity(etime));
1273 sm.setdata()
1274 sm.setvp()
1276 eta_p, eta_s, eta_b, eta_t, eta_q, t_rx = self.noisetemp(telescope=telescope,freq=freq)
1277 eta_a = eta_p * eta_s * eta_b * eta_t
1278 if self.verbose:
1279 self.msg('antenna efficiency = ' + str(eta_a),origin="noise")
1280 self.msg('spillover efficiency = ' + str(eta_s),origin="noise")
1281 self.msg('correlator efficiency = ' + str(eta_q),origin="noise")
1283 if pwv==None:
1284 # RI TODO choose based on freq octile
1285 pwv=2.0
1287 # things hardcoded in ALMA etimecalculator, & defaults in simulator.xml
1288 t_ground=270.
1289 t_cmb=2.725
1290 # eta_q = 0.88
1291 # eta_a = 0.95*0.8*eta_s
1293 if telescope=='ALMA' and (qa.convert(freq,"GHz")['value'])>600.:
1294 rxtype=1 # DSB
1296 if method=="tsys-atm":
1297 sm.setnoise(spillefficiency=eta_s,correfficiency=eta_q,
1298 antefficiency=eta_a,trx=t_rx,
1299 tground=t_ground,tcmb=t_cmb,pwv=str(pwv)+"mm",
1300 mode="tsys-atm",table=tmpname,rxtype=rxtype)
1301 else:
1302 if method=="tsys-manual":
1303 if not t_sky:
1304 t_sky=200.
1305 self.msg("Warning: Sky brightness temp not set, using 200K",origin="noise",priority="warn")
1306 sm.setnoise(spillefficiency=eta_s,correfficiency=eta_q,
1307 antefficiency=eta_a,trx=t_rx,tatmos=t_sky,
1308 tground=t_ground,tcmb=t_cmb,tau=tau0,
1309 mode="tsys-manual",table=tmpname,rxtype=rxtype)
1310 else:
1311 self.msg("Unknown calculation method "+method,priority="error")
1312 return False
1314 if doimnoise:
1315 sm.corrupt()
1317 sm.done()
1320 if doimnoise:
1321 cellsize=qa.quantity(6.e3/250./qa.convert(model_start,"GHz")["value"],"arcsec") # need better cell determination - 250m?!
1322 cellsize=[cellsize,cellsize]
1323 # very light clean - its an empty image!
1324 self.imclean(tmpname+".ms",tmpname,
1325 "csclean",cellsize,[128,128],
1326 "J2000 00:00:00.00 "+qa.angle(dec)[0],
1327 False,100,"0.01mJy","natural",[],True,"I")
1329 ia.open(tmpname+".image")
1330 stats= ia.statistics(robust=True, verbose=False,list=False)
1331 ia.done()
1332 imnoise=(stats["rms"][0])
1333 else:
1334 imnoise=0.
1336 nint = qa.convert(etime,'s')['value'] / qa.convert(integration,'s')['value']
1337 nbase= 0.5*nant*(nant-1)
1339 if os.path.exists(tmpname+".T.cal"):
1340 tb.open(tmpname+".T.cal")
1341 gain=tb.getcol("CPARAM")
1342 flag=tb.getcol("FLAG")
1343 # RI TODO average instead of first?
1344 tb.done()
1345 # gain is per ANT so square for per baseline;
1346 # pick a gain from about the middle of the track
1347 # needs to be unflagged!
1348 z=pl.where(flag[0][0]==False)[0]
1349 nunflagged=len(z)
1350# noiseperbase=1./(gain[0][0][0.5*nint*nant].real)**2
1351 noiseperbase=1./(gain[0][0][z[nunflagged//2]].real)**2
1352 else:
1353 noiseperbase=0.
1355 theoreticalnoise=noiseperbase/pl.sqrt(nint*nbase*2) # assume 2-poln
1357 if debug!=True:
1358 xx=glob.glob(tmpname+"*")
1359 for k in range(len(xx)):
1360 if os.path.isdir(xx[k]):
1361 ctsys.removetable(xx[k])
1362 else:
1363 os.remove(xx[k])
1365 if doimnoise:
1366 return theoreticalnoise , imnoise
1367 else:
1368 return theoreticalnoise
1371 def setcfg(self, mysm, telescope, x, y, z, diam,
1372 padnames, posobs, mounttype=None):
1373 """
1374 Sets the antenna positions for the mysm sm instance, which should have
1375 already opened an MS, given
1376 telescope - telescope name
1377 x - array of X positions, i.e. stnx from readantenna
1378 y - array of Y positions, i.e. stny from readantenna
1379 z - array of Z positions, i.e. stnz from readantenna
1380 diam - numpy.array of antenna diameters, i.e. from readantenna
1381 padnames - list of pad names
1382 antnames - list of antenna names
1383 posobs - The observatory position as a measure.
1385 Optional:
1386 mounttype (Defaults to a guess based on telescope.)
1388 Returns the mounttype that it uses.
1390 Closes mysm and throws a ValueError if it can't set the configuration.
1391 """
1392 if not mounttype:
1393 mounttype = 'alt-az'
1394 if telescope.upper() in ('ASKAP', # Effectively, if not BIZARRE.
1395 'DRAO',
1396 'WSRT'):
1397 mounttype = 'EQUATORIAL'
1398 elif telescope.upper() in ('LOFAR', 'LWA', 'MOST'):
1399 # And other long wavelength arrays...like that one in WA that
1400 # has 3 orthogonal feeds so they can go to the horizon.
1401 #
1402 # The SKA should go here too once people accept
1403 # that it will have to be correlated as stations.
1404 mounttype = 'BIZARRE' # Ideally it would be 'FIXED' or
1405 # 'GROUND'.
1407 if not mysm.setconfig(telescopename=telescope, x=x, y=y, z=z,
1408 dishdiameter=diam.tolist(),
1409 mount=[mounttype], padname=padnames,
1410 coordsystem='global', referencelocation=posobs):
1411 mysm.close()
1412 raise ValueError("Error setting the configuration")
1413 return mounttype
1418 ###########################################################
1419 #===================== ephemeris ==========================
1422 def ephemeris(self, date, usehourangle=True, direction=None, telescope=None, ms=None, cofa=None):
1424 if direction==None: direction=self.direction
1425 if telescope==None: telescope=self.telescopename
1427 # right now, simdata centers at the transit; when that changes,
1428 # or when that becomes optional, then that option needs to be
1429 # stored in the simutil object and used here, to either
1430 # center the plot at transit or not.
1431 #
1432 # direction="J2000 18h00m00.03s -45d59m59.6s"
1433 # refdate="2012/06/21/03:25:00"
1435 # useHourAngle_p means simulate at transit
1436 # TODO: put in reftime parameter, parse 2012/05/21, 2012/05/21/transit,
1437 # and 2012/05/21/22:05:00 separately.
1439 ds=self.direction_splitter(direction) # if list, returns average
1440 src=me.direction(ds[0],ds[1],ds[2])
1442 me.done()
1443 posobs=me.observatory(telescope)
1444 if len(posobs)<=0:
1445 if cofa==None:
1446 self.msg("simutil::ephemeris needs either a known observatory or an explicitly specified cofa measure",priority="error")
1447 posobs=cofa
1448 me.doframe(posobs)
1450 time=me.epoch('TAI',date)
1451 me.doframe(time)
1453 # what is HA of source at refdate?
1454 offset_ha=qa.convert((me.measure(src,'hadec'))['m0'],'h')
1455 peak=me.epoch("TAI",qa.add(date,qa.mul(-1,offset_ha)))
1456 peaktime_float=peak['m0']['value']
1457 if usehourangle:
1458 # offset the reftime to be at transit:
1459 time=peak
1460 me.doframe(time)
1462 reftime_float=time['m0']['value']
1463 reftime_floor=pl.floor(time['m0']['value'])
1464 refdate_str=qa.time(qa.totime(str(reftime_floor)+'d'),form='dmy')[0]
1466 timeinc='15min' # for plotting
1467 timeinc=qa.convert(qa.time(timeinc)[0],'d')['value']
1468 ntime=int(1./timeinc)
1470 # check for circumpolar:
1471 rset = me.riseset(src)
1472 rise = rset['rise']
1473 if rise == 'above':
1474 rise = time
1475 rise['m0']['value'] = rise['m0']['value'] - 0.5
1476 settime = time
1477 settime['m0']['value'] = settime['m0']['value'] + 0.5
1478 elif rise == 'below':
1479 raise ValueError(direction + ' is not visible from ' + telescope)
1480 else:
1481 settime = rset['set']
1482 rise = me.measure(rise['utc'],'tai')
1483 settime = me.measure(settime['utc'],'tai')
1485 # where to start plotting?
1486 offset=-0.5
1487 # this assumes that the values can be compared directly - already assumed in code above
1488 if settime['m0']['value'] < time['m0']['value']: offset -= 0.5
1489 if rise['m0']['value'] > time['m0']['value']: offset+=0.5
1490 time['m0']['value']+=offset
1492 times=[]
1493 az=[]
1494 el=[]
1496 for i in range(ntime):
1497 times.append(time['m0']['value'])
1498 me.doframe(time)
1499 azel=me.measure(src,'azel')
1500 az.append(qa.convert(azel['m0'],'deg')['value'])
1501 el.append(qa.convert(azel['m1'],'deg')['value'])
1502 time['m0']['value']+=timeinc
1504# self.msg(" ref="+date,origin='ephemeris')
1505# self.msg("rise="+qa.time(rise['m0'],form='dmy')[0],origin='ephemeris')
1506# self.msg(" set="+qa.time(settime['m0'],form='dmy')[0],origin='ephemeris')
1508 pl.plot((pl.array(times)-reftime_floor)*24,el)
1509# peak=(rise['m0']['value']+settime['m0']['value'])/2
1510# self.msg("peak="+qa.time('%fd' % peak,form='dmy'),origin='ephemeris')
1511 self.msg("peak="+qa.time('%fd' % reftime_float,form='dmy')[0],origin='ephemeris')
1513 relpeak=peaktime_float-reftime_floor
1514 pl.plot(pl.array([1,1])*24*relpeak,[0,90])
1516 # if theres an ms, figure out the entire range of observation
1517 if ms:
1518 tb.open(ms+"/OBSERVATION")
1519 timerange=tb.getcol("TIME_RANGE")
1520 tb.done()
1521 obsstart=min(timerange.flat)
1522 obsend=max(timerange.flat)
1523 relstart=me.epoch("UTC",str(obsstart)+"s")['m0']['value']-reftime_floor
1524 relend=me.epoch("UTC",str(obsend)+"s")['m0']['value']-reftime_floor
1525 pl.plot([relstart*24,relend*24],[89,89],'r',linewidth=3)
1526 else:
1527 if self.totaltime>0:
1528 etimeh=qa.convert(self.totaltime,'h')['value']
1529 pl.plot(pl.array([-0.5,0.5])*etimeh+relpeak*24,[80,80],'r')
1531 pl.xlabel("hours relative to "+refdate_str,fontsize='x-small')
1532 pl.ylabel("elevation",fontsize='x-small')
1533 ax=pl.gca()
1534 l=ax.get_xticklabels()
1535 pl.setp(l,fontsize="x-small")
1536 l=ax.get_yticklabels()
1537 pl.setp(l,fontsize="x-small")
1540 pl.ylim([0,90])
1541 pl.xlim(pl.array([-12,12])+24*(reftime_float-reftime_floor))
1561 ###########################################################
1562 #==========================================================
1564 def readantenna(self, antab=None):
1565 """
1566 Helper function to read antenna configuration file; example:
1567 # observatory=ALMA
1568 # COFA=-67.75,-23.02
1569 # coordsys=LOC (local tangent plane)
1570 # uid___A002_Xdb6217_X55ec_target.ms
1571 # x y z diam station ant
1572 -5.850273514 -125.9985379 -1.590364043 12. A058 DA41
1573 -19.90369337 52.82680653 -1.892119601 12. A023 DA42
1574 13.45860758 -5.790196849 -2.087805181 12. A035 DA43
1575 5.606192499 7.646657746 -2.087775605 12. A001 DA44
1576 24.10057423 -25.95933768 -2.08466565 12. A036 DA45
1577 lines beginning with "#" will be interpreted as header key=value
1578 pairs if they contain "=", and as comments otherwise
1579 for the observatory name, one can check the known observatories list
1580 me.obslist
1581 if an unknown observatory is specified, then one must either
1582 use absolute positions (coordsys XYZ,UTM), or
1583 specify COFA= lon,lat
1584 coordsys can be XYZ=earth-centered,
1585 UTM=easting,northing,altitude, or LOC=xoffset,yoffset,height
1587 returns: earth-centered x,y,z, diameter, pad_name, antenna_name, observatory_name, observatory_measure_dictionary
1589 """
1590 # all coord systems should have x,y,z,diam,pid, where xyz varies
1591 params={} # to store key=value "header" parameters
1592 inx=[]
1593 iny=[] # to store antenna positions
1594 inz=[]
1595 ind=[] # antenna diameter
1596 pid=[] # pad id
1597 aid=[] # antenna names
1598 nant=0 # counter for input lines that meet format requirements
1600 self.msg("Reading antenna positions from '%s'" % antab,
1601 origin="readantenna")
1602 with open(antab) as f:
1603 try: # to read the file
1604 for line in f:
1605 cols = line.split()
1606 if len(cols) == 0:
1607 pass # ignore empty rows
1608 if line.startswith('#'): # line is header
1609 paramlist = line[1:].split('=')
1610 if len(paramlist)==2: # key=value pair found
1611 # remove extra octothorpes then clean up and assign
1612 key = paramlist[0].replace('#','').strip()
1613 val = paramlist[1].replace('#','').strip()
1614 try:
1615 params[key]=val
1616 except TypeError:
1617 # params = None somehow? Legacy...
1618 params={key:val}
1619 else:
1620 # no key=value pair found, so ignore
1621 pass
1622 else: # otherwise octothorpe starts comments
1623 # take only what's in front of the first one
1624 line = line.partition('#')[0]
1625 cols = line.split()
1626 # now check for data
1627 if len(cols)>3:
1628 # assign first four columns, assuming order
1629 inx.append(float(cols[0]))
1630 iny.append(float(cols[1]))
1631 inz.append(float(cols[2]))
1632 ind.append(float(cols[3]))
1633 if len(cols) > 5:
1634 # Found unique pad and antenna names
1635 pid.append(cols[4])
1636 aid.append(cols[5])
1637 elif len(cols) > 4:
1638 # Found pad names, but not antenna names, so dupl.
1639 pid.append(cols[4])
1640 aid.append(cols[4])
1641 else:
1642 # No antenna or pad names found in antab so default
1643 pid.append('A%02d'%nant)
1644 aid.append('A%02d'%nant)
1645 nant+=1
1646 except IOError:
1647 self.msg("Could not read file: '{}'".format(antab),
1648 origin='readantenna', priority='error')
1649 except ValueError:
1650 self.msg("Could not read file: '{}'".format(antab),
1651 origin='readantenna', priority='error')
1653 if "coordsys" not in params:
1654 self.msg("Must specify XYZ, UTM or LOC coordinate system"\
1655 " in antenna file header",
1656 origin="readantenna",priority="error")
1657 return -1
1658 else:
1659 self.coordsys=params["coordsys"]
1661 if "observatory" in params:
1662 self.telescopename=params["observatory"]
1663 else:
1664 self.telescopename="SIMULATED"
1665 if self.verbose:
1666 self.msg("Using observatory= %s" % self.telescopename,
1667 origin="readantenna")
1669 known, found_match, found_partial_match = False, False, False
1670 # case insensitive check if telescopename is known
1671 obslist_lower = [obs.lower() for obs in me.obslist()]
1672 if self.telescopename.lower() in obslist_lower:
1673 found_match = True
1674 else:
1675 # e.g., aca.tp.cfg, where a substring of a known obsname (ALMASD) is known
1676 for obsname in obslist_lower:
1677 if obsname in self.telescopename.lower():
1678 found_partial_match = True
1680 if found_match == True or found_partial_match == True:
1681 t = self.telescopename
1682 known = True
1683 posobs=me.measure(me.observatory(t),'WGS84')
1685 if "COFA" in params:
1686 if known:
1687 self.msg("antenna config file specifies COFA for a known observatory "+
1688 self.telescopename+", overriding with specified COFA.",priority="warn")
1689 obs_latlon=params["COFA"].split(",")
1690 cofa_lon=float(obs_latlon[0])
1691 cofa_lat=float(obs_latlon[1])
1692 cofa_alt=0.
1693 posobs=me.position("WGS84",qa.quantity(cofa_lon,"deg"),qa.quantity(cofa_lat,"deg"),qa.quantity(cofa_alt,"m"))
1694 elif not known:
1695 if params["coordsys"].upper()[0:3]=="LOC":
1696 self.msg("To use local coordinates in the antenna position file, you must either use a known observatory name, or provide the COFA explicitly",priority="error")
1697 return -1
1698 else:
1699 # we have absolute coords, so can create the posobs from their
1700 # average at the end
1701 posobs={}
1705 if (params["coordsys"].upper()=="XYZ"):
1706 ### earth-centered XYZ i.e. ITRF in casa
1707 stnx=inx
1708 stny=iny
1709 stnz=inz
1710 else:
1711 stnx=[]
1712 stny=[]
1713 stnz=[]
1714 if (params["coordsys"].upper()=="UTM"):
1715 ### expect easting, northing, altitude in m
1716 self.msg("Antenna locations in UTM; will read "\
1717 "from file easting, northing, elevation in m",
1718 origin="readantenna")
1719 if "zone" in params:
1720 zone=params["zone"]
1721 else:
1722 self.msg("You must specify zone=NN in your antenna file",
1723 origin="readantenna",priority="error")
1724 return -1
1725 if "datum" in params:
1726 datum=params["datum"]
1727 else:
1728 self.msg("You must specify datum in your antenna file",
1729 origin="readantenna",priority="error")
1730 return -1
1731 if "hemisphere" in params:
1732 nors=params["hemisphere"]
1733 nors=nors[0].upper()
1734 else:
1735 self.msg("You must specify hemisphere=N|S in your antenna file",origin="readantenna",priority="error")
1736 return -1
1738 vsave=self.verbose
1739 for i in range(len(inx)):
1740 x,y,z = self.utm2xyz(inx[i],iny[i],inz[i],int(zone),datum,nors)
1741 if i==1:
1742 self.verbose=False
1743 stnx.append(x)
1744 stny.append(y)
1745 stnz.append(z)
1746 self.verbose=vsave
1747 else:
1748 if (params["coordsys"].upper()[0:3]=="LOC"):
1749 obslat=qa.convert(posobs['m1'],'deg')['value']
1750 obslon=qa.convert(posobs['m0'],'deg')['value']
1751 obsalt=qa.convert(posobs['m2'],'m')['value']
1752 if self.verbose:
1753 self.msg("converting local tangent plane coordinates"\
1754 "to ITRF using observatory position"\
1755 "= %f %f " % (obslat,obslon),
1756 origin="readantenna")
1757 for i in range(len(inx)):
1758 x,y,z = self.locxyz2itrf(obslat,obslon,obsalt,
1759 inx[i],iny[i],inz[i])
1760 stnx.append(x)
1761 stny.append(y)
1762 stnz.append(z)
1764 if len(posobs)<=0:
1765 # we had absolute coords but unknown telescope and no COFA
1766 cofa_x=pl.average(x)
1767 cofa_y=pl.average(y)
1768 cofa_z=pl.average(z)
1769 cofa_lat,cofa_lon,cofa_alt=self.xyz2long(cofa_x,cofa_y,cofa_z,
1770 'WGS84')
1771 posobs=me.position("WGS84",qa.quantity(cofa_lon,"rad"),
1772 qa.quantity(cofa_lat,"rad"),
1773 qa.quantity(cofa_alt,"m"))
1775 return (stnx, stny, stnz, pl.array(ind), pid, aid, self.telescopename, posobs)
1794 ###########################################################
1795 #==================== geodesy =============================
1798 def tmgeod(self,n,e,eps,cm,fe,sf,so,r,v0,v2,v4,v6,fn,er,esq):
1799 """
1800 Transverse Mercator Projection
1801 conversion of grid coords n,e to geodetic coords
1802 revised subroutine of t. vincenty feb. 25, 1985
1803 orig. source: https://www.ngs.noaa.gov/TOOLS/program_descriptions.html
1804 converted from Fortran R Indebetouw Jan 2009
1805 ********** symbols and definitions ***********************
1806 latitude positive north, longitude positive west.
1807 all angles are in radian measure.
1809 input:
1810 n,e are northing and easting coordinates respectively
1811 er is the semi-major axis of the ellipsoid
1812 esq is the square of the 1st eccentricity
1813 cm is the central meridian of the projection zone
1814 fe is the false easting value at the cm
1815 eps is the square of the 2nd eccentricity
1816 sf is the scale factor at the cm
1817 so is the meridional distance (times the sf) from the
1818 equator to southernmost parallel of lat. for the zone
1819 r is the radius of the rectifying sphere
1820 u0,u2,u4,u6,v0,v2,v4,v6 are precomputed constants for
1821 determination of meridianal dist. from latitude
1822 output:
1823 lat,lon are lat. and long. respectively
1824 conv is convergence
1825 kp point scale factor
1827 the formula used in this subroutine gives geodetic accuracy
1828 within zones of 7 degrees in east-west extent. within state
1829 transverse mercator projection zones, several minor terms of
1830 the equations may be omitted (see a separate ngs publication).
1831 if programmed in full, the subroutine can be used for
1832 computations in surveys extending over two zones.
1833 """
1835 om=(n-fn+so)/(r*sf) # (northing - flag_north + distance_from_equator)
1836 cosom=pl.cos(om)
1837 foot=om+pl.sin(om)*cosom*(v0+v2*cosom*cosom+v4*cosom**4+v6*cosom**6)
1838 sinf=pl.sin(foot)
1839 cosf=pl.cos(foot)
1840 tn=sinf/cosf
1841 ts=tn*tn
1842 ets=eps*cosf*cosf
1843 rn=er*sf/pl.sqrt(1.-esq*sinf*sinf)
1844 q=(e-fe)/rn
1845 qs=q*q
1846 b2=-tn*(1.+ets)/2.
1847 b4=-(5.+3.*ts+ets*(1.-9.*ts)-4.*ets*ets)/12.
1848 b6=(61.+45.*ts*(2.+ts)+ets*(46.-252.*ts-60.*ts*ts))/360.
1849 b1=1.
1850 b3=-(1.+ts+ts+ets)/6.
1851 b5=(5.+ts*(28.+24.*ts)+ets*(6.+8.*ts))/120.
1852 b7=-(61.+662.*ts+1320.*ts*ts+720.*ts**3)/5040.
1853 lat=foot+b2*qs*(1.+qs*(b4+b6*qs))
1854 l=b1*q*(1.+qs*(b3+qs*(b5+b7*qs)))
1855 lon=-l/cosf+cm
1857 # compute scale factor
1858 fi=lat
1859 lam = lon
1860 sinfi=pl.sin(fi)
1861 cosfi=pl.cos(fi)
1862 l1=(lam-cm)*cosfi
1863 ls=l1*l1
1865 tn=sinfi/cosfi
1866 ts=tn*tn
1868 # convergence
1869 c1=-tn
1870 c3=(1.+3.*ets+2.*ets**2)/3.
1871 c5=(2.-ts)/15.
1872 conv=c1*l1*(1.+ls*(c3+c5*ls))
1874 # point scale factor
1875 f2=(1.+ets)/2.
1876 f4=(5.-4.*ts+ets*( 9.-24.*ts))/12.
1877 kp=sf*(1.+f2*ls*(1.+f4*ls))
1879 return lat,lon,conv,kp
1884 def tconpc(self,sf,orlim,er,esq,rf):
1886 """
1887 transverse mercator projection ***
1888 conversion of grid coords to geodetic coords
1889 revised subroutine of t. vincenty feb. 25, 1985
1890 converted from fortran r. indebetouw jan 2009
1891 ********** symbols and definitions ***********************
1892 input:
1893 rf is the reciprocal flattening of ellipsoid
1894 esq = e squared (eccentricity?)
1895 er is the semi-major axis for grs-80
1896 sf is the scale factor at the cm
1897 orlim is the southernmost parallel of latitude for which the
1898 northing coord is zero at the cm
1900 output:
1901 eps
1902 so is the meridional distance (times the sf) from the
1903 equator to southernmost parallel of lat. for the zone
1904 r is the radius of the rectifying sphere
1905 u0,u2,u4,u6,v0,v2,v4,v6 are precomputed constants for
1906 determination of meridional dist. from latitude
1907 **********************************************************
1908 """
1910 f=1./rf
1911 eps=esq/(1.-esq)
1912 pr=(1.-f)*er
1913 en=(er-pr)/(er+pr)
1914 en2=en*en
1915 en3=en*en*en
1916 en4=en2*en2
1918 c2=-3.*en/2.+9.*en3/16.
1919 c4=15.*en2/16.-15.*en4/32.
1920 c6=-35.*en3/48.
1921 c8=315.*en4/512.
1922 u0=2.*(c2-2.*c4+3.*c6-4.*c8)
1923 u2=8.*(c4-4.*c6+10.*c8)
1924 u4=32.*(c6-6.*c8)
1925 u6=128.*c8
1927 c2=3.*en/2.-27.*en3/32.
1928 c4=21.*en2/16.-55.*en4/32.
1929 c6=151.*en3/96.
1930 c8=1097.*en4/512.
1931 v0=2.*(c2-2.*c4+3.*c6-4.*c8)
1932 v2=8.*(c4-4.*c6+10.*c8)
1933 v4=32.*(c6-6.*c8)
1934 v6=128.*c8
1936 r=er*(1.-en)*(1.-en*en)*(1.+2.25*en*en+(225./64.)*en4)
1937 cosor=pl.cos(orlim)
1938 omo=orlim+pl.sin(orlim)*cosor*(u0+u2*cosor*cosor+u4*cosor**4+u6*cosor**6)
1939 so=sf*r*omo
1941 return eps,r,so,v0,v2,v4,v6
1947 def getdatum(self,datumcode,verbose=False):
1948 """
1949 local datums and ellipsoids;
1950 input: datam code e.g. 'WGS84','SAM56'
1951 output:
1952 dx, dy, dz [m] - offsets from ITRF x,y,z i.e. one would take local earth-centered xyz coordinates, add dx,dy,dz to get wgs84 earth-centered
1954 er = equatorial radius of the ellipsoid (semi-major axis) [m]
1955 rf = reciprocal of flatting of the ellipsoid
1956 """
1957 # set equatorial radius and inverse flattening
1959 ellipsoids={'AW':[6377563.396,299.3249647 ,'Airy 1830' ],
1960 'AM':[6377340.189,299.3249647 ,'Modified Airy' ],
1961 'AN':[6378160.0 ,298.25 ,'Australian National' ],
1962 'BR':[6377397.155,299.1528128 ,'Bessel 1841' ],
1963 'CC':[6378206.4 ,294.9786982 ,'Clarke 1866' ],
1964 'CD':[6378249.145,293.465 ,'Clarke 1880' ],
1965 'EA':[6377276.345,300.8017 ,'Everest (India 1830)' ],
1966 'RF':[6378137.0 ,298.257222101,'Geodetic Reference System 1980'],
1967 'HE':[6378200.0 ,298.30 ,'Helmert 1906' ],
1968 'HO':[6378270.0 ,297.00 ,'Hough 1960' ],
1969 'IN':[6378388.0 ,297.00 ,'International 1924' ],
1970 'SA':[6378160.0 ,298.25 ,'South American 1969' ],
1971 'WD':[6378135.0 ,298.26 ,'World Geodetic System 1972' ],
1972 'WE':[6378137.0 ,298.257223563,'World Geodetic System 1984' ]}
1974 datums={
1975 'AGD66' :[-133, -48, 148,'AN','Australian Geodetic Datum 1966'],
1976 'AGD84' :[-134, -48, 149,'AN','Australian Geodetic Datum 1984'],
1977 'ASTRO' :[-104,-129, 239,'IN','Camp Area Astro (Antarctica)' ],
1978 'CAPE' :[-136,-108,-292,'CD','CAPE (South Africa)' ],
1979 'ED50' :[ -87, -98,-121,'IN','European 1950' ],
1980 'ED79' :[ -86, -98,-119,'IN','European 1979' ],
1981 'GRB36' :[ 375,-111, 431,'AW','Great Britain 1936' ],
1982 'HAWAII':[ 89,-279,-183,'IN','Hawaiian Hawaii (Old)' ],
1983 'KAUAI' :[ 45,-290,-172,'IN','Hawaiian Kauai (Old)' ],
1984 'MAUI' :[ 65,-290,-190,'IN','Hawaiian Maui (Old)' ],
1985 'OAHU' :[ 56,-284,-181,'IN','Hawaiian Oahu (Old)' ],
1986 'INDIA' :[ 289, 734, 257,'EA','Indian' ],
1987 'NAD83' :[ 0, 0, 0,'RF','N. American 1983' ],
1988 'CANADA':[ -10, 158, 187,'CC','N. American Canada 1927' ],
1989 'ALASKA':[ -5, 135, 172,'CC','N. American Alaska 1927' ],
1990 'NAD27' :[ -8, 160, 176,'CC','N. American Conus 1927' ],
1991 'CARIBB':[ -7, 152, 178,'CC','N. American Caribbean' ],
1992 'MEXICO':[ -12, 130, 190,'CC','N. American Mexico' ],
1993 'CAMER' :[ 0, 125, 194,'CC','N. American Central America' ],
1994 'SAM56' :[-288, 175,-376,'IN','South American (Provisional1956)'],
1995 'SAM69' :[ -57, 1 , -41,'SA','South American 1969' ],
1996 'CAMPO' :[-148, 136, 90,'IN','S. American Campo Inchauspe (Argentina)'],
1997 'WGS72' :[ 0, 0 , 4.5,'WD','World Geodetic System - 72' ],
1998 'WGS84' :[ 0, 0 , 0,'WE','World Geodetic System - 84' ]}
2000 if datumcode not in datums:
2001 self.msg("unknown datum %s" % datumcode,priority="error")
2002 return -1
2004 datum=datums[datumcode]
2005 ellipsoid=datum[3]
2007 if ellipsoid not in ellipsoids:
2008 self.msg("unknown ellipsoid %s" % ellipsoid,priority="error")
2009 return -1
2011 if verbose:
2012 self.msg("Using %s datum with %s ellipsoid" % (datum[4],ellipsoids[ellipsoid][2]),origin="getdatum")
2013 return datum[0],datum[1],datum[2],ellipsoids[ellipsoid][0],ellipsoids[ellipsoid][1]
2019 def utm2long(self,east,north,zone,datum,nors):
2020 """
2021 this program converts universal transverse meractor coordinates to gps
2022 longitude and latitude (in radians)
2024 converted from Fortran by R. Indebetouw Jan 2009.
2025 orig. source: https://www.ngs.noaa.gov/TOOLS/program_descriptions.html
2026 ri also added other datums and ellipsoids in a helper function
2028 header from original UTMS fortran program:
2029 * this program was originally written by e. carlson
2030 * subroutines tmgrid, tconst, tmgeod, tconpc,
2031 * were written by t. vincenty, ngs, in july 1984 .
2032 * the orginal program was written in september of 1988.
2033 *
2034 * this program was updated on febuary 16, 1989. the update was
2035 * having the option of saving and *81* record file.
2036 *
2037 *
2038 * this program was updated on april 3, 1990. the following update
2039 * were made:
2040 * 1. change from just a choice of nad27 of nad83 reference
2041 * ellipsoids to; clarke 1866, grs80/wgs84, international, and
2042 * allow for user defined other.
2043 * 2. allow use of latitudes in southern hemisphere and longitudes
2044 * up to 360 degrees west.
2045 *
2046 * this program was updated on december 1, 1993. the following update
2047 * was made:
2048 * 1. the northings would compute the right latitude in the southern
2049 * hemishpere.
2050 * 2. the computed latitude on longidutes would be either in e or w.
2051 *
2052 ***************************************************************** *
2053 * disclaimer *
2054 * *
2055 * this program and supporting information is furnished by the *
2056 * government of the united states of america, and is accepted and *
2057 * used by the recipient with the understanding that the united states *
2058 * government makes no warranties, express or implied, concerning the *
2059 * accuracy, completeness, reliability, or suitability of this *
2060 * program, of its constituent parts, or of any supporting data. *
2061 * *
2062 * the government of the united states of america shall be under no *
2063 * liability whatsoever resulting from any use of this program. this *
2064 * program should not be relied upon as the sole basis for solving a *
2065 * problem whose incorrect solution could result in injury to person *
2066 * or property. *
2067 * *
2068 * this program is property of the government of the united states *
2069 * of america. therefore, the recipient further agrees not to assert *
2070 * proprietary rights therein and not to represent this program to *
2071 * anyone as being other than a government program. *
2072 * *
2073 ***********************************************************************
2075 this is the driver program to compute latitudes and longitudes from
2076 the utms for each zone
2078 input:
2079 northing, easting
2080 zone, datum
2081 nors=N/S
2083 determined according to datum:
2084 er = equatorial radius of the ellipsoid (semi-major axis) [m]
2085 rf = reciprocal of flatting of the ellipsod [unitless]
2086 f =
2087 esq= e squared
2089 calculated according to longitude and zone:
2090 rad = radian conversion factor
2091 cm = central meridian ( computed using the longitude)
2092 sf = scale factor of central meridian ( always .9996 for utm)
2093 orlim = southernmost parallel of latitude ( always zero for utm)
2094 r, a, b, c, u, v, w = ellipsoid constants used for computing
2095 meridional distance from latitude
2096 so = meridional distance (multiplied by scale factor )
2097 from the equator to the southernmost parallel of latitude
2098 ( always zero for utm)
2100 """
2102 rad=180./pl.pi
2104 offx,offy,offz,er,rf = self.getdatum(datum,verbose=self.verbose)
2106 f=1./rf
2107 esq=(2*f-f*f)
2109 # find the central meridian if the zone number is less than 30
2111 if zone < 30 : # ie W long - this code uses W=positive lon
2112 iz=zone
2113 icm=(183-(6*iz))
2114 cm=float(icm)/rad
2115 ucm=(icm+3)/rad
2116 lcm=(icm-3)/rad
2117 else:
2118 iz=zone
2119 icm=(543 - (6*iz))
2120 cm= float(icm)/rad
2121 ucm=(icm+3)/rad
2122 lcm=(icm-3)/rad
2124 tol=(5./60.)/rad
2126 if nors == 'S':
2127 fn= 10000000.
2128 else:
2129 fn=0.
2131 fe=500000.0
2132 sf=0.9996
2133 orlim=0.0
2135 found=0
2137 # precompute parameters for this zone:
2138 eps,r,so,v0,v2,v4,v6 = self.tconpc(sf,orlim,er,esq,rf)
2140 # compute the latitudes and longitudes:
2141 lat,lon,conv,kp = self.tmgeod(north,east,eps,cm,fe,sf,so,r,v0,v2,v4,v6,fn,er,esq)
2143 # do the test to see if the longitude is within 5 minutes
2144 # of the boundaries for the zone and if so compute the
2145 # north and easting for the adjacent zone
2147 # if abs(ucm-lam) < tol:
2148 # cm=float(icm+6)/rad
2149 # iz=iz-1
2150 # if iz == 0:
2151 # iz=60
2152 # found=found+1
2153 # lat,lon,conv,kp = tmgeod(n,e,eps,cm,fe,sf,so,r,v0,v2,v4,v6,fn,er,esq)
2154 #
2155 # if abs(lcm-lam) < tol:
2156 # cm=float(icm-6)/rad
2157 # iz=iz+1
2158 # if iz == 61:
2159 # iz=1
2160 # lat,lon,conv,kp = tmgeod(n,e,eps,cm,fe,sf,so,r,v0,v2,v4,v6,fn,er,esq)
2161 # found=found+1
2163 # *** convert to more usual convention of negative lon = W
2164 lon=-lon
2166 if self.verbose:
2167 self.msg(" longitude, latitude = %s %s; conv,kp = %f,%f" % (qa.angle(qa.quantity(lon,"rad"),prec=8)[0],qa.angle(qa.quantity(lat,"rad"),prec=8)[0],conv,kp),origin="utm2long")
2169 return lon,lat
2174 def long2xyz(self,lon,lat,elevation,datum):
2175 """
2176 Returns the nominal ITRF (X, Y, Z) coordinates [m] for a point at
2177 geodetic latitude and longitude [radians] and elevation [m].
2178 The ITRF frame used is not the official ITRF, just a right
2179 handed Cartesian system with X going through 0 latitude and 0 longitude,
2180 and Z going through the north pole.
2181 orig. source: http://www.oc.nps.edu/oc2902w/coord/llhxyz.htm
2182 """
2184 # geodesy/NGS/XYZWIN/
2185 # http://www.oc.nps.edu/oc2902w/coord/llhxyz.htm
2186 dx,dy,dz,er,rf = self.getdatum(datum,verbose=False)
2188 f=1./rf
2189 esq=2*f-f**2
2190 nu=er/pl.sqrt(1.-esq*(pl.sin(lat))**2)
2192 x=(nu+elevation)*pl.cos(lat)*pl.cos(lon) +dx
2193 y=(nu+elevation)*pl.cos(lat)*pl.sin(lon) +dy
2194 z=((1.-esq)*nu+elevation)*pl.sin(lat) +dz
2196 # https://www.ngs.noaa.gov/cgi-bin/xyz_getxyz.prl
2197 # S231200.0 W0674800.0 0.0000
2198 # > 2216194.4188 -5430618.6472 -2497092.5213 GRS80
2200 return x,y,z
2203 def xyz2long(self,x,y,z,datum):
2204 """
2205 Given ITRF Earth-centered (X, Y, Z) coordinates [m] for a point,
2206 returns geodetic latitude and longitude [radians] and elevation [m]
2207 The ITRF frame used is not the official ITRF, just a right
2208 handed Cartesian system with X going through 0 latitude and 0 longitude,
2209 and Z going through the north pole.
2210 Elevation is measured relative to the closest point to
2211 the (latitude, longitude) on the WGS84 reference
2212 ellipsoid.
2213 orig. source: http://www.iausofa.org/2013_1202_C/sofa/gc2gde.html
2214 """
2215 # http://www.iausofa.org/
2216 # http://www.iausofa.org/2013_1202_C/sofa/gc2gde.html
2218 dx,dy,dz,er,rf = self.getdatum(datum,verbose=False)
2219 f=1./rf
2221 # Validate ellipsoid parameters.
2222 if ( f < 0.0 or f >= 1.0 ): return -1,-1,-1
2223 if ( er <= 0.0 ): return -1,-1,-1
2225 #Functions of ellipsoid parameters (with further validation of f).
2226 e2 = (2.0 - f) * f
2227 e4t = e2*e2 * 1.5
2228 ec2 = 1.0 - e2 # = 1-esq = (1-f)**2
2229 if ( ec2 <= 0.0 ): return -1,-1,-1
2230 ec = pl.sqrt(ec2) # =1-f
2231 b = er * ec
2233 # Distance from polar axis
2234 r = pl.sqrt( (x-dx)**2 + (y-dy)**2 )
2236 # Longitude.
2237 if r>0.:
2238 lon = pl.arctan2(y-dx, x-dx)
2239 else:
2240 lon = 0.
2242 # Prepare Newton correction factors.
2243 s0 = abs(z-dz) / er
2244 c0 = ec * r / er
2246 a0 = pl.sqrt( c0**2 + s0**2 )
2247 d0 = ec* s0* a0**3 + e2* s0**3
2248 f0 = r/er* a0**3 - e2* c0**3
2250 # Prepare Halley correction factor.
2251 b0 = e4t * s0**2 * c0**2 * r/er * (a0 - ec)
2252 s1 = d0*f0 - b0*s0
2253 cc = ec * (f0*f0 - b0*c0)
2255 # Evaluate latitude and height. */
2256 lat = pl.arctan2(s1,cc)
2257 height = (r*cc + abs(z-dz)*s1 - er*pl.sqrt(ec2*s1**2 + cc**2))/pl.sqrt(s1**2 + cc**2)
2259 # Restore sign of latitude.
2260 if (z-dz) < 0: lat = -lat
2262 return lon,lat,height
2265 def xyz2long_old(self,x,y,z,datum):
2267 dx,dy,dz,er,rf = self.getdatum(datum,verbose=False)
2269 f=1./rf
2271 b= ((x-dx)**2 + (y-dy)**2) / er**2
2272 c= (z-dx)**2 / er**2
2274 esq=2*f-f**2 # (a2-b2)/a2
2276 a0=c*(1-esq)
2277 a1=2*a0
2278 efth=esq**2
2279 a2=a0+b-efth
2280 a3=-2.*efth
2281 a4=-efth
2283 b0=4.*a0
2284 b1=3.*a1
2285 b2=2.*a2
2286 b3=a3
2288 # refine/calculate esq
2289 nlqk=esq
2290 for i in range(3):
2291 nlqks = nlqk * nlqk
2292 nlqkc = nlqk * nlqks
2293 nlf = a0*nlqks*nlqks + a1*nlqkc + a2*nlqks + a3*nlqk + a4
2294 nlfprm = b0*nlqkc + b1*nlqks + b2*nlqk + b3
2295 nlqk = nlqk - (nlf / nlfprm)
2297 y0 = (1.+nlqk)*(z-dz)
2298 x0 = pl.sqrt((x-dx)**2 + (y-dy)**2)
2299 lat=pl.arctan2(y0,x0)
2300 lon=pl.arctan2(y-dy,x-dx)
2301 #casalog.post... x-dx,y-dy,z-dz,x0,y0
2303 return lon,lat
2306 def utm2xyz(self,easting,northing,elevation,zone,datum,nors):
2307 """
2308 Returns the nominal ITRF (X, Y, Z) coordinates [m] for a point at
2309 UTM easting, northing, elevation [m], and zone of a given
2310 datum (e.g. 'WGS84') and north/south flag ("N" or "S").
2311 The ITRF frame used is not the official ITRF, just a right
2312 handed Cartesian system with X going through 0 latitude and 0 longitude,
2313 and Z going through the north pole.
2314 """
2316 lon,lat = self.utm2long(easting,northing,zone,datum,nors)
2317 x,y,z = self.long2xyz(lon,lat,elevation,datum)
2319 return x,y,z
2322 def locxyz2itrf(self, lat, longitude, alt, locx=0.0, locy=0.0, locz=0.0):
2323 """
2324 Returns the nominal ITRF (X, Y, Z) coordinates [m] for a point at
2325 "local" (x, y, z) [m] measured at geodetic latitude lat and longitude
2326 longitude (degrees) and altitude of the reference point of alt.
2327 The ITRF frame used is not the official ITRF, just a right
2328 handed Cartesian system with X going through 0 latitude and 0 longitude,
2329 and Z going through the north pole. The "local" (x, y, z) are measured
2330 relative to the closest point to (lat, longitude) on the WGS84 reference
2331 ellipsoid, with z normal to the ellipsoid and y pointing north.
2332 """
2333 # from Rob Reid; need to generalize to use any datum...
2334 phi, lmbda = map(pl.radians, (lat, longitude))
2335 sphi = pl.sin(phi)
2336 a = 6378137.0 # WGS84 equatorial semimajor axis
2337 b = 6356752.3142 # WGS84 polar semimajor axis
2338 ae = pl.arccos(b / a)
2339 N = a / pl.sqrt(1.0 - (pl.sin(ae) * sphi)**2)
2341 # Now you see the connection between the Old Ones and Antarctica...
2342 Nploczcphimlocysphi = (N + locz+alt) * pl.cos(phi) - locy * sphi
2344 clmb = pl.cos(lmbda)
2345 slmb = pl.sin(lmbda)
2347 x = Nploczcphimlocysphi * clmb - locx * slmb
2348 y = Nploczcphimlocysphi * slmb + locx * clmb
2349 z = (N * (b / a)**2 + locz+alt) * sphi + locy * pl.cos(phi)
2351 return x, y, z
2356 def itrf2loc(self, x,y,z, cx,cy,cz):
2357 """
2358 Given Earth-centered ITRF (X, Y, Z) coordinates [m]
2359 and the Earth-centered coords of the center of array,
2360 Returns local (x, y, z) [m] relative to the center of the array,
2361 oriented with x and y tangent to the closest point
2362 at the COFA (latitude, longitude) on the WGS84 reference
2363 ellipsoid, with z normal to the ellipsoid and y pointing north.
2364 """
2365 clon,clat,h = self.xyz2long(cx,cy,cz,'WGS84')
2366 ccoslon=pl.cos(clon)
2367 csinlon=pl.sin(clon)
2368 csinlat=pl.sin(clat)
2369 ccoslat=pl.cos(clat)
2371 if isinstance(x,float): # weak
2372 x=[x]
2373 y=[y]
2374 z=[z]
2375 n=x.__len__()
2376 lat=pl.zeros(n)
2377 lon=pl.zeros(n)
2378 el=pl.zeros(n)
2380 # do like MsPlotConvert
2381 for i in range(n):
2382 # translate w/o rotating:
2383 xtrans=x[i]-cx
2384 ytrans=y[i]-cy
2385 ztrans=z[i]-cz
2386 # rotate
2387 lat[i] = (-csinlon*xtrans) + (ccoslon*ytrans)
2388 lon[i] = (-csinlat*ccoslon*xtrans) - (csinlat*csinlon*ytrans) + ccoslat*ztrans
2389 el[i] = (ccoslat*ccoslon*xtrans) + (ccoslat*csinlon*ytrans) + csinlat*ztrans
2391 return lat,lon,el
2397 def itrf2locname(self, x,y,z, obsname):
2398 """
2399 Given Earth-centered ITRF (X, Y, Z) coordinates [m]
2400 and the name of an known array (see me.obslist()),
2401 Returns local (x, y, z) [m] relative to the center of the array,
2402 oriented with x and y tangent to the closest point
2403 at the COFA (latitude, longitude) on the WGS84 reference
2404 ellipsoid, with z normal to the ellipsoid and y pointing north.
2405 """
2406 cofa=me.measure(me.observatory(obsname),'WGS84')
2407 cx,cy,cz=self.long2xyz(cofa['m0']['value'],cofa['m1']['value'],cofa['m2']['value'],cofa['refer'])
2409 return self.itrf2loc(x,y,z,cx,cy,cz)
2423 ###########################################################
2425 def plotants(self,x,y,z,d,name):
2426 # given globals
2428 cx=pl.mean(x)
2429 cy=pl.mean(y)
2430 cz=pl.mean(z)
2431 lat,lon,el = self.itrf2loc(x,y,z,cx,cy,cz)
2432 n=lat.__len__()
2434 dolam=0
2435 # TODO convert to klam: (d too)
2436 ###
2438 rg=max(lat)-min(lat)
2439 r2=max(lon)-min(lon)
2440 if r2>rg:
2441 rg=r2
2442 if max(d)>0.01*rg:
2443 pl.plot(lat,lon,',')
2444 #casalog.post(max(d),ra)
2445 for i in range(n):
2446 pl.gca().add_patch(pl.Circle((lat[i],lon[i]),radius=0.5*d[i],fc="#dddd66"))
2447 if n<10:
2448 pl.text(lat[i],lon[i],name[i],horizontalalignment='center',verticalalignment='center')
2449 else:
2450 pl.plot(lat,lon,'o',c="#dddd66")
2451 if n<10:
2452 for i in range(n):
2453 pl.text(lat[i],lon[i],name[i],horizontalalignment='center',fontsize=8)
2455 pl.axis('equal')
2456 #if dolam:
2457 # pl.xlabel("kilolamda")
2458 # pl.ylabel("kilolamda")
2478 ######################################################
2479 # helper function to get the pixel size from an image (positive increments)
2481 def cellsize(self,image):
2482 ia.open(image)
2483 mycs=ia.coordsys()
2484 ia.done()
2485 increments=mycs.increment(type="direction")['numeric']
2486 cellx=qa.quantity(abs(increments[0]),mycs.units(type="direction")[0])
2487 celly=qa.quantity(abs(increments[1]),mycs.units(type="direction")[1])
2488 xform=mycs.lineartransform(type="direction")
2489 offdiag=max(abs(xform[0,1]),abs(xform[1,0]))
2490 if offdiag > 1e-4:
2491 self.msg("Your image is rotated with respect to Lat/Lon. I can't cope with that yet",priority="error")
2492 cellx=qa.mul(cellx,abs(xform[0,0]))
2493 celly=qa.mul(celly,abs(xform[1,1]))
2494 return [cellx,celly]
2496 ######################################################
2497 # helper function to get the spectral size from an image
2499 def spectral(self,image):
2500 ia.open(image)
2501 cs=ia.coordsys()
2502 sh=ia.shape()
2503 ia.done()
2504 spc = cs.findcoordinate("spectral")
2505 if not spc['return']:
2506 return 0.0
2507 model_width=str(cs.increment(type="spectral")['numeric'][0])+cs.units(type="spectral")[0]
2508 model_nchan=sh[spc['pixel']]
2509 return model_nchan,model_width
2511 ######################################################
2513 def is4d(self, image):
2514 ia.open(image)
2515 s=ia.shape()
2516 if len(s)!=4: return False
2517 cs=ia.coordsys()
2518 ia.done()
2519 dir=cs.findcoordinate("direction")
2520 spc=cs.findcoordinate("spectral")
2521 stk=cs.findcoordinate("stokes")
2522 if not (dir['return'] and spc['return'] and stk['return']): return False
2523 if dir['pixel'].__len__() != 2: return False
2524 if spc['pixel'].__len__() != 1: return False
2525 if stk['pixel'].__len__() != 1: return False
2526 # they have to be in the correct order too
2527 if stk['pixel']!=2: return False
2528 if spc['pixel']!=3: return False
2529 if dir['pixel'][0]!=0: return False
2530 if dir['pixel'][1]!=1: return False
2531 cs.done()
2532 return True
2534 ##################################################################
2535 # fit modelimage into a 4 coordinate image defined by the parameters
2537 # TODO spectral extrapolation and regridding using innchan ****
2539 def modifymodel(self, inimage, outimage,
2540 inbright,indirection,incell,incenter,inwidth,innchan,
2541 flatimage=False): # if nonzero, create mom -1 image
2543 # new ia tool
2544 in_ia=ia.newimagefromfile(inimage)
2545 in_shape=in_ia.shape()
2546 in_csys=in_ia.coordsys()
2548 # pull data first, since ia.stats doesn't work w/o a CS:
2549 if outimage!=inimage:
2550 if self.verbose: self.msg("rearranging input data (may take some time for large cubes)",origin="setup model")
2551 arr=in_ia.getchunk()
2552 else:
2553 # TODO move rearrange to inside ia tool, and at least don't do this:
2554 arr=pl.zeros(in_shape)
2555 axmap=[-1,-1,-1,-1]
2556 axassigned=[-1,-1,-1,-1]
2560 # brightness scaling
2561 if (inbright=="unchanged") or (inbright==""):
2562 scalefactor=1.
2563 else:
2564 if self.isquantity(inbright,halt=False):
2565 qinb=qa.quantity(inbright)
2566 if qinb['unit']!='':
2567 # qa doesn't deal universally well with pixels and beams
2568 # so this may fail:
2569 try:
2570 inb=qa.convert(qinb,"Jy/pixel")['value']
2571 except:
2572 inb=qinb['value']
2573 self.msg("assuming inbright="+str(inbright)+" means "+str(inb)+" Jy/pixel",priority="warn")
2574 inbright=inb
2575 try:
2576 scalefactor=float(inbright)/pl.nanmax(arr)
2577 except Exception:
2578 in_ia.close()
2579 raise
2581 # check shape characteristics of the input;
2582 # add degenerate axes as neeed:
2584 in_dir=in_csys.findcoordinate("direction")
2585 in_spc=in_csys.findcoordinate("spectral")
2586 in_stk=in_csys.findcoordinate("stokes")
2589 # first check number of pixel axes and raise to 4 if required
2590 in_nax=arr.shape.__len__()
2591 if in_nax<2:
2592 in_ia.close()
2593 self.msg("Your input model has fewer than 2 dimensions. Can't proceed",priority="error")
2594 return False
2595 if in_nax==2:
2596 arr=arr.reshape([arr.shape[0],arr.shape[1],1])
2597 in_shape=arr.shape
2598 in_nax=in_shape.__len__() # which should be 3
2599 if in_nax==3:
2600 arr=arr.reshape([arr.shape[0],arr.shape[1],arr.shape[2],1])
2601 in_shape=arr.shape
2602 in_nax=in_shape.__len__() # which should be 4
2603 if in_nax>4:
2604 in_ia.close()
2605 self.msg("model image has more than 4 dimensions. Not sure how to proceed",priority="error")
2606 return False
2609 # make incell a list if not already
2610 try:
2611 if type(incell) == type([]):
2612 incell = map(qa.convert,incell,['arcsec','arcsec'])
2613 else:
2614 incell = qa.abs(qa.convert(incell,'arcsec'))
2615 # incell[0]<0 for RA increasing left
2616 incell = [qa.mul(incell,-1),incell]
2617 except Exception:
2618 # invalid incell
2619 in_ia.close()
2620 raise
2621 # later, we can test validity with qa.compare()
2624 # now parse coordsys:
2625 model_refdir=""
2626 model_cell=""
2627 # look for direction coordinate, with two pixel axes:
2628 if in_dir['return']:
2629 in_ndir = in_dir['pixel'].__len__()
2630 if in_ndir != 2:
2631 self.msg("Mal-formed direction coordinates in modelimage. Discarding and using first two pixel axes for RA and Dec.")
2632 axmap[0]=0 # direction in first two pixel axes
2633 axmap[1]=1
2634 axassigned[0]=0 # coordinate corresponding to first 2 pixel axes
2635 axassigned[1]=0
2636 else:
2637 # we've found direction axes, and may change their increments and direction or not.
2638 dirax=in_dir['pixel']
2639 axmap[0]=dirax[0]
2640 axmap[1]=dirax[1]
2641 axassigned[dirax[0]]=0
2642 axassigned[dirax[1]]=0
2643 if self.verbose: self.msg("Direction coordinate (%i,%i) parsed" % (axmap[0],axmap[1]),origin="setup model")
2644 # move model_refdir to center of image
2645 model_refpix=[0.5*in_shape[axmap[0]],0.5*in_shape[axmap[1]]]
2646 ra = in_ia.toworld(model_refpix,'q')['quantity']['*'+str(axmap[0]+1)]
2647 dec = in_ia.toworld(model_refpix,'q')['quantity']['*'+str(axmap[1]+1)]
2648 model_refdir= in_csys.referencecode(type="direction",list=False)[0]+" "+qa.formxxx(ra,format='hms',prec=5)+" "+qa.formxxx(dec,format='dms',prec=5)
2649 model_projection=in_csys.projection()['type']
2650 model_projpars=in_csys.projection()['parameters']
2652 # cell size from image
2653 increments=in_csys.increment(type="direction")['numeric']
2654 incellx=qa.quantity(increments[0],in_csys.units(type="direction")[0])
2655 incelly=qa.quantity(increments[1],in_csys.units(type="direction")[1])
2656 xform=in_csys.lineartransform(type="direction")
2657 offdiag=max(abs(xform[0,1]),abs(xform[1,0]))
2658 if offdiag > 1e-4:
2659 self.msg("Your image is rotated with respect to Lat/Lon. I can't cope with that yet",priority="error")
2660 incellx=qa.mul(incellx,xform[0,0])
2661 incelly=qa.mul(incelly,xform[1,1])
2663 # preserve sign in case input image is backwards (RA increasing right)
2664 model_cell = [qa.convert(incellx,'arcsec'),qa.convert(incelly,'arcsec')]
2666 else: # no valid direction coordinate
2667 axmap[0]=0 # assign direction to first two pixel axes
2668 axmap[1]=1
2669 axassigned[0]=0 # assign coordinate corresponding to first 2 pixel axes
2670 axassigned[1]=0
2673 # try to parse direction using splitter function and override model_refdir
2674 if type(indirection)==type([]):
2675 if len(indirection) > 1:
2676 in_ia.close()
2677 self.msg("error parsing indirection "+str(indirection)+" -- should be a single direction string")
2678 return False
2679 else:
2680 indirection=indirection[0]
2681 if self.isdirection(indirection,halt=False):
2682 # lon/lat = ra/dec if J2000, =glon/glat if galactic
2683 epoch, lon, lat = self.direction_splitter(indirection)
2685 model_refdir=epoch+qa.formxxx(lon,format='hms',prec=5)+" "+qa.formxxx(lat,format='dms',prec=5)
2686 model_refpix=[0.5*in_shape[axmap[0]],0.5*in_shape[axmap[1]]]
2687 model_projection="SIN" # for indirection we default to SIN.
2688 model_projpars=pl.array([0.,0])
2689 if self.verbose: self.msg("setting model image direction to indirection = "+model_refdir,origin="setup model")
2690 else:
2691 # indirection is not set - is there a direction in the model already?
2692 if not self.isdirection(model_refdir,halt=False):
2693 in_ia.close()
2694 self.msg("Cannot determine reference direction in model image. Valid 'indirection' parameter must be provided.",priority="error")
2695 return False
2698 # override model_cell?
2699 cell_replaced=False
2700 if self.isquantity(incell[0],halt=False):
2701 if qa.compare(incell[0],"1arcsec"):
2702 model_cell=incell
2703 cell_replaced=True
2704 if self.verbose: self.msg("replacing existing model cell size with incell",origin="setup model")
2705 valid_modcell=False
2706 if not cell_replaced:
2707 if self.isquantity(model_cell[0],halt=False):
2708 if qa.compare(model_cell[0],"1arcsec"):
2709 valid_modcell=True
2710 if not valid_modcell:
2711 in_ia.close()
2712 self.msg("Unable to determine model cell size. Valid 'incell' parameter must be provided.",priority="error")
2713 return False
2717 if self.verbose:
2718 self.msg("model image shape="+str(in_shape),origin="setup model")
2719 self.msg("model pixel = %8.2e x %8.2e arcsec" % (model_cell[0]['value'],model_cell[1]['value']),origin="setup model")
2727 # we've now found or assigned two direction axes, and changed direction and cell if required
2728 # next, work on spectral axis:
2730 model_specrefval=""
2731 model_width=""
2732 # look for a spectral axis:
2733 if in_spc['return']:
2734 #if type(in_spc[1]) == type(1) :
2735 # # should not come here after SWIG switch over
2736 # foo=in_spc[1]
2737 #else:
2738 foo=in_spc['pixel'][0]
2739 if in_spc['pixel'].__len__() > 1:
2740 self.msg("you seem to have two spectral axes",priority="warn")
2741 model_nchan=arr.shape[foo]
2742 axmap[3]=foo
2743 axassigned[foo]=3
2744 model_restfreq=in_csys.restfrequency()
2745 model_specrefpix=in_csys.referencepixel(type="spectral")['numeric'][0]
2746 model_width =in_csys.increment(type="spectral")['numeric'][0]
2747 model_specrefval=in_csys.referencevalue(type="spectral")['numeric'][0]
2748 # make quantities
2749 model_width=qa.quantity(model_width,in_csys.units(type="spectral")[0])
2750 model_specrefval=qa.quantity(model_specrefval,in_csys.units(type="spectral")[0])
2751 add_spectral_coord=False
2752 if self.verbose: self.msg("Spectral Coordinate %i parsed" % axmap[3],origin="setup model")
2753 else:
2754 # need to add one to the coordsys
2755 add_spectral_coord=True
2758 if add_spectral_coord:
2759 # find first unused axis - probably at end, but just in case its not:
2760 i=0
2761 extra_axis=-1
2762 while extra_axis<0 and i<4:
2763 if axassigned[i]<0: extra_axis=i
2764 i+=1
2765 if extra_axis<0:
2766 in_ia.close()
2767 self.msg("I can't find an unused axis to make Spectral [%i %i %i %i] " % (axassigned[0],axassigned[1],axassigned[2],axassigned[3]),priority="error",origin="setup model")
2768 return False
2770 axmap[3]=extra_axis
2771 axassigned[extra_axis]=3
2772 model_nchan=arr.shape[extra_axis]
2776 # override specrefval?
2777 specref_replaced=False
2778 if self.isquantity(incenter,halt=False):
2779 if qa.compare(incenter,"1Hz"):
2780 if (qa.quantity(incenter))['value']>=0:
2781 model_specrefval=incenter
2782 model_specrefpix=pl.floor(model_nchan*0.5)
2783 model_restfreq=incenter
2784 specref_replaced=True
2785 if self.verbose: self.msg("setting central frequency to "+incenter,origin="setup model")
2786 valid_modspec=False
2787 if not specref_replaced:
2788 if self.isquantity(model_specrefval,halt=False):
2789 if qa.compare(model_specrefval,"1Hz"):
2790 valid_modspec=True
2791 if not valid_modspec:
2792 in_ia.close()
2793 self.msg("Unable to determine model frequency. Valid 'incenter' parameter must be provided.",priority="error")
2794 return False
2796 # override inwidth?
2797 width_replaced=False
2798 if self.isquantity(inwidth,halt=False):
2799 if qa.compare(inwidth,"1Hz"):
2800 if (qa.quantity(inwidth))['value']>=0:
2801 model_width=inwidth
2802 width_replaced=True
2803 if self.verbose: self.msg("setting channel width to "+inwidth,origin="setup model")
2804 valid_modwidth=False
2805 if not width_replaced:
2806 if self.isquantity(model_width,halt=False):
2807 if qa.compare(model_width,"1Hz"):
2808 valid_modwidth=True
2809 if not valid_modwidth:
2810 in_ia.close()
2811 self.msg("Unable to determine model channel or bandwidth. Valid 'inwidth' parameter must be provided.",priority="error")
2812 return False
2818 model_stokes=""
2819 # look for a stokes axis
2820 if in_stk['return']:
2821 model_stokes=in_csys.stokes()
2822 foo=model_stokes[0]
2823 out_nstk=model_stokes.__len__()
2824 for i in range(out_nstk-1):
2825 foo=foo+model_stokes[i+1]
2826 model_stokes=foo
2827 #if type(in_stk[1]) == type(1):
2828 # # should not come here after SWIG switch over
2829 # foo=in_stk[1]
2830 #else:
2831 foo=in_stk['pixel'][0]
2832 if in_stk['pixel'].__len__() > 1:
2833 self.msg("you seem to have two stokes axes",priority="warn")
2834 axmap[2]=foo
2835 axassigned[foo]=2
2836 if in_shape[foo]>4:
2837 in_ia.close()
2838 self.msg("you appear to have more than 4 Stokes components - please edit your header and/or parameters",priority="error")
2839 return False
2840 add_stokes_coord=False
2841 if self.verbose: self.msg("Stokes Coordinate %i parsed" % axmap[2],origin="setup model")
2842 else:
2843 # need to add one to the coordsys
2844 add_stokes_coord=True
2848 if add_stokes_coord:
2849 # find first unused axis - probably at end, but just in case its not:
2850 i=0
2851 extra_axis=-1
2852 while extra_axis<0 and i<4:
2853 if axassigned[i]<0: extra_axis=i
2854 i+=1
2855 if extra_axis<0:
2856 in_ia.close()
2857 self.msg("I can't find an unused axis to make Stokes [%i %i %i %i] " % (axassigned[0],axassigned[1],axassigned[2],axassigned[3]),priority="error",origin="setup model")
2858 return False
2859 axmap[2]=extra_axis
2860 axassigned[extra_axis]=2
2862 if arr.shape[extra_axis]>4:
2863 in_ia.close()
2864 self.msg("you have %i Stokes parameters in your potential Stokes axis %i. something is wrong." % (arr.shape[extra_axis],extra_axis),priority="error")
2865 return False
2866 if self.verbose: self.msg("Adding Stokes Coordinate",origin="setup model")
2867 if arr.shape[extra_axis]==4:
2868 model_stokes="IQUV"
2869 if arr.shape[extra_axis]==3:
2870 model_stokes="IQV"
2871 self.msg("setting IQV Stokes parameters from the 4th axis of you model. If that's not what you want, then edit the header",origin="setup model",priority="warn")
2872 if arr.shape[extra_axis]==2:
2873 model_stokes="IQ"
2874 self.msg("setting IQ Stokes parameters from the 4th axis of you model. If that's not what you want, then edit the header",origin="setup model",priority="warn")
2875 if arr.shape[extra_axis]<=1:
2876 model_stokes="I"
2881 # ========================================
2882 # now we should have 4 assigned pixel axes, and model_cell, model_refdir, model_nchan,
2883 # model_stokes all set
2884 out_nstk=len(model_stokes)
2885 if self.verbose:
2886 self.msg("axis map for model image = %i %i %i %i" %
2887 (axmap[0],axmap[1],axmap[2],axmap[3]),origin="setup model")
2889 modelshape=[in_shape[axmap[0]], in_shape[axmap[1]],out_nstk,model_nchan]
2890 if outimage!=inimage:
2891 ia.fromshape(outimage,modelshape,overwrite=True)
2892 modelcsys=ia.coordsys()
2893 else:
2894 modelcsys=in_ia.coordsys()
2895 modelcsys.setunits(['rad','rad','','Hz'])
2896 modelcsys.setincrement([qa.convert(model_cell[0],modelcsys.units()[0])['value'], # already -1
2897 qa.convert(model_cell[1],modelcsys.units()[1])['value']],
2898 type="direction")
2900 dirm=self.dir_s2m(model_refdir)
2901 lonq=dirm['m0']
2902 latq=dirm['m1']
2903 modelcsys.setreferencecode(dirm['refer'],type="direction")
2904 if len(model_projpars)>0:
2905 modelcsys.setprojection(parameters=model_projpars.tolist(),type=model_projection)
2906 else:
2907 modelcsys.setprojection(type=model_projection)
2908 modelcsys.setreferencevalue(
2909 [qa.convert(lonq,modelcsys.units()[0])['value'],
2910 qa.convert(latq,modelcsys.units()[1])['value']],
2911 type="direction")
2912 modelcsys.setreferencepixel(model_refpix,"direction")
2913 if self.verbose:
2914 self.msg("sky model image direction = "+model_refdir,origin="setup model")
2915 self.msg("sky model image increment = "+str(model_cell[0]),origin="setup model")
2917 modelcsys.setspectral(refcode="LSRK",restfreq=model_restfreq)
2918 modelcsys.setreferencevalue(qa.convert(model_specrefval,modelcsys.units()[3])['value'],type="spectral")
2919 modelcsys.setreferencepixel(model_specrefpix,type="spectral") # but not half-pixel
2920 modelcsys.setincrement(qa.convert(model_width,modelcsys.units()[3])['value'],type="spectral")
2923 # first assure that the csys has the expected order
2924 expected=['Direction', 'Direction', 'Stokes', 'Spectral']
2925 if modelcsys.axiscoordinatetypes() != expected:
2926 in_ia.close()
2927 self.msg("internal error with coordinate axis order created by Imager",priority="error")
2928 self.msg(modelcsys.axiscoordinatetypes().__str__(),priority="error")
2929 return False
2931 # more checks:
2932 foo=pl.array(modelshape)
2933 if not (pl.array(arr.shape) == pl.array(foo.take(axmap).tolist())).all():
2934 in_ia.close()
2935 self.msg("internal error: I'm confused about the shape if your model data cube",priority="error")
2936 self.msg("have "+foo.take(axmap).__str__()+", want "+in_shape.__str__(),priority="error")
2937 return False
2939 if outimage!=inimage:
2940 ia.setcoordsys(modelcsys.torecord())
2941 ia.done()
2942 ia.open(outimage)
2944 # now rearrange the pixel axes if ness.
2945 for ax in range(4):
2946 if axmap[ax] != ax:
2947 if self.verbose: self.msg("swapping input axes %i with %i" % (ax,axmap[ax]),origin="setup model")
2948 arr=arr.swapaxes(ax,axmap[ax])
2949 tmp=axmap[ax]
2950 axmap[ax]=ax
2951 axmap[tmp]=tmp
2954 # there's got to be a better way to remove NaNs:
2955 if outimage!=inimage:
2956 for i0 in range(arr.shape[0]):
2957 for i1 in range(arr.shape[1]):
2958 for i2 in range(arr.shape[2]):
2959 for i3 in range(arr.shape[3]):
2960 foo=arr[i0,i1,i2,i3]
2961 if foo!=foo: arr[i0,i1,i2,i3]=0.0
2963 if self.verbose and outimage!=inimage:
2964 self.msg("model array minmax= %e %e" % (arr.min(),arr.max()),origin="setup model")
2965 self.msg("scaling model brightness by a factor of %f" % scalefactor,origin="setup model")
2966 self.msg("image channel width = %8.2e GHz" % qa.convert(model_width,'GHz')['value'],origin="setup model")
2967 if arr.nbytes > 5e7:
2968 self.msg("your model is large - predicting visibilities may take a while.",priority="warn")
2970 if outimage!=inimage:
2971 ia.putchunk(arr*scalefactor)
2972 ia.setbrightnessunit("Jy/pixel")
2973 ia.close()
2974 in_ia.close()
2975 # outimage should now have correct Coordsys and shape
2977 # make a moment 0 image
2978 if flatimage and outimage!=inimage:
2979 self.flatimage(outimage,verbose=self.verbose)
2981 # returning to the outside world we'll return a positive cell:
2982 model_cell=[qa.abs(model_cell[0]),qa.abs(model_cell[1])]
2983 model_size=[qa.mul(modelshape[0],model_cell[0]),
2984 qa.mul(modelshape[1],model_cell[1])]
2986 if self.verbose: self.msg(" ") # add a line after my spewage
2988 return model_refdir,model_cell,model_size,model_nchan,model_specrefval,model_specrefpix,model_width,model_stokes
2994 ##################################################################
2995 # image/clean subtask
2997 def imclean(self,mstoimage,imagename,
2998 cleanmode,psfmode,cell,imsize,imcenter,
2999 interactive,niter,threshold,weighting,
3000 outertaper,pbcor,stokes,sourcefieldlist="",
3001 modelimage="",mask=[],dryrun=False):
3002 """
3003 Wrapper function to call CASA imaging task 'clean' on a MeasurementSet
3004 mstoimage parameter expects the path to a MeasurementSet
3005 imsize parameter expects a length-2 list of integers
3006 cell parameter expects a length-2 list containing qa.quantity objects
3007 interactive and dryrun parameters expect boolean type input
3008 Other parameters expect input in format compatible with 'clean'
3010 No return
3012 Creates a template '[imagename].clean.last' file in addition to
3013 outputs of task 'clean'
3014 """
3016 # determine channelization from (first) ms:
3017 if is_array_type(mstoimage):
3018 ms0=mstoimage[0]
3019 if len(mstoimage)==1:
3020 mstoimage=mstoimage[0]
3021 else:
3022 ms0=mstoimage
3024 if os.path.exists(ms0):
3025 tb.open(ms0+"/SPECTRAL_WINDOW")
3026 if tb.nrows() > 1:
3027 self.msg("determining output cube parameters from FIRST of several SPW in MS "+ms0)
3028 freq=tb.getvarcol("CHAN_FREQ")['r1']
3029 nchan=freq.size
3030 tb.done()
3031 elif dryrun:
3032 nchan=1 # May be wrong
3034 if nchan==1:
3035 chanmode="mfs"
3036 else:
3037 chanmode="channel"
3039 psfmode="clark"
3040 ftmachine="ft"
3041 imagermode="clark" # set default to prevent UnboundLocalError
3043 if cleanmode=="csclean":
3044 imagermode='csclean'
3045 #if cleanmode=="clark":
3046 # imagermode=""
3047 if cleanmode=="mosaic":
3048 imagermode="mosaic"
3049 ftmachine="mosaic"
3051 # in 3.4 clean doesn't accept just any imsize
3052 optsize=[0,0]
3053 optsize[0]=_su.getOptimumSize(imsize[0])
3054 nksize=len(imsize)
3055 if nksize==1: # imsize can be a single element or array
3056 optsize[1]=optsize[0]
3057 else:
3058 optsize[1]=_su.getOptimumSize(imsize[1])
3059 if((optsize[0] != imsize[0]) or (nksize!=1 and optsize[1] != imsize[1])):
3060 self.msg(str(imsize)+' is not an acceptable imagesize, will use '+str(optsize)+" instead",priority="warn")
3061 imsize=optsize
3063 if not interactive:
3064 interactive=False
3065 # print clean inputs no matter what, so user can use them.
3066 # and write a clean.last file
3067 cleanlast=open(imagename+".clean.last",'w')
3068 cleanlast.write('taskname = "clean"\n')
3070 #self.msg("clean inputs:")
3071 if self.verbose: self.msg(" ")
3072 if type(mstoimage)==type([]):
3073 cleanstr="clean(vis="+str(mstoimage)+",imagename='"+imagename+"'"
3074 cleanlast.write('vis = '+str(mstoimage)+'\n')
3075 else:
3076 cleanstr="clean(vis='"+str(mstoimage)+"',imagename='"+imagename+"'"
3077 cleanlast.write('vis = "'+str(mstoimage)+'"\n')
3078 cleanlast.write('imagename = "'+imagename+'"\n')
3079 cleanlast.write('outlierfile = ""\n')
3080 cleanlast.write('field = "'+sourcefieldlist+'"\n')
3081 cleanlast.write('spw = ""\n')
3082 cleanlast.write('selectdata = False\n')
3083 cleanlast.write('timerange = ""\n')
3084 cleanlast.write('uvrange = ""\n')
3085 cleanlast.write('antenna = ""\n')
3086 cleanlast.write('scan = ""\n')
3087 if nchan>1:
3088 cleanstr=cleanstr+",mode='"+chanmode+"',nchan="+str(nchan)
3089 cleanlast.write('mode = "'+chanmode+'"\n')
3090 cleanlast.write('nchan = '+str(nchan)+'\n')
3091 else:
3092 cleanlast.write('mode = "mfs"\n')
3093 cleanlast.write('nchan = -1\n')
3094 cleanlast.write('gridmode = ""\n')
3095 cleanlast.write('wprojplanes = 1\n')
3096 cleanlast.write('facets = 1\n')
3097 cleanlast.write('cfcache = "cfcache.dir"\n')
3098 cleanlast.write('painc = 360.0\n')
3099 cleanlast.write('epjtable = ""\n')
3100 #cleanstr=cleanstr+",interpolation='nearest'" # default change 20100518
3101 cleanlast.write('interpolation = "linear"\n')
3102 cleanstr=cleanstr+",niter="+str(niter)
3103 cleanlast.write('niter = '+str(niter)+'\n')
3104 cleanlast.write('gain = 0.1\n')
3105 cleanstr=cleanstr+",threshold='"+str(threshold)+"'"
3106 cleanlast.write('threshold = "'+str(threshold)+'"\n')
3107 cleanstr=cleanstr+",psfmode='"+psfmode+"'"
3108 cleanlast.write('psfmode = "'+psfmode+'"\n')
3109 if imagermode != "":
3110 cleanstr=cleanstr+",imagermode='"+imagermode+"'"
3111 cleanlast.write('imagermode = "'+imagermode+'"\n')
3112 cleanstr=cleanstr+",ftmachine='"+ftmachine+"'"
3113 cleanlast.write('ftmachine = "'+ftmachine+'"\n')
3114 cleanlast.write('mosweight = False\n')
3115 cleanlast.write('scaletype = "SAULT"\n')
3116 cleanlast.write('multiscale = []\n')
3117 cleanlast.write('negcomponent = -1\n')
3118 cleanlast.write('smallscalebias = 0.0\n')
3119 cleanlast.write('interactive = '+str(interactive)+'\n')
3120 if interactive:
3121 cleanstr=cleanstr+",interactive=True"
3122 if type(mask)==type(" "):
3123 cleanlast.write('mask = "'+mask+'"\n')
3124 cleanstr=cleanstr+",mask='"+mask+"'"
3125 else:
3126 cleanlast.write('mask = '+str(mask)+'\n')
3127 cleanstr=cleanstr+",mask="+str(mask)
3128 cleanlast.write('start = 0\n')
3129 cleanlast.write('width = 1\n')
3130 cleanlast.write('outframe = ""\n')
3131 cleanlast.write('veltype = "radio"\n')
3132 cellstr="['"+str(cell[0]['value'])+str(cell[0]['unit'])+"','"+str(cell[1]['value'])+str(cell[1]['unit'])+"']"
3133 cleanstr=cleanstr+",imsize="+str(imsize)+",cell="+cellstr+",phasecenter='"+str(imcenter)+"'"
3134 cleanlast.write('imsize = '+str(imsize)+'\n');
3135 cleanlast.write('cell = '+cellstr+'\n');
3136 cleanlast.write('phasecenter = "'+str(imcenter)+'"\n');
3137 cleanlast.write('restfreq = ""\n');
3138 if stokes != "I":
3139 cleanstr=cleanstr+",stokes='"+stokes+"'"
3140 cleanlast.write('stokes = "'+stokes+'"\n');
3141 cleanlast.write('weighting = "'+weighting+'"\n');
3142 cleanstr=cleanstr+",weighting='"+weighting+"'"
3143 if weighting == "briggs":
3144 cleanstr=cleanstr+",robust=0.5"
3145 cleanlast.write('robust = 0.5\n');
3146 robust=0.5
3147 else:
3148 cleanlast.write('robust = 0.0\n');
3149 robust=0.
3151 taper=False
3152 if len(outertaper) >0:
3153 taper=True
3154 if type(outertaper) == type([]):
3155 if len(outertaper[0])==0:
3156 taper=False
3157 if taper:
3158 uvtaper=True
3159 cleanlast.write('uvtaper = True\n');
3160 cleanlast.write('outertaper = "'+str(outertaper)+'"\n');
3161 cleanstr=cleanstr+",uvtaper=True,outertaper="+str(outertaper)+",innertaper=[]"
3162 else:
3163 uvtaper=False
3164 cleanlast.write('uvtaper = False\n');
3165 cleanlast.write('outertaper = []\n');
3166 cleanstr=cleanstr+",uvtaper=False"
3167 cleanlast.write('innertaper = []\n');
3168 if os.path.exists(modelimage):
3169 cleanstr=cleanstr+",modelimage='"+str(modelimage)+"'"
3170 cleanlast.write('modelimage = "'+str(modelimage)+'"\n');
3171 else:
3172 cleanlast.write('modelimage = ""\n');
3173 cleanlast.write("restoringbeam = ['']\n");
3174 cleanstr=cleanstr+",pbcor="+str(pbcor)
3175 cleanlast.write("pbcor = "+str(pbcor)+"\n");
3176 cleanlast.write("minpb = 0.2\n");
3177 cleanlast.write("calready = True\n");
3178 cleanlast.write('noise = ""\n');
3179 cleanlast.write('npixels = 0\n');
3180 cleanlast.write('npercycle = 100\n');
3181 cleanlast.write('cyclefactor = 1.5\n');
3182 cleanlast.write('cyclespeedup = -1\n');
3183 cleanlast.write('nterms = 1\n');
3184 cleanlast.write('reffreq = ""\n');
3185 cleanlast.write('chaniter = False\n');
3186 cleanstr=cleanstr+")"
3187 if self.verbose:
3188 self.msg(cleanstr,priority="warn",origin="simutil")
3189 else:
3190 self.msg(cleanstr,priority="info",origin="simutil")
3191 cleanlast.write("#"+cleanstr+"\n")
3192 cleanlast.close()
3194 if not dryrun:
3195 casalog.filter("ERROR")
3196 clean(vis=mstoimage, imagename=imagename, mode=chanmode,
3197 niter=niter, threshold=threshold, selectdata=False, nchan=nchan,
3198 psfmode=psfmode, imagermode=imagermode, ftmachine=ftmachine,
3199 imsize=imsize, cell=[str(cell[0]['value'])+str(cell[0]['unit']),str(cell[1]['value'])+str(cell[1]['unit'])], phasecenter=imcenter,
3200 stokes=stokes, weighting=weighting, robust=robust,
3201 interactive=interactive,
3202 uvtaper=uvtaper,outertaper=outertaper, pbcor=True, mask=mask,
3203 modelimage=modelimage)
3204 casalog.filter()
3206 del freq,nchan # something is holding onto the ms in table cache
3212 ##################################################################
3213 # image/tclean subtask
3215 def imtclean(self, mstoimage, imagename,
3216 gridder, deconvolver,
3217 cell, imsize, imdirection,
3218 interactive, niter, threshold, weighting,
3219 outertaper, pbcor, stokes,
3220 modelimage="", mask=[], dryrun=False):
3221 """
3222 Wrapper function for Radio Interferometric Image Reconstruction from
3223 input MeasurementSet using standard CASA imaging task ('tclean').
3225 Duplicates the method "imclean" but with non-deprecated task call.
3226 Selecting individual fields for imaging is not supported
3228 mstoimage parameter expects the path to a MeasurementSet,
3229 or list of MeasurementSets
3231 imagename parameter expects string for output image file
3233 imsize parameter expects a length-1 or length-2 list of integers
3235 cell parameter expects a length-2 list containing qa.quantity objects
3237 Just like imclean, does not yield return
3238 Creates a template '[imagename].tclean.last' file in addition to
3239 normal tclean task outputs
3240 """
3242 invocation_parameters = OrderedDict( )
3244 # use the first provided MS to determine channelization for output
3245 if is_array_type(mstoimage):
3246 ms0 = mstoimage[0]
3247 if len(mstoimage) == 1:
3248 mstoimage = mstoimage[0]
3249 else:
3250 ms0 = mstoimage
3252 if os.path.exists(ms0):
3253 tb.open(ms0 + "/SPECTRAL_WINDOW")
3254 if tb.nrows() > 1:
3255 self.msg("Detected more than one SPW in " + ms0,
3256 priority="info", origin="simutil")
3257 self.msg("Determining output cube parameters using " +
3258 "first SPW present in " + ms0,
3259 priority="info", origin="simutil")
3260 freq=tb.getvarcol("CHAN_FREQ")['r1']
3261 nchan=freq.size
3262 tb.done()
3263 elif dryrun:
3264 nchan=1 # duplicate imclean method
3265 self.msg("nchan > 1 is not supported for dryrun = True",
3266 priority="info", origin="simutil")
3268 if nchan == 1:
3269 chanmode = 'mfs'
3270 else:
3271 chanmode = 'cube'
3273 # next, define tclean call defaults
3275 # legacy comparison of image size input against heuristic
3276 optsize = [0,0]
3277 optsize[0] = _su.getOptimumSize(imsize[0])
3278 if len(imsize) == 1: # user expects square output images
3279 optsize[1]=optsize[0]
3280 else:
3281 optsize[1]=_su.getOptimumSize(imsize[1])
3282 if((optsize[0] != imsize[0]) or
3283 (len(imsize) != 1 and optsize[1] != imsize[1])):
3284 imsize = optsize
3285 self.msg(str(imsize)+" is not an acceptable imagesize, " +
3286 " using imsize=" + str(optsize) + " instead",
3287 priority="warn", origin="simutil")
3289 # the cell parameter expects a list of qa.quantity objects,
3290 formatted_correctly = [qa.isquantity(cell[i]) for i in range(len(cell))]
3291 assert False not in formatted_correctly, "simutil function imtclean expects cell parameter input to be comprised of quantity objects"
3293 # convert the first two elements for storage in the tclean.last file
3294 cellparam = [str(cell[0]['value']) + str(cell[0]['unit']),
3295 str(cell[1]['value']) + str(cell[1]['unit'])]
3297 if os.path.exists(modelimage):
3298 pass
3299 elif len(modelimage) > 0:
3300 modelimage = ""
3301 self.msg("Could not find modelimage, proceeding without one",
3302 priority="warn", origin="simutil")
3304 # set tclean top-level parameters (no parent nodes)
3305 invocation_parameters['vis'] = mstoimage
3306 invocation_parameters['selectdata'] = False
3307 invocation_parameters['imagename'] = imagename
3308 invocation_parameters['imsize'] = imsize
3309 invocation_parameters['cell'] = cellparam
3310 invocation_parameters['phasecenter'] = imdirection
3311 invocation_parameters['stokes'] = stokes
3312 invocation_parameters['startmodel'] = modelimage
3313 invocation_parameters['specmode'] = chanmode
3314 invocation_parameters['gridder'] = gridder
3315 invocation_parameters['deconvolver'] = deconvolver
3316 invocation_parameters['restoration'] = True
3317 invocation_parameters['outlierfile'] = ''
3318 invocation_parameters['weighting'] = weighting
3319 invocation_parameters['niter'] = niter
3320 invocation_parameters['usemask'] = 'user'
3321 invocation_parameters['fastnoise'] = True
3322 invocation_parameters['restart'] = True
3323 invocation_parameters['savemodel'] = 'none'
3324 invocation_parameters['calcres'] = True
3325 invocation_parameters['calcpsf'] = True
3326 invocation_parameters['parallel'] = False
3328 # subparameters
3329 invocation_parameters['restoringbeam'] = 'common'
3330 invocation_parameters['pbcor'] = pbcor
3331 invocation_parameters['uvtaper'] = outertaper
3332 if niter > 0:
3333 invocation_parameters['threshold'] = threshold
3334 invocation_parameters['interactive'] = interactive
3335 else:
3336 invocation_parameters['threshold'] = ''
3337 invocation_parameters['interactive'] = False
3338 invocation_parameters['mask'] = mask
3339 invocation_parameters['pbmask'] = 0.0
3341 # write the tclean.last file (template in case of dryrun=True)
3342 #
3343 # it would be preferable to have a helper function do _this_,
3344 # then just call tclean right from simanalyze, but precedent.
3346 filename = os.path.join(os.getcwd(),imagename+'.tclean.last')
3347 if os.path.isfile(filename):
3348 self.msg("Overwriting existing 'tclean.last' file",
3349 priority="info",origin="simutil")
3350 os.remove(filename)
3351 else:
3352 with open(filename, 'w'): pass
3354 # fill in the tclean.last file
3355 #
3356 # the sections of code that handle this for tasks are
3357 # autogenerated during the CASAshell build process. See:
3358 # [unpacked_build]/lib/py/lib/python3.6/site-packages/casashell/private
3360 with open(filename,'w') as f:
3361 # fill in the used parameters first
3362 for i in invocation_parameters:
3363 # catch None type objects returned by repr function
3364 if i.startswith('<') and s.endswith('>'):
3365 f.write("{:<20}={!s}\n".format(i, None))
3366 else:
3367 f.write("{:<20}={!r}\n".format(i,
3368 invocation_parameters[i]))
3369 # next, open and fill the task call
3370 f.write("#tclean( ")
3371 count = 0
3372 for i in invocation_parameters:
3373 if i.startswith('<') and s.endswith('>'):
3374 f.write("{!s}={!s}".format(i, None))
3375 else:
3376 f.write("{!s}={!r}".format(i,
3377 invocation_parameters[i]))
3378 count += 1
3379 if count < len(invocation_parameters): f.write(",")
3380 # finally close the task call
3381 f.write(" )\n")
3383 # gather tclean task call by attempting to parse the file just created
3384 with open(filename, 'r') as _f:
3385 for line in _f:
3386 if line.startswith('#tclean( '):
3387 task_call = line[1:-1]
3388 if self.verbose:
3389 self.msg(task_call, priority="warn", origin="simutil")
3390 else:
3391 self.msg(task_call, priority="info", origin="simutil")
3393 # now that the tclean call is specified, it may be executed
3394 if not dryrun:
3395 casalog.filter("ERROR")
3396 tclean( vis = invocation_parameters['vis'],
3397 selectdata = invocation_parameters['selectdata'],
3398 imagename=invocation_parameters['imagename'],
3399 imsize=invocation_parameters['imsize'],
3400 cell=invocation_parameters['cell'],
3401 phasecenter=invocation_parameters['phasecenter'],
3402 stokes=invocation_parameters['stokes'],
3403 startmodel=invocation_parameters['startmodel'],
3404 specmode=invocation_parameters['specmode'],
3405 gridder=invocation_parameters['gridder'],
3406 deconvolver=invocation_parameters['deconvolver'],
3407 restoration=invocation_parameters['restoration'],
3408 outlierfile=invocation_parameters['outlierfile'],
3409 weighting=invocation_parameters['weighting'],
3410 niter=invocation_parameters['niter'],
3411 usemask=invocation_parameters['usemask'],
3412 fastnoise=invocation_parameters['fastnoise'],
3413 restart=invocation_parameters['restart'],
3414 savemodel=invocation_parameters['savemodel'],
3415 calcres=invocation_parameters['calcres'],
3416 calcpsf=invocation_parameters['calcpsf'],
3417 parallel=invocation_parameters['parallel'],
3418 restoringbeam=invocation_parameters['restoringbeam'],
3419 pbcor=invocation_parameters['pbcor'],
3420 uvtaper=invocation_parameters['uvtaper'],
3421 threshold=invocation_parameters['threshold'],
3422 interactive=invocation_parameters['interactive'],
3423 mask=invocation_parameters['mask'],
3424 pbmask=invocation_parameters['pbmask'] )
3425 casalog.filter()
3431 ###################################################
3433 def flatimage(self,image,complist="",verbose=False):
3434 # flat output
3435 ia.open(image)
3436 imsize=ia.shape()
3437 imcsys=ia.coordsys()
3438 ia.done()
3439 spectax=imcsys.findcoordinate('spectral')['pixel']
3440 nchan=imsize[spectax]
3441 stokesax=imcsys.findcoordinate('stokes')['pixel']
3442 nstokes=imsize[stokesax]
3444 flat=image+".flat"
3445 if nchan>1:
3446 if verbose: self.msg("creating moment zero image "+flat,origin="flatimage")
3447 ia.open(image)
3448 flat_ia = ia.moments(moments=[-1],outfile=flat,overwrite=True)
3449 ia.done()
3450 flat_ia.close()
3451 del flat_ia
3452 else:
3453 if verbose: self.msg("removing degenerate image axes in "+flat,origin="flatimage")
3454 # just remove degenerate axes from image
3455 flat_ia = ia.newimagefromimage(infile=image,outfile=flat,dropdeg=True,overwrite=True)
3456 flat_ia.close()
3458 # seems no way to just drop the spectral and keep the stokes.
3459 if nstokes<=1:
3460 os.rename(flat,flat+".tmp")
3461 ia.open(flat+".tmp")
3462 flat_ia = ia.adddegaxes(outfile=flat,stokes='I',overwrite=True)
3463 ia.done()
3464 flat_ia.close()
3465 shutil.rmtree(flat+".tmp")
3466 del flat_ia
3468 if nstokes>1:
3469 os.rename(flat,flat+".tmp")
3470 po.open(flat+".tmp")
3471 foo=po.stokesi(outfile=flat)
3472 foo.done()
3473 po.done()
3474 shutil.rmtree(flat+".tmp")
3476 imcsys.done()
3477 del imcsys
3479 # add components
3480 if len(complist)>0:
3481 ia.open(flat)
3482 if not os.path.exists(complist):
3483 self.msg("sky component list "+str(complist)+" not found in flatimage",priority="error")
3484 cl.open(complist)
3485 ia.modify(cl.torecord(),subtract=False)
3486 cl.done()
3487 ia.done()
3492 ###################################################
3494 def convimage(self,modelflat,outflat,complist=""):
3495 # regrid flat input to flat output shape and convolve
3496 modelregrid = modelflat+".regrid"
3497 # get outflatcoordsys from outflat
3498 ia.open(outflat)
3499 outflatcs=ia.coordsys()
3500 outflatshape=ia.shape()
3501 # and beam TODO is beam the same in flat as a cube?
3502 beam=ia.restoringbeam()
3503 ia.done()
3505 ia.open(modelflat)
3506 modelflatcs=ia.coordsys()
3507 modelflatshape=ia.shape()
3508 ia.setrestoringbeam(beam=beam)
3509 tmpxx=ia.regrid(outfile=modelregrid+'.tmp', overwrite=True,
3510 csys=outflatcs.torecord(),shape=outflatshape, asvelocity=False)
3511 # ia.regrid assumes a surface brightness, or more accurately doesnt
3512 # pay attention to units at all, so we now have to scale
3513 # by the pixel size to have the right values in jy/pixel,
3515 # get pixel size from model image coordsys
3516 tmpxx.done()
3517 increments=outflatcs.increment(type="direction")['numeric']
3518 incellx=qa.quantity(abs(increments[0]),outflatcs.units(type="direction")[0])
3519 incelly=qa.quantity(abs(increments[1]),outflatcs.units(type="direction")[1])
3520 xform=outflatcs.lineartransform(type="direction")
3521 offdiag=max(abs(xform[0,1]),abs(xform[1,0]))
3522 if offdiag > 1e-4:
3523 self.msg("Your image is rotated with respect to Lat/Lon. I can't cope with that yet",priority="error")
3524 incellx=qa.mul(incellx,abs(xform[0,0]))
3525 incelly=qa.mul(incelly,abs(xform[1,1]))
3526 model_cell = [qa.convert(incellx,'arcsec'),qa.convert(incelly,'arcsec')]
3528 # and from outflat (the clean image)
3529 increments=outflatcs.increment(type="direction")['numeric']
3530 incellx=qa.quantity(abs(increments[0]),outflatcs.units(type="direction")[0])
3531 incelly=qa.quantity(abs(increments[1]),outflatcs.units(type="direction")[1])
3532 xform=outflatcs.lineartransform(type="direction")
3533 offdiag=max(abs(xform[0,1]),abs(xform[1,0]))
3534 if offdiag > 1e-4:
3535 self.msg("Your image is rotated with respect to Lat/Lon. I can't cope with that yet",priority="error")
3536 incellx=qa.mul(incellx,abs(xform[0,0]))
3537 incelly=qa.mul(incelly,abs(xform[1,1]))
3538 cell = [qa.convert(incellx,'arcsec'),qa.convert(incelly,'arcsec')]
3540 # image scaling
3541 factor = (qa.convert(cell[0],"arcsec")['value'])
3542 factor *= (qa.convert(cell[1],"arcsec")['value'])
3543 factor /= (qa.convert(model_cell[0],"arcsec")['value'])
3544 factor /= (qa.convert(model_cell[1],"arcsec")['value'])
3545 imrr = ia.imagecalc(modelregrid,
3546 "'%s'*%g" % (modelregrid+'.tmp',factor),
3547 overwrite = True)
3548 shutil.rmtree(modelregrid+".tmp")
3549 if self.verbose:
3550 self.msg("scaling model by pixel area ratio %g" % factor,origin='convimage')
3552 # add unresolved components in Jy/pix
3553 # don't do this if you've already done it in flatimage()!
3554 if (os.path.exists(complist)):
3555 cl.open(complist)
3556 imrr.modify(cl.torecord(),subtract=False)
3557 cl.done()
3559 imrr.done()
3560 ia.done()
3561 del imrr
3563 # Convolve model with beam.
3564 convolved = modelregrid + '.conv'
3565 ia.open(modelregrid)
3566 ia.setbrightnessunit("Jy/pixel")
3567 tmpcnv=ia.convolve2d(convolved,major=beam['major'],minor=beam['minor'],
3568 pa=beam['positionangle'],overwrite=True)
3569 ia.done()
3571 #tmpcnv.open(convolved)
3572 tmpcnv.setbrightnessunit("Jy/beam")
3573 tmpcnv.setrestoringbeam(beam=beam)
3574 tmpcnv.done()
3577 def bandname(self, freq):
3578 """
3579 Given freq in GHz, return the silly and in some cases deliberately
3580 confusing band name that some people insist on using.
3581 """
3582 # TODO: add a trivia argument to optionally provide the historical
3583 # radar band info from Wikipedia.
3584 band = ''
3585 if freq > 90: # Use the ALMA system.
3586 band = 'band%.0f' % (0.01 * freq) # () are mandatory!
3587 # Now switch to radar band names. Above 40 GHz different groups used
3588 # different names. Wikipedia uses ones from Baytron, a now defunct company
3589 # that made test equipment.
3590 elif freq > 75: # used as a visual sensor for experimental autonomous vehicles
3591 band = 'W'
3592 elif freq > 50: # Very strongly absorbed by atmospheric O2,
3593 band = 'V' # which resonates at 60 GHz.
3594 elif freq >= 40:
3595 band = 'Q'
3596 elif freq >= 26.5: # mapping, short range, airport surveillance;
3597 band = 'Ka' # frequency just above K band (hence 'a')
3598 # Photo radar operates at 34.300 +/- 0.100 GHz.
3599 elif freq >= 18:
3600 band = 'K' # from German kurz, meaning 'short'; limited use due to
3601 # absorption by water vapor, so Ku and Ka were used
3602 # instead for surveillance. Used for detecting clouds
3603 # and at 24.150 +/- 0.100 GHz for speeding motorists.
3604 elif freq >= 12:
3605 band = 'U' # or Ku
3606 elif freq >= 8: # Missile guidance, weather, medium-resolution mapping and ground
3607 band = 'X' # surveillance; in the USA the narrow range 10.525 GHz +/- 25 MHz
3608 # is used for airport radar and short range tracking. Named X
3609 # band because the frequency was a secret during WW2.
3610 elif freq >= 4:
3611 band = 'C' # Satellite transponders; a compromise (hence 'C') between X
3612 # and S bands; weather; long range tracking
3613 elif freq >= 2:
3614 band = 'S' # Moderate range surveillance, air traffic control,
3615 # long-range weather; 'S' for 'short'
3616 elif freq >= 1:
3617 band = 'L' # Long range air traffic control and surveillance; 'L' for 'long'
3618 elif freq >= 0.3:
3619 band = 'UHF'
3620 else:
3621 band = 'P' # for 'previous', applied retrospectively to early radar systems
3622 # Includes HF and VHF. Worse, it leaves no space for me
3623 # to put in rock band easter eggs.
3624 return band
3627 def polsettings(self, telescope, poltype=None, listall=False):
3628 """
3629 Returns stokes parameters (for use as stokes in sm.setspwindow)
3630 and feed type (for use as mode in sm.setfeed).
3632 If poltype is not specified or recognized, a guess is made using
3633 telescope. Defaults to 'XX YY', 'perfect X Y'
3635 If listall is True, return the options instead.
3636 """
3637 psets = {'circular': ('RR LL', 'perfect R L'),
3638 'linear': ('XX YY', 'perfect X Y'),
3639 'RR': ('RR', 'perfect R L')}
3640 if listall:
3641 return psets
3642 if poltype not in psets:
3643 poltype = 'linear'
3644 for circ in ('VLA', 'DRAO'): # Done this way to
3645 if circ in telescope.upper(): # catch EVLA.
3646 poltype = 'circular'
3647 return psets[poltype]
3649#######################################
3650# ALMA calcpointings
3652 def applyRotate(self, x, y, tcos, tsin):
3653 return tcos*x-tsin*y, tsin*x+tcos*y
3655 def isEven(self, num):
3656 return (num % 2) == 0
3658 # this was the only algorithm in Cycle 0 - for Cycle 1 it was
3659 # recast as BaseTriangularTiling.getPointings(), the width>height
3660 # case statement was removed, and BaseTriangularTiling was supplemented by
3661 # ShiftTriangularTiling.getPointings()
3662 def getTrianglePoints(self, width, height, angle, spacing):
3663 tcos = pl.cos(angle*pl.pi/180)
3664 tsin = pl.sin(angle*pl.pi/180)
3666 xx=[]
3667 yy=[]
3669 if (width >= height):
3670 wSpacing = spacing
3671 hSpacing = (pl.sqrt(3.) / 2) * spacing
3673 nwEven = int(pl.floor((width / 2) / wSpacing))
3674 nwOdd = int(pl.floor((width / 2) / wSpacing + 0.5))
3675 nh = int(pl.floor((height / 2) / hSpacing))
3677 for ih in pl.array(range(nh*2+1))-nh:
3678 if (self.isEven(ih)):
3679 for iw in pl.array(range(nwEven*2+1))-nwEven:
3680 x,y = self.applyRotate(iw*wSpacing, ih*hSpacing, tcos, tsin)
3681 xx.append(x)
3682 yy.append(y)
3683 else:
3684 for iw in pl.array(range(nwOdd*2+1))-nwOdd:
3685 x,y = self.applyRotate((iw+0.5)*wSpacing, ih*hSpacing, tcos, tsin)
3686 xx.append(x)
3687 yy.append(y)
3688 else:
3689 wSpacing = (pl.sqrt(3) / 2) * spacing
3690 hSpacing = spacing
3692 nw = int(pl.floor((width / 2) / wSpacing))
3693 nhEven = int(pl.floor((height / 2) / hSpacing))
3694 nhOdd = int(pl.floor((height / 2) / hSpacing + 0.5))
3696 for iw in pl.array(range(nw*2+1))-nw:
3697 if (self.isEven(iw)):
3698 for ih in pl.array(range(nhEven*2+1))-nhEven:
3699 x,y = self.applyRotate(iw*wSpacing, ih*hSpacing, tcos, tsin)
3700 xx.append(x)
3701 yy.append(y)
3702 else:
3703 for ih in pl.array(range(nhOdd*2+1))-nhOdd:
3704 x,y = self.applyRotate(iw*wSpacing, (ih+0.5)*hSpacing, tcos, tsin)
3705 xx.append(x)
3706 yy.append(y)
3707 return xx,yy
3711 def getTriangularTiling(self, longlength, latlength, angle, spacing, pb):
3713 # OT if isLandscape, shortside=Latlength
3714 # else isLandscape=false, shortside=Longlength
3716 if longlength > latlength: # OT isLandscape=True (OT uses >= )
3717 width=longlength # arcsec
3718 height=latlength # arcsec
3719 shortside=height
3720 else:
3721 width=latlength
3722 height=longlength
3723 shortside=height
3724 angle=angle+90
3726 # which tiling? Base or Shifted (OT getTiling)
3727 vSpacing = (pl.sqrt(3) / 2) * spacing
3728 n = pl.ceil(shortside / vSpacing)
3730 if n % 2 == 0:
3731 return self.getShiftTriangularTiling(width, height, angle, spacing, pb)
3732 else:
3733 return self.getBaseTriangularTiling(width, height, angle, spacing, pb)
3735 def needsFiller(self, length, spacing, bs, npoints):
3736 if length > spacing * (npoints - 1) + bs:
3737 return 1
3738 else:
3739 return 0
3741 def getBaseTriangularTiling(self, width, height, angle, spacing, pb):
3742 tcos = pl.cos(angle*pl.pi/180)
3743 tsin = pl.sin(angle*pl.pi/180)
3745 xx=[]
3746 yy=[]
3748 wSpacing = spacing
3749 hSpacing = (pl.sqrt(3.) / 2) * spacing
3751 nwEven = int(pl.floor((width / 2) / wSpacing))
3752 nwEven += self.needsFiller(width, wSpacing, pb, nwEven*2+1)
3754 nwOdd = int(pl.floor((width / 2) / wSpacing + 0.5))
3755 nwOdd += self.needsFiller(width, wSpacing, pb, nwOdd*2)
3757 nh = int(pl.floor((height / 2) / hSpacing))
3758 nh += self.needsFiller(height, hSpacing, pb, nh*2+1)
3760 for ih in pl.array(range(nh*2+1))-nh:
3761 if (self.isEven(ih)):
3762 for iw in pl.array(range(nwEven*2+1))-nwEven:
3763 x,y = self.applyRotate(iw*wSpacing, ih*hSpacing, tcos, tsin)
3764 xx.append(x)
3765 yy.append(-y) # will require additional testing @ angle>0
3766 else:
3767 for iw in pl.array(range(nwOdd*2))-nwOdd:
3768 x,y = self.applyRotate((iw+0.5)*wSpacing, ih*hSpacing, tcos, tsin)
3769 xx.append(x)
3770 yy.append(-y)
3772 return xx,yy
3777 def getShiftTriangularTiling(self, width, height, angle, spacing, pb):
3778 tcos = pl.cos(angle*pl.pi/180)
3779 tsin = pl.sin(angle*pl.pi/180)
3781 xx=[]
3782 yy=[]
3784 wSpacing = spacing
3785 hSpacing = (pl.sqrt(3.) / 2) * spacing
3787 nwEven = int(pl.floor((width / 2) / wSpacing + 0.5))
3788 nwEven += self.needsFiller(width, wSpacing, pb, nwEven*2)
3790 nwOdd = int(pl.floor((width / 2) / wSpacing ))
3791 nwOdd += self.needsFiller(width, wSpacing, pb, nwOdd*2+1)
3793 nh = int(pl.floor((height - hSpacing) / 2 / hSpacing +1))
3794 nh += self.needsFiller(height, hSpacing, pb, nh*2)
3796 for ih in pl.array(range(nh*2))-nh:
3797 if (self.isEven(ih)):
3798 for iw in pl.array(range(nwEven*2))-nwEven:
3799 x,y = self.applyRotate((iw+0.5)*wSpacing, (ih+0.5)*hSpacing, tcos, tsin)
3800 xx.append(x)
3801 yy.append(-y)
3802 else:
3803 for iw in pl.array(range(nwOdd*2+1))-nwOdd:
3804 x,y = self.applyRotate(iw*wSpacing, (ih+0.5)*hSpacing, tcos, tsin)
3805 xx.append(x)
3806 yy.append(-y)
3808 return xx,yy
3813######################################
3814 # adapted from aU.getBaselineStats
3815 def baselineLengths(self, configfile):
3816 stnx, stny, stnz, stnd, padnames, antnames, telescopename, posobs = self.readantenna(configfile)
3818 # use mean position, not official COFA=posobs
3819 cx=pl.mean(stnx)
3820 cy=pl.mean(stny)
3821 cz=pl.mean(stnz)
3822 nAntennas=len(stnx)
3823 lat,lon,el = self.itrf2loc(stnx,stny,stnz,cx,cy,cz)
3825 #l = {}
3826 #neighborlist = []
3827 maxlength = 0
3828 minlength = 1e9
3829 #mylengths = pl.zeros([nAntennas,nAntennas])
3830 mylengths=pl.zeros(nAntennas*(nAntennas-1)//2)
3831 k=0
3833 for i in range(nAntennas):
3834 for j in range(i+1,nAntennas):
3835 x = lat[i]-lat[j]
3836 y = lon[i]-lon[j]
3837 z = el[i]- el[j]
3838 #mylengths[i][j] = (x**2 + y**2 + z**2)**0.5
3839 mylengths[k]=(x**2 + y**2 + z**2)**0.5
3840 k=k+1
3842 return mylengths
3845######################
3846 def approxBeam(self,configfile,freq):
3847 # freq must be in GHz
3848 mylengths=self.baselineLengths(configfile)
3849 rmslength = pl.sqrt(pl.mean(mylengths.flatten()**2))
3850 ninety = scoreatpercentile(mylengths, 90)
3852 return 0.2997924/freq/ninety*3600.*180/pl.pi # lambda/b converted to arcsec
3855######################
3856 def sfBeam1d(self,beam="", cell="", convsupport=-1, sampling=""):
3857 """
3858 Calculates PSF of image generated by gridfunction 'SF'.
3859 Migrated codes from gjinc.sfBeam() in analysisUtil module
3860 by Todd Hunter.
3861 Note: this is simplified version of the function.
3863 Paramters:
3864 beam : Antenna primary beam (quantity)
3865 cell : Cell size of image (quantity)
3866 convsupport : convolution support used for imaging (integer)
3867 sampling : pointing spacing of observation (quantity).
3868 If not defined, comvolution of sampling kernel is
3869 skipped with warning.
3870 Returns:
3871 Estimated PSF of image (quantity).
3872 """
3874 if not qa.compare(beam, "rad"):
3875 raise ValueError("beam should be a quantity of antenna primary beam size (angle)")
3876 if not qa.compare(cell, "rad"):
3877 raise ValueError("cell should be a quantity of image pixel size (angle)")
3878 if len(sampling) > 0 and not qa.compare(sampling, "rad"):
3879 raise ValueError("sampling should be a quantity of pointing spacing (angle)")
3880 if (convsupport < -1):
3881 convsupport = 3
3883 supportwidth = (convsupport*2 + 1)
3884 c = supportwidth * pl.pi/2.
3885 pb_asec = qa.getvalue(qa.convert(beam, "arcsec"))
3886 cell_asec = qa.getvalue(qa.convert(cell, "arcsec"))
3887 samp_asec = -1.
3888 if len(sampling) > 0:
3889 samp_asec = qa.getvalue(qa.convert(sampling, "arcsec"))
3890 # Define kernel array
3891 myxaxis = pl.arange(-3*pb_asec,3*pb_asec+0.1,0.2)
3892 # Gaussian func of PB
3893 scale = 0.5 / ( (pb_asec/2.3548201)**2 )
3894 mygaussian = pl.exp(- scale * myxaxis**2 ) ### exp(x**2/(2*sigma**2))
3895 # Spheroidal kernel func
3896 sfaxis = myxaxis/float((supportwidth-1)*cell_asec/2.0)
3897 indices = pl.where(abs(sfaxis)<1)[0]
3898 centralRegion = sfaxis[indices]
3899 centralRegionY = spspec.pro_ang1_cv(0, 0, c, spspec.pro_cv(0,0,c),
3900 centralRegion)[0]
3901 mysf = pl.zeros(len(myxaxis))
3902 mysf[indices] += centralRegionY/max(centralRegionY)
3903 # Convolution of Gaussian PB and SF
3904 result = spsig.convolve(mysf, mygaussian, mode='same')
3905 del mygaussian, sfaxis, indices, centralRegion, centralRegionY
3906 # Sampling function
3907 if samp_asec > 0:
3908 myboxcar = pl.zeros(len(myxaxis))
3909 indices = pl.where(abs(myxaxis)<samp_asec/2.)[0]
3910 myboxcar[indices] = 1
3911 result = spsig.convolve(result, myboxcar, mode='same')
3912 else:
3913 self.msg("Pointing spacing is not specified. Calculated PSF could be less accurate.", priority="warn", origin="sfBeam1d")
3914 # Calculate FWHM
3915 result /= max(result)
3916 halfmax = max(result)*0.5
3917 spline = spintrp.UnivariateSpline(myxaxis, result-halfmax, s=0)
3918 x0, x1 = spline.roots()
3919 del result, myxaxis, spline
3920 return qa.quantity(abs(x1-x0), "arcsec")