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

273 statements  

« prev     ^ index     » next       coverage.py v7.6.4, created at 2024-11-01 07:19 +0000

1############################################### 

2## To plot stuff in weather tables, saved to MSname+.plotWX.png 

3## and estimate zenith opacity per spw, returned as a list named myTau 

4## 

5## 

6## J. Marvil 2.6.12 

7## revised 4.27.12 to add support for missing/empty weather table 

8## revised 11.05.12 to address CASA 4.0 changes 

9############################################### 

10 

11from casatools import atmosphere, table, ms, quanta, measures 

12from casatasks import casalog 

13 

14mytb = table( ) 

15myms = ms( ) 

16myqa = quanta( ) 

17myme = measures( ) 

18 

19import pylab as pl 

20import numpy as np 

21from math import pi,floor 

22import os.path as osp 

23 

24def _find(condition): 

25 """Returns indices where ravel(a) is true. 

26 Private implementation of deprecated matplotlib.mlab.find 

27 Thanks to: https://github.com/python-control/python-control/pull/262/files 

28 """ 

29 return np.nonzero(np.ravel(condition))[0] 

30 

31############### 

32## hides the extreme Y-axis ticks, helps stack plots close together without labels overlaping 

33def jm_clip_Yticks(): 

34 xa=pl.gca() 

35 nlabels=0 

36 for label in xa.yaxis.get_ticklabels(): 

37 nlabels+=1 

38 thislabel=0 

39 if nlabels>3: 

40 for label in xa.yaxis.get_ticklabels(): 

41 if thislabel==0: label.set_alpha(0) 

42 if thislabel==nlabels-1: label.set_alpha(0) 

43 thislabel+=1 

44 

45############## 

46## sets the position of the y-axis label to the right side of the plot, can also move up/down 

47def jm_set_Ylabel_pos(pos=(0.5,0.5)): 

48 ax=pl.gca(); 

49 ax.yaxis.set_label_position('right') 

50 ax.yaxis.label.set_rotation(270) 

51 ax.yaxis.label.set_position(pos) 

52 

53 

54############### 

55## fixed y-ticks, from myMin to myMax 

56def jm_set_Ylim_ticks(myMin=-1,myMax=1): 

57 myYlocs=pl.linspace(round(myMin,1),round(myMax,1),5) 

58 myLocator = pl.FixedLocator(myYlocs) 

59 ax=pl.gca() 

60 ax.yaxis.set_major_locator( myLocator ) 

61 pl.ylim(myMin,myMax) 

62 jm_clip_Yticks() 

63 

64############### 

65## variable y-ticks, but not more than 1+ this argument 

66def jm_set_Yvar_ticks(myScale=4): 

67 xa=pl.gca() 

68 xa.yaxis.set_major_locator(pl.MaxNLocator(myScale)) 

69 jm_clip_Yticks() 

70 

71############### 

72## calculates K-band zenith opacity from temperature and dewpoint 

73def Tau_K_Calc(D,T,day, weights=(.5,.5)): 

74 P = pl.exp(1.81+(17.27*D)/(D+237.3)) # water vapor partial pressure 

75 h = 324.7*P/(T+273.15) # PWV in mm 

76 tau_w = 3.8 + 0.23*h + 0.065*h**2 # tau from weather, in %, at 22GHz 

77 if day > 199: day = day - 365. 

78 m = day + 165. # modified day of the year 

79 tau_d = 22.1 - 0.178*m + 0.00044*m**2 # tau from seaonal model, in % 

80 tau_k = weights[0]*tau_w + weights[1]*tau_d # the average, with equal weights (as in the AIPS default) 

81 return tau_k, h 

82 

83################ 

84## calculates elevation of the sun 

85def jm_sunEL(mytime): 

86 myme.doframe(myme.observatory('VLA')) 

87 myme.doframe(myme.epoch('utc',mytime)) 

88 mysun=myme.measure(myme.direction('SUN'),'AZELGEO') 

89 return mysun['m1']['value'] 

90 

91################ 

92## gets and plots data from the weather table of the given MS 

93def plotweather(vis='', seasonal_weight=0.5, doPlot=True, plotName = ''): 

94 myMS=vis 

95 if plotName == '': 

96 if myMS.endswith("/"): 

97 plotName = myMS + myMS.rstrip("/") + '.plotweather.png' 

98 else: 

99 plotName = myMS + '.plotweather.png' 

100 

101 # check for weather table 

102 

103 if osp.isdir(myMS+'/WEATHER'): 

104 

105 try: 

106 mytb.open(myMS+'/WEATHER') 

107 firstTime = mytb.getcol('TIME')[0] 

108 mytb.close() 

109 WEATHER_table_exists = True 

110 

111 except: 

112 casalog.post('could not open weather table, using seasonal model only and turning off plots') 

113 mytb.close() 

114 WEATHER_table_exists = False 

115 doPlot=False 

116 seasonal_weight = 1.0 

117 else: 

118 casalog.post('could not find a weather table, using seasonal model only and turning off plots') 

119 WEATHER_table_exists = False 

120 doPlot=False 

121 seasonal_weight = 1.0 

122 

123 

124 

125 ##retrieve center frequency for each sub-band 

126 mytb.open(myMS+'/SPECTRAL_WINDOW') 

127 spwFreqs=mytb.getcol('REF_FREQUENCY') * 1e-9 

128 mytb.close() 

129 

130 ##retrieve stuff from weather table, if exists 

131 

132 if WEATHER_table_exists: 

133 mytb.open(myMS+'/WEATHER') 

134 mytime=mytb.getcol('TIME') 

135 mytemp=mytb.getcol('TEMPERATURE') - 273.15 

136 mydew=mytb.getcol('DEW_POINT') - 273.15 

137 mywinds=mytb.getcol('WIND_SPEED') 

138 # Text starts at 90 degrees, whereas the wind direction starts at 0 

139 # Hence the wind direction is adjusted 90 degrees counterclockwise 

140 # to make the arrows point to right direction 

141 mywindd=270-mytb.getcol('WIND_DIRECTION')*(180.0/pi) 

142 mypres=mytb.getcol('PRESSURE') 

143 myhum=mytb.getcol('REL_HUMIDITY') 

144 mytb.close() 

145 

146 else: 

147 myms.open(myMS) 

148 mytime_range = myms.range(["time"]) 

149 mytime = [mytime_range['time'][0]] 

150 myms.close() 

151 

152 ##calculate the elevation of the sun 

153 sunEL=[] 

154 for time in mytime: 

155 t1= myqa.quantity(time,'s') 

156 myday=myqa.convert(t1,'d') 

157 sunEL1=jm_sunEL(myday) 

158 sunEL.append(sunEL1) 

159 

160 ##convert time to a string of date/time 

161 myTimestr = [] 

162 for time in mytime: 

163 q1=myqa.quantity(time,'s') 

164 time1=myqa.time(q1,form='ymd')[0] 

165 if '24:00:00' in time1: 

166 q1=myqa.quantity(time + 0.5,'s') 

167 time1=myqa.time(q1,form='ymd')[0] 

168 myTimestr.append(time1) 

169 

170 ##convert time to a decimal 

171 numtime=pl.datestr2num(myTimestr) 

172 

173 ##### calculate opacity as in EVLA memo 143 

174 thisday= 30*(float(myTimestr[0][5:7])-1)+float(myTimestr[0][8:10]) 

175 thisday=thisday + 5 * (thisday / 365.) 

176 

177 

178 if WEATHER_table_exists: 

179 # get 22 GHz zenith opacity and pwv estimate from weatherstation (myPWV) 

180 myTauZ, myPWV1 = Tau_K_Calc(mydew,mytemp,thisday) 

181 myTauZ1, myPWV = Tau_K_Calc(mydew,mytemp,thisday, weights=(0,1.0)) 

182 myTauZ2, myPWV = Tau_K_Calc(mydew,mytemp,thisday, weights=(1.0,0)) 

183 

184 # estimate pwv from seasonal model zenith opacity 

185 myPWV2 = -1.71 + 1.3647*myTauZ1 

186 myPWV = (1-seasonal_weight)*myPWV1 + seasonal_weight*myPWV2 

187 

188 else: 

189 day = thisday*1.0 

190 if day > 199: day = day - 365. 

191 m = day + 165. # modified day of the year 

192 myTauZ = 22.1 - 0.178*m + 0.00044*m**2 # tau from seaonal model, in % 

193 myPWV = -1.71 + 1.3647*myTauZ 

194 myPWV1, myPWV2 = myPWV, myPWV 

195 myTauZ1, myTauZ2 = myTauZ, myTauZ 

196 

197 tmp = myqa.quantity(270.0,'K') 

198 pre = myqa.quantity(790.0,'mbar') 

199 alt = myqa.quantity(2125,'m') 

200 h0 = myqa.quantity(2.0,'km') 

201 wvl = myqa.quantity(-5.6, 'K/km') 

202 mxA = myqa.quantity(48,'km') 

203 dpr = myqa.quantity(10.0,'mbar') 

204 dpm = 1.2 

205 att = 1 

206 nb = 1 

207 

208 fC=myqa.quantity(25.0,'GHz') 

209 fW=myqa.quantity(50.,'GHz') 

210 fR=myqa.quantity(0.025,'GHz') 

211 

212 at=atmosphere() 

213 hum=20.0 

214 

215 myatm=at.initAtmProfile(alt,tmp,pre,mxA,hum,wvl,dpr,dpm,h0,att) 

216 

217 at.initSpectralWindow(nb,fC,fW,fR) 

218 sg=at.getSpectralWindow() 

219 mysg = sg['value'] 

220 

221 nstep = 20 

222 pwv = [] 

223 opac = pl.zeros((len(mysg),nstep)) 

224 

225 for i in range(nstep): 

226 hum = 20.0*(i+1) 

227 myatm = at.initAtmProfile(alt,tmp,pre,mxA,hum,wvl,dpr,dpm,h0,att) 

228 w=at.getGroundWH2O() 

229 pwv.append(w['value'][0]) 

230 at.initSpectralWindow(nb,fC,fW,fR) 

231 at.setUserWH2O(w) 

232 sg=at.getSpectralWindow() 

233 mysg = sg['value'] 

234 sdry=at.getDryOpacitySpec() 

235 swet=at.getWetOpacitySpec() 

236 sd=sdry[1] 

237 sw=swet[1]['value'] 

238 stot = pl.array(sd)+pl.array(sw) 

239 opac[:,i]=stot 

240 

241 pwv_coef=pl.zeros((len(mysg),2)) 

242 for i in range(len(mysg)): 

243 myfit=pl.polyfit(pwv,opac[i,:],1) 

244 pwv_coef[i,:]=myfit 

245 

246 freqs=pl.array(mysg)/1e9 

247 coef0=pwv_coef[:,1]/1e-3 

248 coef1=pwv_coef[:,0]/1e-3 

249 

250 

251 #interpolate between nearest table entries for each spw center frequency 

252 meanTau=[] 

253 

254 for i in range(len(spwFreqs)): 

255 mysearch=(pl.array(freqs)-spwFreqs[i])**2 

256 hits=_find(mysearch == min(mysearch)) 

257 # Fix deprecation warning: could be array of 1 

258 #if len(hits) > 1: hits=hits[0] 

259 if not isinstance(hits, int): 

260 hits = hits[0] 

261 tau_interp = (pl.array(coef0[hits-2:hits+2])+pl.array(coef1[hits-2:hits+2])*pl.mean(myPWV)) * 1e-1 #percent 

262 tau_F = pl.interp(spwFreqs[i],freqs[hits-2:hits+2],tau_interp) 

263 meanTau.append(pl.mean(tau_F*.01)) #nepers 

264 

265 

266 tau_allF = (pl.array(coef0) + pl.array(coef1)*pl.mean(myPWV)) * 1e-1 #percent 

267 tau_allF1 = (pl.array(coef0) + pl.array(coef1)*pl.mean(myPWV1)) *1e-1 

268 tau_allF2 = (pl.array(coef0) + pl.array(coef1)*pl.mean(myPWV2)) *1e-1 

269 

270 casalog.post('SPW : Frequency (GHz) : Zenith opacity (nepers)') 

271 for i in range(len(meanTau)): 

272 myStr = str(i).rjust(3) + ' : ' 

273 myStr2 = '%.3f'%(spwFreqs[i]) 

274 myStr += myStr2.rjust(7) + ' : ' +str(round(meanTau[i], 3)) 

275 casalog.post(myStr) 

276 

277 

278 

279 ##make the plots 

280 

281 if doPlot==False: 

282 return meanTau 

283 

284 pl.ioff() 

285 myColor2='#A6A6A6' 

286 myColorW='#92B5F2' 

287 myColor1='#4D4DFF' 

288 myOrangeColor='#FF6600' 

289 myYellowColor='#FFCC00' 

290 myWeirdColor='#006666' 

291 myLightBrown='#996633' 

292 myDarkGreay='#333333' 

293 

294 thisfig=pl.figure(1) 

295 thisfig.clf() 

296 thisfig.set_size_inches(8.5,10) 

297 

298 Xrange=numtime[-1]-numtime[0] 

299 Yrange=max(mywinds)-min(mywinds) 

300 Xtextoffset=-Xrange*.01 

301 Ytextoffset=-Yrange*.08 

302 Xplotpad=Xrange*.03 

303 Yplotpad=Yrange*.03 

304 

305 sp1=thisfig.add_axes([.13,.8,.8,.15]) 

306 pl.ylabel('solar el') 

307 nsuns=30 

308 myj=pl.array(pl.linspace(0,len(sunEL)-1,nsuns),dtype='int') 

309 for i in myj: 

310 if sunEL[i]<0: pl.plot([numtime[i],numtime[i]],[(180/pi)*sunEL[i],(180/pi)*sunEL[i]],'kH') 

311 else: pl.plot([numtime[i],numtime[i]],[(180/pi)*sunEL[i],(180/pi)*sunEL[i]],'H',color=myYellowColor) 

312 pl.plot([numtime[0],numtime[-1]],[0,0],'-',color='brown') 

313 xa=pl.gca(); xa.set_xticklabels('') 

314 jm_set_Ylim_ticks(myMin=-90,myMax=90) 

315 jm_set_Ylabel_pos(pos=(0,.5)) 

316 pl.title('Weather Summary for '+myMS) 

317 pl.xlim(numtime[0]-Xplotpad,numtime[-1]+Xplotpad) 

318 xa.set_xticks(pl.linspace(min(numtime),max(numtime),3)) 

319 

320 sp2=thisfig.add_axes([.13,.65,.8,.15]) 

321 pl.ylabel('wind (m/s)') 

322 nwind=60 

323 myj=pl.array(pl.linspace(0,len(mywinds)-1,nwind),dtype='int') 

324 for i in myj: 

325 pl.text(numtime[i]+Xtextoffset,Ytextoffset+mywinds[i],'-->',rotation=mywindd[i], alpha=1,color='purple',fontsize=12) 

326 

327 pl.plot(numtime, .3+mywinds,'.', color='black', ms=2, alpha=0) 

328 jm_set_Ylabel_pos(pos=(0,.5)) 

329 jm_set_Yvar_ticks(5) 

330 xa=pl.gca(); xa.set_xticklabels('') 

331 pl.xlim(numtime[0]-Xplotpad,numtime[-1]+Xplotpad) 

332 pl.ylim(min(mywinds)-Yplotpad,max(mywinds)+Yplotpad) 

333 xa.set_xticks(pl.linspace(min(numtime),max(numtime),3)) 

334 

335 sp4=thisfig.add_axes([.13,.5,.8,.15]) 

336 pl.plot(numtime, mytemp,'-', color=myOrangeColor,lw=2) 

337 pl.plot(numtime, mydew,'-', color=myWeirdColor,lw=2) 

338 pl.ylabel('temp,dew') 

339 jm_set_Ylabel_pos(pos=(0, .5)) 

340 xa=pl.gca(); xa.set_xticklabels('') 

341 jm_set_Yvar_ticks(5) 

342 pl.xlim(numtime[0]-Xplotpad,numtime[-1]+Xplotpad) 

343 xa.set_xticks(pl.linspace(min(numtime),max(numtime),3)) 

344 

345 sp7=thisfig.add_axes([.13,.35,.8,.15]) 

346 pl.ylabel('PWV (mm)') 

347 pl.plot(numtime, myPWV2, color=myColor2, lw=2, label='seasonal model') 

348 pl.plot(numtime, myPWV1, color=myColor1, lw=2, label='weather station') 

349 pl.plot(numtime, myPWV, color=myColorW,lw=2, label='weighted') 

350 

351 thismin=min([min(myPWV),min(myPWV1),min(myPWV2)]) 

352 thismax=max([max(myPWV),max(myPWV1),max(myPWV2)]) 

353 pl.ylim(.8*thismin,1.2*thismax) 

354 jm_set_Ylabel_pos(pos=(0,.5)) 

355 jm_set_Yvar_ticks(5) 

356 xa=pl.gca(); xa.set_xticklabels('') 

357 pl.xlim(numtime[0]-Xplotpad,numtime[-1]+Xplotpad) 

358 

359 middletimei=int(floor(len(myTimestr)/2.)) 

360 middletimes=str(myTimestr[middletimei])[11:] 

361 endtimes=myTimestr[-1][11:] 

362 ax=pl.gca() 

363 axt=ax.get_xticks() 

364 ax.set_xticks(pl.linspace(min(numtime),max(numtime),3)) 

365 ax.set_xticklabels([myTimestr[0],middletimes,endtimes ]) 

366 

367 sp8=thisfig.add_axes([.13,.1,.8,.2]) 

368 pl.plot(freqs,.01*tau_allF2,'-', color=myColor2, lw=2, label='seasonal model') 

369 pl.plot(freqs,.01*tau_allF1,'-', color=myColor1, lw=2, label='weather station') 

370 pl.plot(freqs,.01*tau_allF,'-', color=myColorW, lw=2,label='weighted') 

371 

372 

373 sp8.legend(loc=2, borderaxespad=0) 

374 pl.ylabel('Tau_Z (nepers)') 

375 pl.xlabel('Frequency (GHz)') 

376 pl.ylim(0,.25) 

377 jm_set_Yvar_ticks(6) 

378 jm_set_Ylabel_pos(pos=(0,.5)) 

379 pl.savefig( plotName, dpi=150) 

380 pl.close() 

381 

382 casalog.post('wrote weather figure: '+plotName) 

383 return meanTau