Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/task_plotweather.py: 7%
273 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-10-31 18:48 +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###############################################
11from casatools import atmosphere, table, ms, quanta, measures
12from casatasks import casalog
14mytb = table( )
15myms = ms( )
16myqa = quanta( )
17myme = measures( )
19import pylab as pl
20import numpy as np
21from math import pi,floor
22import os.path as osp
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]
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
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)
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()
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()
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
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']
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'
101 # check for weather table
103 if osp.isdir(myMS+'/WEATHER'):
105 try:
106 mytb.open(myMS+'/WEATHER')
107 firstTime = mytb.getcol('TIME')[0]
108 mytb.close()
109 WEATHER_table_exists = True
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
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()
130 ##retrieve stuff from weather table, if exists
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()
146 else:
147 myms.open(myMS)
148 mytime_range = myms.range(["time"])
149 mytime = [mytime_range['time'][0]]
150 myms.close()
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)
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)
170 ##convert time to a decimal
171 numtime=pl.datestr2num(myTimestr)
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.)
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))
184 # estimate pwv from seasonal model zenith opacity
185 myPWV2 = -1.71 + 1.3647*myTauZ1
186 myPWV = (1-seasonal_weight)*myPWV1 + seasonal_weight*myPWV2
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
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
208 fC=myqa.quantity(25.0,'GHz')
209 fW=myqa.quantity(50.,'GHz')
210 fR=myqa.quantity(0.025,'GHz')
212 at=atmosphere()
213 hum=20.0
215 myatm=at.initAtmProfile(alt,tmp,pre,mxA,hum,wvl,dpr,dpm,h0,att)
217 at.initSpectralWindow(nb,fC,fW,fR)
218 sg=at.getSpectralWindow()
219 mysg = sg['value']
221 nstep = 20
222 pwv = []
223 opac = pl.zeros((len(mysg),nstep))
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
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
246 freqs=pl.array(mysg)/1e9
247 coef0=pwv_coef[:,1]/1e-3
248 coef1=pwv_coef[:,0]/1e-3
251 #interpolate between nearest table entries for each spw center frequency
252 meanTau=[]
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
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
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)
279 ##make the plots
281 if doPlot==False:
282 return meanTau
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'
294 thisfig=pl.figure(1)
295 thisfig.clf()
296 thisfig.set_size_inches(8.5,10)
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
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))
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)
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))
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))
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')
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)
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 ])
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')
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()
382 casalog.post('wrote weather figure: '+plotName)
383 return meanTau