Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/simutil.py: 61%

2174 statements  

« prev     ^ index     » next       coverage.py v7.6.4, created at 2024-11-01 07:19 +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 

14 

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( ) 

29 

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 

37 

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) 

53 

54 return ( out.rstrip(sep) + ")" ) 

55 

56def __get_str(paramval): 

57 if type(paramval) == str: 

58 return ("'%s'" % paramval) 

59 # else 

60 return str(paramval) 

61 

62 

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] 

91 

92 

93 

94 

95 

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 

118 

119 

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() 

126 

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) 

135 

136 

137 ########################################################### 

138 

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 

151 

152 ########################################################### 

153 

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() 

165 

166 

167 

168 ########################################################### 

169 

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 

177 

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 

188 

189 clr="" 

190 if self.verbose: 

191 toterm=True 

192 else: 

193 toterm=False 

194 

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" 

213 

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") 

220 

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] 

227 

228 if origin: 

229 casalog.post(clr+"["+origin+"] "+bw+s) 

230 else: 

231 casalog.post(s) 

232 

233 

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) 

240 

241 

242 ########################################################### 

243 

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 

253 

254 

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") 

264 

265 

266 def closereport(self): 

267 self.report.close() 

268 

269 

270 

271 

272 

273 ########################################################### 

274 

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 

286 

287 ########################################################### 

288 

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 

306 

307 ########################################################### 

308 

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 

336 

337 

338 

339 ########################################################### 

340 # plot an image (optionally), and calculate its statistics 

341 

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() 

388 

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. 

398 

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 

439 

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 

461 

462 

463 

464 

465 

466 

467 

468 ########################################################### 

469 

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. 

474  

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) 

487 

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) 

496 

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] 

503 

504 

505 # haveing eliminated other options, we need to calculate: 

506 epoch, centx, centy = self.direction_splitter() 

507 

508 pointings = [] 

509 

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) 

518 

519 xsize=qa.quantity(size[0]) 

520 ysize=qa.quantity(size[1]) 

521 

522 nrows = 1+ int(pl.floor(qa.convert(qa.div(ysize, yspacing), '')['value'] 

523 - 2.309401077 * relmargin)) 

524 

525 availcols = 1 + qa.convert(qa.div(xsize, spacing), 

526 '')['value'] - 2.0 * relmargin 

527 ncols = int(pl.floor(availcols)) 

528 

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 = [] 

538 

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) 

544 

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 

562 

563 xsize=qa.quantity(size[0]) 

564 ysize=qa.quantity(size[1]) 

565 

566 nrows = 1+ int(pl.floor(qa.convert(qa.div(ysize, yspacing), '')['value'] 

567 - 2.0 * relmargin)) 

568 

569 availcols = 1 + qa.convert(qa.div(xsize, spacing), '')['value'] \ 

570 - 2.0 * relmargin 

571 ncols = int(pl.floor(availcols)) 

572 

573 

574 ncolstomin = 0.5 * (ncols - 1) # O O O O 

575 pointings = [] # O O O O 

576 

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) 

582 

583 xmin = qa.sub(centx, qa.mul(ncolstomin, xspacing)) 

584 stopcolp1 = ncols 

585 

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]) 

597 

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 

602 

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) 

609 

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

619 

620 

621 # if could not fit any pointings, then return single pointing 

622 if(len(pointings)==0): 

623 pointings.append(direction) 

624 

625 self.msg("using %i generated pointing(s)" % len(pointings),origin='calc_pointings') 

626 self.pointings=pointings 

627 return pointings 

628 

629 

630 

631 

632 

633 

634 

635 

636 

637 

638 ########################################################### 

639 

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 

646  

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.  

654  

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 

660  

661 """ 

662 f=open(filename) 

663 line= ' ' 

664 time=[] 

665 pointings=[] 

666 

667 # add option of different epoch in a header line like read_antenna? 

668 

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() 

705 

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 

713 

714 return len(pointings), pointings, time 

715 

716 

717 

718 

719 

720 

721 

722 ########################################################### 

723 

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

728  

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 

734  

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

745 

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") 

754 

755 f.close() 

756 return 

757 

758 

759 

760 ########################################################### 

761 

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 

805 

806 

807 

808 ########################################################### 

809 

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 

857 

858 

859 

860 ########################################################### 

861 

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) 

891 

892 return epoch, qa.convert(x, 'deg'), qa.convert(y, 'deg') 

893 

894 

895 ########################################################### 

896 

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

915 

916 

917 ########################################################### 

918 

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) 

929 

930 ########################################################### 

931 

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 

943 

944 

945 

946 

947 

948 

949 

950 

951 

952 

953 ########################################################### 

954 #========================== tsys ========================== 

955 

956 def noisetemp(self, telescope=None, freq=None, 

957 diam=None, epsilon=None): 

958 

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

976 

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) 

982 

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 

987 

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 

993 

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 

999 

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 ;) 

1005 

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") 

1009 

1010 # blockage efficiency.  

1011 eta_b = 1.-(diam_subreflector/diam)**2 

1012 

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 

1019 

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") 

1029 

1030 # antenna efficiency 

1031 # eta_a = eta_p*eta_s*eta_b*eta_t 

1032 

1033 # correlator quantization efficiency.  

1034 eta_q = cq[iobs] 

1035 

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/ 

1040 

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] 

1057 

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] 

1069 

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] 

1099 

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

1103 

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] 

1115 

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") 

1123 

1124 return eta_p, eta_s, eta_b, eta_t, eta_q, t_rx 

1125 

1126 

1127 

1128 

1129 

1130 

1131 

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): 

1138 

1139 if qa.convert(elevation,"deg")['value']<3: 

1140 self.msg("Elevation < ALMA limit of 3 deg",priority="error") 

1141 return False 

1142 

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]) 

1155 

1156 msfile=tmpname+".ms" 

1157 sm.open(msfile) 

1158 

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 

1170 

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 

1182 

1183 if found_match == True or found_partial_match == True: 

1184 known = True 

1185 

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 

1192 

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] 

1201 

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) 

1215 

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 

1225 

1226 nant=len(stnx) 

1227 # diam is only used as a test below, not quantitatively 

1228 diam = pl.average(stnd) 

1229 

1230 

1231 if (telescope==None or diam==None): 

1232 self.msg("Telescope name or antenna diameter have not been set.",priority="error") 

1233 return False 

1234 

1235 # copied from task_simdata: 

1236 

1237 self.setcfg(sm, telescope, stnx, stny, stnz, stnd, padnames, posobs) 

1238 

1239 model_nchan=1 

1240 # RI TODO isquantity checks 

1241 model_width=qa.quantity(bandwidth) # note: ATM uses band center 

1242 

1243 # start is center of first channel. for nch=1, that equals center 

1244 model_start=qa.quantity(freq) 

1245 

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=['']) 

1252 

1253 sm.setlimits(shadowlimit=0.01, elevationlimit='3deg') 

1254 sm.setauto(0.0) 

1255 

1256 obslat=qa.convert(posobs['m1'],'deg') 

1257 dec=qa.add(obslat, qa.add(qa.quantity("90deg"),qa.mul(elevation,-1))) 

1258 

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) 

1268 

1269 sm.observe(sourcename="src1", spwname="band1", 

1270 starttime=qa.quantity(0, "s"), 

1271 stoptime=qa.quantity(etime)); 

1272 

1273 sm.setdata() 

1274 sm.setvp() 

1275 

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") 

1282 

1283 if pwv==None: 

1284 # RI TODO choose based on freq octile 

1285 pwv=2.0 

1286 

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 

1292 

1293 if telescope=='ALMA' and (qa.convert(freq,"GHz")['value'])>600.: 

1294 rxtype=1 # DSB 

1295 

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 

1313 

1314 if doimnoise: 

1315 sm.corrupt() 

1316 

1317 sm.done() 

1318 

1319 

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") 

1328 

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. 

1335 

1336 nint = qa.convert(etime,'s')['value'] / qa.convert(integration,'s')['value'] 

1337 nbase= 0.5*nant*(nant-1) 

1338 

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. 

1354 

1355 theoreticalnoise=noiseperbase/pl.sqrt(nint*nbase*2) # assume 2-poln 

1356 

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]) 

1364 

1365 if doimnoise: 

1366 return theoreticalnoise , imnoise 

1367 else: 

1368 return theoreticalnoise 

1369 

1370 

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. 

1384 

1385 Optional: 

1386 mounttype (Defaults to a guess based on telescope.) 

1387 

1388 Returns the mounttype that it uses. 

1389 

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

1406 

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 

1414 

1415 

1416 

1417 

1418 ########################################################### 

1419 #===================== ephemeris ========================== 

1420 

1421 

1422 def ephemeris(self, date, usehourangle=True, direction=None, telescope=None, ms=None, cofa=None): 

1423 

1424 if direction==None: direction=self.direction 

1425 if telescope==None: telescope=self.telescopename 

1426 

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" 

1434 

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. 

1438 

1439 ds=self.direction_splitter(direction) # if list, returns average 

1440 src=me.direction(ds[0],ds[1],ds[2]) 

1441 

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) 

1449 

1450 time=me.epoch('TAI',date) 

1451 me.doframe(time) 

1452 

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) 

1461 

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] 

1465 

1466 timeinc='15min' # for plotting 

1467 timeinc=qa.convert(qa.time(timeinc)[0],'d')['value'] 

1468 ntime=int(1./timeinc) 

1469 

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') 

1484 

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 

1491 

1492 times=[] 

1493 az=[] 

1494 el=[] 

1495 

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 

1503 

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') 

1507 

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') 

1512 

1513 relpeak=peaktime_float-reftime_floor 

1514 pl.plot(pl.array([1,1])*24*relpeak,[0,90]) 

1515 

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') 

1530 

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") 

1538 

1539 

1540 pl.ylim([0,90]) 

1541 pl.xlim(pl.array([-12,12])+24*(reftime_float-reftime_floor)) 

1542 

1543 

1544 

1545 

1546 

1547 

1548 

1549 

1550 

1551 

1552 

1553 

1554 

1555 

1556 

1557 

1558 

1559 

1560 

1561 ########################################################### 

1562 #========================================================== 

1563 

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 

1586  

1587 returns: earth-centered x,y,z, diameter, pad_name, antenna_name, observatory_name, observatory_measure_dictionary 

1588 

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 

1599 

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') 

1652 

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"] 

1660 

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") 

1668 

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 

1679 

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') 

1684 

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={} 

1702 

1703 

1704 

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 

1737 

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) 

1763 

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")) 

1774 

1775 return (stnx, stny, stnz, pl.array(ind), pid, aid, self.telescopename, posobs) 

1776 

1777 

1778 

1779 

1780 

1781 

1782 

1783 

1784 

1785 

1786 

1787 

1788 

1789 

1790 

1791 

1792 

1793 

1794 ########################################################### 

1795 #==================== geodesy ============================= 

1796 

1797 

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. 

1808  

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 

1826  

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

1834 

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 

1856 

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 

1864 

1865 tn=sinfi/cosfi 

1866 ts=tn*tn 

1867 

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

1873 

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

1878 

1879 return lat,lon,conv,kp 

1880 

1881 

1882 

1883 

1884 def tconpc(self,sf,orlim,er,esq,rf): 

1885 

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 

1899  

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

1909 

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 

1917 

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 

1926 

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 

1935 

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 

1940 

1941 return eps,r,so,v0,v2,v4,v6 

1942 

1943 

1944 

1945 

1946 

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 

1953  

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 

1958 

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' ]} 

1973 

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' ]} 

1999 

2000 if datumcode not in datums: 

2001 self.msg("unknown datum %s" % datumcode,priority="error") 

2002 return -1 

2003 

2004 datum=datums[datumcode] 

2005 ellipsoid=datum[3] 

2006 

2007 if ellipsoid not in ellipsoids: 

2008 self.msg("unknown ellipsoid %s" % ellipsoid,priority="error") 

2009 return -1 

2010 

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] 

2014 

2015 

2016 

2017 

2018 

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) 

2023 

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 

2027  

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

2074  

2075 this is the driver program to compute latitudes and longitudes from 

2076 the utms for each zone 

2077  

2078 input: 

2079 northing, easting 

2080 zone, datum 

2081 nors=N/S 

2082  

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 

2088  

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) 

2099  

2100 """ 

2101 

2102 rad=180./pl.pi 

2103 

2104 offx,offy,offz,er,rf = self.getdatum(datum,verbose=self.verbose) 

2105 

2106 f=1./rf 

2107 esq=(2*f-f*f) 

2108 

2109 # find the central meridian if the zone number is less than 30 

2110 

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 

2123 

2124 tol=(5./60.)/rad 

2125 

2126 if nors == 'S': 

2127 fn= 10000000. 

2128 else: 

2129 fn=0. 

2130 

2131 fe=500000.0 

2132 sf=0.9996 

2133 orlim=0.0 

2134 

2135 found=0 

2136 

2137 # precompute parameters for this zone: 

2138 eps,r,so,v0,v2,v4,v6 = self.tconpc(sf,orlim,er,esq,rf) 

2139 

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) 

2142 

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 

2146 

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 

2162 

2163 # *** convert to more usual convention of negative lon = W 

2164 lon=-lon 

2165 

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") 

2168 

2169 return lon,lat 

2170 

2171 

2172 

2173 

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

2183 

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) 

2187 

2188 f=1./rf 

2189 esq=2*f-f**2 

2190 nu=er/pl.sqrt(1.-esq*(pl.sin(lat))**2) 

2191 

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 

2195 

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 

2199 

2200 return x,y,z 

2201 

2202 

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 

2217 

2218 dx,dy,dz,er,rf = self.getdatum(datum,verbose=False) 

2219 f=1./rf 

2220 

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 

2224 

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 

2232 

2233 # Distance from polar axis 

2234 r = pl.sqrt( (x-dx)**2 + (y-dy)**2 ) 

2235 

2236 # Longitude. 

2237 if r>0.: 

2238 lon = pl.arctan2(y-dx, x-dx) 

2239 else: 

2240 lon = 0. 

2241 

2242 # Prepare Newton correction factors. 

2243 s0 = abs(z-dz) / er 

2244 c0 = ec * r / er 

2245 

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 

2249 

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) 

2254 

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) 

2258 

2259 # Restore sign of latitude.  

2260 if (z-dz) < 0: lat = -lat 

2261 

2262 return lon,lat,height 

2263 

2264 

2265 def xyz2long_old(self,x,y,z,datum): 

2266 

2267 dx,dy,dz,er,rf = self.getdatum(datum,verbose=False) 

2268 

2269 f=1./rf 

2270 

2271 b= ((x-dx)**2 + (y-dy)**2) / er**2 

2272 c= (z-dx)**2 / er**2 

2273 

2274 esq=2*f-f**2 # (a2-b2)/a2 

2275 

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 

2282 

2283 b0=4.*a0 

2284 b1=3.*a1 

2285 b2=2.*a2 

2286 b3=a3 

2287 

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) 

2296 

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 

2302 

2303 return lon,lat 

2304 

2305 

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

2315 

2316 lon,lat = self.utm2long(easting,northing,zone,datum,nors) 

2317 x,y,z = self.long2xyz(lon,lat,elevation,datum) 

2318 

2319 return x,y,z 

2320 

2321 

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) 

2340 

2341 # Now you see the connection between the Old Ones and Antarctica... 

2342 Nploczcphimlocysphi = (N + locz+alt) * pl.cos(phi) - locy * sphi 

2343 

2344 clmb = pl.cos(lmbda) 

2345 slmb = pl.sin(lmbda) 

2346 

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) 

2350 

2351 return x, y, z 

2352 

2353 

2354 

2355 

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) 

2370 

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) 

2379 

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 

2390 

2391 return lat,lon,el 

2392 

2393 

2394 

2395 

2396 

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']) 

2408 

2409 return self.itrf2loc(x,y,z,cx,cy,cz) 

2410 

2411 

2412 

2413 

2414 

2415 

2416 

2417 

2418 

2419 

2420 

2421 

2422 

2423 ########################################################### 

2424 

2425 def plotants(self,x,y,z,d,name): 

2426 # given globals 

2427 

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__() 

2433 

2434 dolam=0 

2435 # TODO convert to klam: (d too) 

2436 ### 

2437 

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) 

2454 

2455 pl.axis('equal') 

2456 #if dolam: 

2457 # pl.xlabel("kilolamda") 

2458 # pl.ylabel("kilolamda") 

2459 

2460 

2461 

2462 

2463 

2464 

2465 

2466 

2467 

2468 

2469 

2470 

2471 

2472 

2473 

2474 

2475 

2476 

2477 

2478 ###################################################### 

2479 # helper function to get the pixel size from an image (positive increments) 

2480 

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] 

2495 

2496 ###################################################### 

2497 # helper function to get the spectral size from an image 

2498 

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 

2510 

2511 ###################################################### 

2512 

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 

2533 

2534 ################################################################## 

2535 # fit modelimage into a 4 coordinate image defined by the parameters 

2536 

2537 # TODO spectral extrapolation and regridding using innchan **** 

2538 

2539 def modifymodel(self, inimage, outimage, 

2540 inbright,indirection,incell,incenter,inwidth,innchan, 

2541 flatimage=False): # if nonzero, create mom -1 image  

2542 

2543 # new ia tool 

2544 in_ia=ia.newimagefromfile(inimage) 

2545 in_shape=in_ia.shape() 

2546 in_csys=in_ia.coordsys() 

2547 

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] 

2557 

2558 

2559 

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 

2580 

2581 # check shape characteristics of the input; 

2582 # add degenerate axes as neeed: 

2583 

2584 in_dir=in_csys.findcoordinate("direction") 

2585 in_spc=in_csys.findcoordinate("spectral") 

2586 in_stk=in_csys.findcoordinate("stokes") 

2587 

2588 

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 

2607 

2608 

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() 

2622 

2623 

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'] 

2651 

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]) 

2662 

2663 # preserve sign in case input image is backwards (RA increasing right) 

2664 model_cell = [qa.convert(incellx,'arcsec'),qa.convert(incelly,'arcsec')] 

2665 

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 

2671 

2672 

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) 

2684 

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 

2696 

2697 

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 

2714 

2715 

2716 

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") 

2720 

2721 

2722 

2723 

2724 

2725 

2726 

2727 # we've now found or assigned two direction axes, and changed direction and cell if required 

2728 # next, work on spectral axis: 

2729 

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 

2756 

2757 

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 

2769 

2770 axmap[3]=extra_axis 

2771 axassigned[extra_axis]=3 

2772 model_nchan=arr.shape[extra_axis] 

2773 

2774 

2775 

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 

2795 

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 

2813 

2814 

2815 

2816 

2817 

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 

2845 

2846 

2847 

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 

2861 

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" 

2877 

2878 

2879 

2880 

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") 

2888 

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") 

2899 

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") 

2916 

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") 

2921 

2922 

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 

2930 

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 

2938 

2939 if outimage!=inimage: 

2940 ia.setcoordsys(modelcsys.torecord()) 

2941 ia.done() 

2942 ia.open(outimage) 

2943 

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 

2952 

2953 

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 

2962 

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") 

2969 

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 

2976 

2977 # make a moment 0 image 

2978 if flatimage and outimage!=inimage: 

2979 self.flatimage(outimage,verbose=self.verbose) 

2980 

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])] 

2985 

2986 if self.verbose: self.msg(" ") # add a line after my spewage 

2987 

2988 return model_refdir,model_cell,model_size,model_nchan,model_specrefval,model_specrefpix,model_width,model_stokes 

2989 

2990 

2991 

2992 

2993 

2994 ################################################################## 

2995 # image/clean subtask 

2996 

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' 

3009 

3010 No return 

3011 

3012 Creates a template '[imagename].clean.last' file in addition to  

3013 outputs of task 'clean' 

3014 """ 

3015 

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 

3023 

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 

3033 

3034 if nchan==1: 

3035 chanmode="mfs" 

3036 else: 

3037 chanmode="channel" 

3038 

3039 psfmode="clark" 

3040 ftmachine="ft" 

3041 imagermode="clark" # set default to prevent UnboundLocalError 

3042 

3043 if cleanmode=="csclean": 

3044 imagermode='csclean' 

3045 #if cleanmode=="clark": 

3046 # imagermode="" 

3047 if cleanmode=="mosaic": 

3048 imagermode="mosaic" 

3049 ftmachine="mosaic" 

3050 

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 

3062 

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') 

3069 

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. 

3150 

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() 

3193 

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() 

3205 

3206 del freq,nchan # something is holding onto the ms in table cache 

3207 

3208 

3209 

3210 

3211 

3212 ################################################################## 

3213 # image/tclean subtask 

3214 

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

3224 

3225 Duplicates the method "imclean" but with non-deprecated task call. 

3226 Selecting individual fields for imaging is not supported 

3227 

3228 mstoimage parameter expects the path to a MeasurementSet,  

3229 or list of MeasurementSets 

3230  

3231 imagename parameter expects string for output image file 

3232 

3233 imsize parameter expects a length-1 or length-2 list of integers 

3234 

3235 cell parameter expects a length-2 list containing qa.quantity objects 

3236 

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

3241 

3242 invocation_parameters = OrderedDict( ) 

3243 

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 

3251 

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") 

3267 

3268 if nchan == 1: 

3269 chanmode = 'mfs' 

3270 else: 

3271 chanmode = 'cube' 

3272 

3273 # next, define tclean call defaults 

3274 

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") 

3288 

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" 

3292 

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'])] 

3296 

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") 

3303 

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 

3327 

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 

3340 

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. 

3345 

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 

3353 

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 

3359 

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") 

3382 

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") 

3392 

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() 

3426 

3427 

3428 

3429 

3430 

3431 ################################################### 

3432 

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] 

3443 

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() 

3457 

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 

3467 

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") 

3475 

3476 imcsys.done() 

3477 del imcsys 

3478 

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() 

3488 

3489 

3490 

3491 

3492 ################################################### 

3493 

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() 

3504 

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,  

3514 

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')] 

3527 

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')] 

3539 

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') 

3551 

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() 

3558 

3559 imrr.done() 

3560 ia.done() 

3561 del imrr 

3562 

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() 

3570 

3571 #tmpcnv.open(convolved) 

3572 tmpcnv.setbrightnessunit("Jy/beam") 

3573 tmpcnv.setrestoringbeam(beam=beam) 

3574 tmpcnv.done() 

3575 

3576 

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 

3625 

3626 

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

3631 

3632 If poltype is not specified or recognized, a guess is made using 

3633 telescope. Defaults to 'XX YY', 'perfect X Y' 

3634 

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] 

3648 

3649####################################### 

3650# ALMA calcpointings 

3651 

3652 def applyRotate(self, x, y, tcos, tsin): 

3653 return tcos*x-tsin*y, tsin*x+tcos*y 

3654 

3655 def isEven(self, num): 

3656 return (num % 2) == 0 

3657 

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) 

3665 

3666 xx=[] 

3667 yy=[] 

3668 

3669 if (width >= height): 

3670 wSpacing = spacing 

3671 hSpacing = (pl.sqrt(3.) / 2) * spacing 

3672 

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

3676 

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 

3691 

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

3695 

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 

3708 

3709 

3710 

3711 def getTriangularTiling(self, longlength, latlength, angle, spacing, pb): 

3712 

3713 # OT if isLandscape, shortside=Latlength 

3714 # else isLandscape=false, shortside=Longlength 

3715 

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 

3725 

3726 # which tiling? Base or Shifted (OT getTiling) 

3727 vSpacing = (pl.sqrt(3) / 2) * spacing 

3728 n = pl.ceil(shortside / vSpacing) 

3729 

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) 

3734 

3735 def needsFiller(self, length, spacing, bs, npoints): 

3736 if length > spacing * (npoints - 1) + bs: 

3737 return 1 

3738 else: 

3739 return 0 

3740 

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) 

3744 

3745 xx=[] 

3746 yy=[] 

3747 

3748 wSpacing = spacing 

3749 hSpacing = (pl.sqrt(3.) / 2) * spacing 

3750 

3751 nwEven = int(pl.floor((width / 2) / wSpacing)) 

3752 nwEven += self.needsFiller(width, wSpacing, pb, nwEven*2+1) 

3753 

3754 nwOdd = int(pl.floor((width / 2) / wSpacing + 0.5)) 

3755 nwOdd += self.needsFiller(width, wSpacing, pb, nwOdd*2) 

3756 

3757 nh = int(pl.floor((height / 2) / hSpacing)) 

3758 nh += self.needsFiller(height, hSpacing, pb, nh*2+1) 

3759 

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) 

3771 

3772 return xx,yy 

3773 

3774 

3775 

3776 

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) 

3780 

3781 xx=[] 

3782 yy=[] 

3783 

3784 wSpacing = spacing 

3785 hSpacing = (pl.sqrt(3.) / 2) * spacing 

3786 

3787 nwEven = int(pl.floor((width / 2) / wSpacing + 0.5)) 

3788 nwEven += self.needsFiller(width, wSpacing, pb, nwEven*2) 

3789 

3790 nwOdd = int(pl.floor((width / 2) / wSpacing )) 

3791 nwOdd += self.needsFiller(width, wSpacing, pb, nwOdd*2+1) 

3792 

3793 nh = int(pl.floor((height - hSpacing) / 2 / hSpacing +1)) 

3794 nh += self.needsFiller(height, hSpacing, pb, nh*2) 

3795 

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) 

3807 

3808 return xx,yy 

3809 

3810 

3811 

3812 

3813###################################### 

3814 # adapted from aU.getBaselineStats 

3815 def baselineLengths(self, configfile): 

3816 stnx, stny, stnz, stnd, padnames, antnames, telescopename, posobs = self.readantenna(configfile) 

3817 

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) 

3824 

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 

3832 

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 

3841 

3842 return mylengths 

3843 

3844 

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) 

3851 

3852 return 0.2997924/freq/ninety*3600.*180/pl.pi # lambda/b converted to arcsec 

3853 

3854 

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. 

3862  

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

3873 

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 

3882 

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")