Source code for pipeline.extern.adopted
import os
import xml.dom.minidom as minidom
import numpy as np
import pipeline.infrastructure.logging as logging
from pipeline.infrastructure import casa_tools
LOG = logging.get_logger(__name__)
qa = casa_tools.quanta
[docs]def getMedianPWV(vis='.', myTimes=[0, 999999999999], asdm='', verbose=False):
"""
Extracts the PWV measurements from the WVR on all antennas for the
specified time range. The time range is input as a two-element list of
MJD seconds (default = all times). First, it tries to find the ASDM_CALWVR
table in the ms. If that fails, it then tries to find the
ASDM_CALATMOSPHERE table in the ms. If that fails, it then tried to find
the CalWVR.xml in the specified ASDM, or failing that, an ASDM of the
same name (-.ms). If neither of these exist, then it tries to find
CalWVR.xml in the present working directory. If it still fails, it looks
for CalWVR.xml in the .ms directory. Thus, you only need to copy this
xml file from the ASDM into your ms, rather than the entire ASDM. Returns
the median and standard deviation in millimeters.
Returns:
The median PWV, and the median absolute deviation (scaled to match rms)
For further help and examples, see https://safe.nrao.edu/wiki/bin/view/ALMA/GetMedianPWV
-- Todd Hunter
"""
tb = casa_tools.table
pwvmean = 0
success = False
if verbose:
print("in getMedianPWV with myTimes = {}".format(myTimes))
try:
if os.path.exists("%s/ASDM_CALWVR" % vis):
tb.open("%s/ASDM_CALWVR" % vis)
pwvtime = tb.getcol('startValidTime') # mjdsec
antenna = tb.getcol('antennaName')
pwv = tb.getcol('water')
tb.close()
success = True
if len(pwv) < 1:
if os.path.exists("%s/ASDM_CALATMOSPHERE" % vis):
pwvtime, antenna, pwv = readPWVFromASDM_CALATMOSPHERE(vis)
success = True
if len(pwv) < 1:
print("Found no data in ASDM_CALWVR nor ASDM_CALATMOSPHERE table")
return 0, -1
else:
if verbose:
print("Did not find ASDM_CALATMOSPHERE in the ms")
return 0, -1
if verbose:
print("Opened ASDM_CALWVR table, len(pwvtime)={}".format(len(pwvtime)))
else:
if verbose:
print("Did not find ASDM_CALWVR table in the ms. Will look for ASDM_CALATMOSPHERE next.")
if os.path.exists("%s/ASDM_CALATMOSPHERE" % vis):
pwvtime, antenna, pwv = readPWVFromASDM_CALATMOSPHERE(vis)
success = True
if len(pwv) < 1:
print("Found no data in ASDM_CALATMOSPHERE table")
return 0, -1
else:
if verbose:
print("Did not find ASDM_CALATMOSPHERE in the ms")
except:
if verbose:
print("Could not open ASDM_CALWVR table in the ms")
finally:
# try to find the ASDM table
if not success:
if len(asdm) > 0:
if not os.path.exists(asdm):
print("Could not open ASDM = {}".format(asdm))
return 0, -1
try:
[pwvtime, pwv, antenna] = readpwv(asdm)
except:
if verbose:
print("Could not open ASDM = {}".format(asdm))
return pwvmean, -1
else:
try:
tryasdm = vis.split('.ms')[0]
if verbose:
print("No ASDM name provided, so I will try this name = {}".format(tryasdm))
[pwvtime, pwv, antenna] = readpwv(tryasdm)
except:
try:
if verbose:
print("Still did not find it. Will look for CalWVR.xml in current directory.")
[pwvtime, pwv, antenna] = readpwv('.')
except:
try:
if verbose:
print("Still did not find it. Will look for CalWVR.xml in the .ms directory.")
[pwvtime, pwv, antenna] = readpwv('%s/' % vis)
except:
if verbose:
print("No CalWVR.xml file found, so no PWV retrieved. Copy it to this directory and try"
" again.")
return pwvmean, -1
try:
matches = np.where(np.array(pwvtime) > myTimes[0])[0]
except:
print("Found no times > {:d}".format(myTimes[0]))
return 0, -1
if len(pwv) < 1:
print("Found no PWV data")
return 0, -1
if verbose:
print("{:d} matches = {}".format(len(matches), matches))
print("{:d} pwv = {}".format(len(pwv), pwv))
ptime = np.array(pwvtime)[matches]
matchedpwv = np.array(pwv)[matches]
matches2 = np.where(ptime <= myTimes[-1])[0]
if verbose:
print("matchedpwv = {}".format(matchedpwv))
print("pwv = {}".format(pwv))
if len(matches2) < 1:
# look for the value with the closest start time
mindiff = 1e12
for i in range(len(pwvtime)):
if abs(myTimes[0] - pwvtime[i]) < mindiff:
mindiff = abs(myTimes[0] - pwvtime[i])
matchedpwv = []
for i in range(len(pwvtime)):
if abs(abs(myTimes[0]-pwvtime[i]) - mindiff) < 1.0:
matchedpwv.append(pwv[i])
pwvmean = 1000 * np.median(matchedpwv)
if verbose:
print("Taking the median of {:d} pwv measurements from all antennas = {:.3f} mm".format(len(matchedpwv),
pwvmean))
pwvstd = 1000 * MAD(matchedpwv)
else:
pwvmean = 1000 * np.median(matchedpwv[matches2])
pwvstd = 1000 * MAD(matchedpwv[matches2])
if verbose:
print("Taking the median of {:d} pwv measurements from all antennas = {:.3f} mm".format(len(matches2),
pwvmean))
return pwvmean, pwvstd
[docs]def readPWVFromASDM_CALATMOSPHERE(vis):
"""
Reads the PWV via the water column of the ASDM_CALATMOSPHERE table.
- Todd Hunter
"""
if not os.path.exists(vis + '/ASDM_CALATMOSPHERE'):
if vis.find('.ms') < 0:
vis += '.ms'
if not os.path.exists(vis):
print("Could not find measurement set = {}".format(vis))
return
elif not os.path.exists(vis + '/ASDM_CALATMOSPHERE'):
print("Could not find ASDM_CALATMOSPHERE in the measurement set")
return
else:
print("Could not find measurement set")
return
mytb = casa_tools.table
mytb.open("%s/ASDM_CALATMOSPHERE" % vis)
pwvtime = mytb.getcol('startValidTime') # mjdsec
antenna = mytb.getcol('antennaName')
pwv = mytb.getcol('water')[0] # There seem to be 2 identical entries per row, so take first one.
mytb.close()
return pwvtime, antenna, pwv
[docs]def readpwv(asdm):
"""
This function assembles the dictionary returned by readwvr() into arrays
containing the PWV measurements written by TelCal into the ASDM.
Units are in meters.
-- Todd Hunter
"""
wvrdict = readwvr(asdm)
bigantlist = []
for entry in wvrdict:
bigantlist.append(wvrdict[entry]['antenna'])
watertime = []
water = []
antenna = []
for entry in wvrdict:
measurements = 1
for i in range(measurements):
watertime.append(wvrdict[entry]['startmjdsec'] + (i * 1.0 / measurements) * wvrdict[entry]['duration'])
water.append(wvrdict[entry]['water'])
antenna.append(wvrdict[entry]['antenna'])
return [watertime, water, antenna]
[docs]def readwvr(sdmfile, verbose=False):
"""
This function reads the CalWVR.xml table from the ASDM and returns a
dictionary containing: 'start', 'end', 'startmjd', 'endmjd',
'startmjdsec', 'endmjdsec',
'timerange', 'antenna', 'water', 'duration'.
'water' is the zenith PWV in meters.
This function is called by readpwv(). -- Todd Hunter
"""
if not os.path.exists(sdmfile):
print("readwvr(): Could not find file = {}".format(sdmfile))
return
xmlscans = minidom.parse(sdmfile + '/CalWVR.xml')
scandict = {}
rowlist = xmlscans.getElementsByTagName("row")
fid = 0
for rownode in rowlist:
rowpwv = rownode.getElementsByTagName("water")
pwv = float(rowpwv[0].childNodes[0].nodeValue)
water = pwv
scandict[fid] = {}
# start and end times in mjd ns
rowstart = rownode.getElementsByTagName("startValidTime")
start = int(rowstart[0].childNodes[0].nodeValue)
startmjd = float(start) * 1.0E-9 / 86400.0
t = qa.quantity(startmjd, 'd')
starttime = call_qa_time(t, form="ymd", prec=8)
rowend = rownode.getElementsByTagName("endValidTime")
end = int(rowend[0].childNodes[0].nodeValue)
endmjd = float(end) * 1.0E-9 / 86400.0
t = qa.quantity(endmjd, 'd')
endtime = call_qa_time(t, form="ymd", prec=8)
# antenna
rowantenna = rownode.getElementsByTagName("antennaName")
antenna = str(rowantenna[0].childNodes[0].nodeValue)
scandict[fid]['start'] = starttime
scandict[fid]['end'] = endtime
scandict[fid]['startmjd'] = startmjd
scandict[fid]['endmjd'] = endmjd
scandict[fid]['startmjdsec'] = startmjd * 86400
scandict[fid]['endmjdsec'] = endmjd * 86400
timestr = starttime + '~' + endtime
scandict[fid]['timerange'] = timestr
scandict[fid]['antenna'] = antenna
scandict[fid]['water'] = water
scandict[fid]['duration'] = (endmjd - startmjd) * 86400
fid += 1
if verbose:
print(' Found {} rows in CalWVR.xml'.format(rowlist.length))
# return the dictionary for later use
return scandict
[docs]def call_qa_time(arg, form='', prec=0, showform=False):
"""
This is a wrapper for qa.time(), which in casa 4.0.0 returns a list
of strings instead of just a scalar string.
- Todd Hunter
"""
if isinstance(arg, dict) and isinstance(arg['value'], (list, np.ndarray)) and len(arg['value']) > 1:
LOG.warning('qa_time() received a dictionary containing a list'
'of length=%d rather than a scalar. Using first '
'value.', len(arg['value']))
arg['value'] = arg['value'][0]
result = qa.time(arg, form=form, prec=prec, showform=showform)
if type(result) in (list, np.ndarray):
return result[0]
else:
return result
[docs]def MAD(a, c=0.6745, axis=0):
"""
Median Absolute Deviation along given axis of an array:
median(abs(a - median(a))) / c
c = 0.6745 is the constant to convert from MAD to std; it is used by
default
"""
a = np.array(a)
good = (a == a)
a = np.asarray(a, np.float64)
if a.ndim == 1:
d = np.median(a[good])
m = np.median(np.fabs(a[good] - d) / c)
else:
d = np.median(a[good], axis=axis)
# I don't want the array to change so I have to copy it?
if axis > 0:
aswp = np.swapaxes(a[good], 0, axis)
else:
aswp = a[good]
m = np.median(np.fabs(aswp - d) / c, axis=0)
return m
[docs]def getMJD():
"""
Returns the current MJD. See also getCurrentMJDSec().
-Todd
"""
myme = casa_tools.measures
mjd = myme.epoch('utc', 'today')['m0']['value']
myme.done()
return mjd