Source
xxxxxxxxxx
17
17
import inspect
18
18
import math
19
19
import os
20
20
import re # used for testing if a string is a float
21
21
import time
22
22
23
23
import matplotlib.transforms
24
24
import numpy as np
25
25
import pylab as pb
26
26
from casatasks import casalog
27
+
from casatasks.private import sdutil, simutil
27
28
from casatools import (atmosphere, ctsys, measures, ms, msmetadata, quanta,
28
29
table)
29
30
from matplotlib.ticker import (FormatStrFormatter, MultipleLocator,
30
31
ScalarFormatter)
31
32
from six.moves import input, range
32
33
33
34
# CAS-13722, CAS-13385
34
35
import warnings
35
36
import matplotlib.cbook
36
37
warnings.filterwarnings("ignore",category=matplotlib.cbook.MatplotlibDeprecationWarning)
5517
5518
else:
5518
5519
try:
5519
5520
pwvmean = float(pwv)
5520
5521
except:
5521
5522
pwvmean = defaultPWV
5522
5523
5523
5524
if (verbose):
5524
5525
print("Using PWV = %.2f mm" % pwvmean)
5525
5526
5526
5527
# default values in case we can't find them below
5528
+
myqa = quanta()
5527
5529
airmass = 1.5
5528
-
P = 563.0
5530
+
P = myqa.quantity(563.0, 'mbar')
5529
5531
H = 20.0
5530
5532
T = 273.0
5531
5533
roundedScanTimes = []
5532
5534
if (type(mymsmd) != str):
5533
5535
if (verbose):
5534
5536
print("Looking for scans for field integer = %d" % (field))
5535
5537
scans = mymsmd.scansforfield(field)
5536
5538
if (verbose):
5537
5539
print(("For field %s, Got scans = " % str(field),scans))
5538
5540
for myscan in scans:
5552
5554
# if (verbose): print("tdiff = %s" % (str(tdiff)))
5553
5555
if (tdiff < mindiff):
5554
5556
bestscan = scans[i]
5555
5557
if (verbose): print("bestscan = %s" % (str(bestscan)))
5556
5558
mindiff = tdiff
5557
5559
if (verbose):
5558
5560
print("For timestamp=%.1f, got closest scan = %d, %.0f sec away" %(timestamp, bestscan,mindiff))
5559
5561
if (verbose): print("Calling getWeather()")
5560
5562
[conditions,myTimes] = getWeather(vis,bestscan,antenna,verbose,mymsmd)
5561
5563
if (verbose): print("Done getWeather()")
5562
-
P = conditions['pressure']
5564
+
P = myqa.quantity(conditions['pressure'], conditions['pressure_unit'])
5563
5565
H = conditions['humidity']
5564
5566
T = conditions['temperature']+273.15
5565
-
if (P <= 0.0):
5566
-
P = 563
5567
+
if (P['value'] <= 0.0):
5568
+
P = myqa.quantity(563, 'mbar')
5567
5569
if (H <= 0.0):
5568
5570
H = 20
5569
5571
else:
5570
5572
conditions = {}
5571
5573
if (type(mymsmd) != str):
5572
5574
if ('elevation' not in list(conditions.keys())):
5573
5575
# Someone cleared the POINTING table, so calculate elevation from Ra/Dec/MJD
5574
5576
# myfieldId = mymsmd.fieldsforname(mymsmd.fieldsforscan(bestscan))
5575
5577
myfieldId = mymsmd.fieldsforscan(bestscan)[0]
5576
5578
myscantime = np.mean(mymsmd.timesforscan(bestscan))
5586
5588
myazel = computeAzElFromRADecMJD(mydirection, myscantime/86400., telescopeName)
5587
5589
conditions['elevation'] = myazel[1] * 180/math.pi
5588
5590
conditions['azimuth'] = myazel[0] * 180/math.pi
5589
5591
if (verbose):
5590
5592
print("Computed elevation = %.1f deg" % (conditions['elevation']))
5591
5593
else:
5592
5594
conditions['elevation'] = 90
5593
5595
bestscan = -1
5594
5596
if (verbose):
5595
5597
print("CalcAtm: found elevation=%f (airmass=%.3f) for scan: %s" % (conditions['elevation'],1/np.sin(conditions['elevation']*np.pi/180.), str(bestscan)))
5596
-
print("P,H,T = %f,%f,%f" % (P,H,T))
5598
+
print("P,H,T = %f,%f,%f" % (P['value'],H,T))
5597
5599
if (conditions['elevation'] <= 3):
5598
5600
print("Using 45 deg elevation instead")
5599
5601
airmass = 1.0/math.cos(45*math.pi/180.)
5600
5602
else:
5601
5603
airmass = 1.0/math.cos((90-conditions['elevation'])*math.pi/180.)
5602
5604
5605
+
geodetic_elevation = 5059
5606
+
with sdutil.table_manager(os.path.join(vis, 'ANTENNA')) as tb:
5607
+
_X, _Y, _Z = (float(i) for i in tb.getcell('POSITION', antenna))
5608
+
_pos = simutil.simutil().xyz2long(_X, _Y, _Z, 'WGS84')
5609
+
geodetic_elevation = _pos[2]
5610
+
5603
5611
tropical = 1
5604
5612
midLatitudeSummer = 2
5605
5613
midLatitudeWinter = 3
5606
5614
numchan = len(freqs)
5607
5615
# Set the reference freq to be the middle of the middle two channels
5608
5616
reffreq=0.5*(freqs[numchan//2-1]+freqs[numchan//2])
5609
5617
originalnumchan = numchan
5610
5618
while (numchan > MAX_ATM_CALC_CHANNELS):
5611
5619
numchan //= 2
5612
5620
# print("Reducing numchan to ", numchan)
5613
5621
chans = list(range(0,originalnumchan,(originalnumchan//numchan)))
5614
5622
5615
5623
chansep = (freqs[-1]-freqs[0])/(numchan-1)
5616
5624
nbands = 1
5617
5625
if (verbose): print("Opening casac.atmosphere()")
5618
5626
myat = atmosphere()
5619
5627
if (verbose): print("Opened")
5620
-
myqa = quanta()
5621
5628
fCenter = myqa.quantity(reffreq,'GHz')
5622
5629
fResolution = myqa.quantity(chansep,'GHz')
5623
5630
fWidth = myqa.quantity(numchan*chansep,'GHz')
5624
-
myat.initAtmProfile(humidity=H,temperature=myqa.quantity(T,"K"),altitude=myqa.quantity(5059,"m"),pressure=myqa.quantity(P,'mbar'),atmType=midLatitudeWinter)
5631
+
myat.initAtmProfile(humidity=H, temperature=myqa.quantity(T, "K"),
5632
+
altitude=myqa.quantity(geodetic_elevation, "m"),
5633
+
pressure=myqa.convert(P, 'mbar'), atmType=midLatitudeWinter)
5625
5634
myat.initSpectralWindow(nbands,fCenter,fWidth,fResolution)
5626
5635
myat.setUserWH2O(myqa.quantity(pwvmean,'mm'))
5627
5636
5628
5637
# myat.setAirMass() # This does not affect the opacity, but it does effect TebbSky, so do it manually.
5629
5638
5630
5639
n = myat.getNumChan()
5631
5640
if (verbose): print("numchan = %s" % (str(n)))
5632
5641
if ctsys.compare_version('<',[4,0,0]):
5633
5642
dry = np.array(myat.getDryOpacitySpec(0)['dryOpacity'])
5634
5643
wet = np.array(myat.getWetOpacitySpec(0)['wetOpacity'].value)
7023
7032
mytb.open("%s/WEATHER" % vis)
7024
7033
except:
7025
7034
print("Could not open the WEATHER table for this ms.")
7026
7035
mytb.done()
7027
7036
return([conditions,myTimes])
7028
7037
if (True):
7029
7038
mjdsec = mytb.getcol('TIME')
7030
7039
indices = np.argsort(mjdsec)
7031
7040
mjd = mjdsec/86400.
7032
7041
pressure = mytb.getcol('PRESSURE')
7042
+
conditions['pressure_unit'] = mytb.getcolkeyword('PRESSURE', 'QuantumUnits')[0]
7033
7043
relativeHumidity = mytb.getcol('REL_HUMIDITY')
7034
7044
temperature = mytb.getcol('TEMPERATURE')
7035
7045
if (np.median(temperature) > 100):
7036
7046
# must be in units of Kelvin, so convert to C
7037
7047
temperature -= 273.15
7038
7048
if 'DEW_POINT' in mytb.colnames():
7039
7049
dewPoint = mytb.getcol('DEW_POINT')
7040
7050
if (np.median(dewPoint) > 100):
7041
7051
# must be in units of Kelvin, so convert to C
7042
7052
dewPoint -= 273.15
7604
7614
agrees with the path.
7605
7615
"""
7606
7616
mypwd = os.getcwd() + '/'
7607
7617
newfilelist = []
7608
7618
for f in filelist:
7609
7619
fstart = 0
7610
7620
if (f.find(mypwd) == 0):
7611
7621
fstart = len(mypwd)
7612
7622
newfilelist.append(f[fstart:])
7613
7623
return(newfilelist)
7614
-