#!/usr/bin/env python
# This module has been derived from Todd Hunter's plotSpectrumFromMask module.
# It is used to create the diagnostic spectra plots for the cube imaging weblog.
import os
import re
from math import degrees
import matplotlib.pyplot as plt
import matplotlib.ticker
import numpy as np
from casatasks import imhead
import pipeline.infrastructure as infrastructure
from pipeline.infrastructure import casa_tools
LOG = infrastructure.get_logger(__name__)
ARCSEC_PER_RAD = 206264.80624709636
c_mks = 2.99792458e8
[docs]def addFrequencyAxisAbove(ax1, firstFreq, lastFreq, freqType='', spw=None,
fontsize=10, showlabel=True, twinx=False,
ylimits=None, xlimits=None):
"""
Given the descriptor returned by plt.subplot, draws a frequency
axis label at the top x-axis.
firstFreq and lastFreq: in Hz
twinx: if True, the create an independent y-axis on right edge
freqType: 'LSRK', 'REST', etc. for axis label purposes only
spw: integer or string (simply to add to the title)
xlimits: tuple or list or array of 2 values (in Hz)
"""
if showlabel:
if (spw != '' and spw is not None):
label = '(Spw %s) %s Frequency (GHz)' % (str(spw), freqType)
else:
label = '%s Frequency (GHz)' % freqType
if twinx:
ax2 = ax1.twiny().twinx()
ax2.set_ylabel('Per-channel noise (mJy/beam)', color='k')
ax2.set_ylim(ylimits)
# ax2.set_xlabel does not work in this case, so use plt.text instead
plt.text(0.5, 1.055, label, ha='center', va='center', transform=plt.gca().transAxes, size=11)
else:
ax2 = ax1.twiny()
ax2.set_xlabel(label, size=fontsize)
if xlimits is not None:
ax2.set_xlim(np.array(xlimits)*1e-9)
else:
ax2.set_xlim(firstFreq*1e-9, lastFreq*1e-9)
freqRange = np.abs(lastFreq-firstFreq)
power = int(np.log10(freqRange))-9
ax2.xaxis.set_major_locator(matplotlib.ticker.MultipleLocator(10**power))
if (len(ax2.get_xticks()) < 2):
ax2.xaxis.set_major_locator(matplotlib.ticker.MultipleLocator(0.5*10**power))
ax2.xaxis.set_minor_locator(matplotlib.ticker.MultipleLocator(0.1*10**power))
ax2.xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter(useOffset=False))
[docs]def numberOfChannelsInCube(img, returnFreqs=False, returnChannelWidth=False, verbose=False):
"""
Finds the number of channels in a CASA image cube, using the ia tool.
returnFreqs: if True, then also return the frequency of the
first and last channel (in Hz)
returnChannelWidth: if True, then also return the channel width (in Hz)
verbose: if True, then print the frequencies of first and last channel
Returns: an int
-Todd Hunter
"""
with casa_tools.ImageReader(img) as image:
lcs = image.coordsys()
try:
axis = lcs.findaxisbyname('spectral')
except:
if image.shape().shape[0] > 3:
LOG.warn("Can't find spectral axis. Assuming it is 3.")
axis = 3
elif image.shape().shape[0] > 2:
LOG.warn("Can't find spectral axis. Assuming it is 2.")
axis = 2
elif image.shape().shape[0] == 2:
LOG.error("No spectral axis found.")
raise Exception('No spectral axis found/')
naxes = image.shape().shape[0]
nchan = image.shape()[axis]
cdelt = lcs.increment()['numeric'][axis]
pixel = [0]*naxes
firstFreq = lcs.toworld(pixel, format='n')['numeric'][axis]
pixel[axis] = nchan-1
lastFreq = lcs.toworld(pixel, format='n')['numeric'][axis]
lcs.done()
nchan = int(nchan)
if returnFreqs:
if returnChannelWidth:
return nchan, firstFreq, lastFreq, cdelt
else:
return nchan, firstFreq, lastFreq
else:
if returnChannelWidth:
return nchan, cdelt
else:
return nchan
# FIXME: function contains unresolved references, clean-up or remove if not needed by Pipeline.
[docs]def cubeLSRKToTopo(img, freqrange='', prec=4, verbose=False,
nchan=None, f0=None, f1=None, chanwidth=None,
vis='', header=''):
"""
Reads the date of observation, central RA and Dec,
and observatory from an image cube header and then calls lsrkToTopo to
return the specified frequency range in TOPO (in Hz).
freqrange: desired range of frequencies (empty string or list = whole cube)
floating point list of two frequencies, or a delimited string
(delimiter = ',', '~' or space)
prec: in fractions of Hz (only used to display the value when verbose=True)
vis: read date of observation from the specified measurement set
header: dictionary output by imheadlist
-Todd Hunter
"""
if header == '':
header = imheadlist(img, omitBeam=True)
if nchan is None or f0 is None or f1 is None or chanwidth is None:
nchan, f0, f1, chanwidth = numberOfChannelsInCube(img, returnFreqs=True, returnChannelWidth=True)
if 'reffreqtype' in header:
if header['reffreqtype'].upper() == 'TOPO':
return np.array([f0, f1])
if len(freqrange) == 0:
startFreq = f0
stopFreq = f1
elif isinstance(freqrange, str):
if freqrange.find(',') > 0:
freqrange = [parseFrequencyArgument(i) for i in freqrange.split(',')]
elif freqrange.find('~') > 0:
freqrange = [parseFrequencyArgument(i) for i in freqrange.split('~')]
else:
freqrange = [parseFrequencyArgument(i) for i in freqrange.split()]
startFreq, stopFreq = freqrange
else:
startFreq, stopFreq = freqrange
ra, dec = rad2radec(header['crval1'], header['crval2'], delimiter=' ', verbose=False).split()
equinox = header['equinox']
observatory = header['telescope']
if vis == '':
datestring = header['date-obs']
else:
# datestring = getObservationStartDate(vis, measuresToolFormat=True)
# get the datetime stamp of the middle of the measurement set
# FIXME: fix unresolved reference and/or clean-up/remove function if not needed by Pipeline.
datestring = mjdsecToUT(np.mean(getObservationMJDSecRange(vis)), measuresToolFormat=True)
f0 = lsrkToTopo(startFreq, datestring, ra, dec, equinox, observatory, prec, verbose)
f1 = lsrkToTopo(stopFreq, datestring, ra, dec, equinox, observatory, prec, verbose)
return np.array([f0, f1])
[docs]def lsrkToTopo(lsrkFrequency, datestring, ra, dec, equinox='J2000', observatory='ALMA', prec=4, verbose=False):
"""
Converts an LSRKfrequency and observing date/direction
to the corresponding frequency in the TOPO frame.
Inputs:
lsrkFrequency: floating point value in Hz or GHz, or a string with units
datestring: "YYYY/MM/DD/HH:MM:SS" (format = image header keyword 'date-obs')
ra: string "HH:MM:SS.SSSS"
dec: string "DD.MM.SS.SSSS" or "DD:MM:SS.SSSS" (colons will be replaced with .)
prec: only used to display the value when verbose=True
Returns: the TOPO frequency in Hz
-Todd Hunter
"""
velocityLSRK = 0 # does not matter what it is, just needs to be same in both calls
restFreqHz = lsrkToRest(lsrkFrequency, velocityLSRK, datestring, ra, dec, equinox,
observatory, prec, verbose)
topoFrequencyHz = restToTopo(restFreqHz, velocityLSRK, datestring, ra, dec, equinox, observatory, verbose=verbose)
return topoFrequencyHz
[docs]def restToTopo(restFrequency, velocityLSRK, datestring, ra, dec, equinox='J2000',
observatory='ALMA', veltype='radio', verbose=False):
"""
Converts a rest frequency, LSRK velocity, and observing date/direction
to the corresponding frequency in the TOPO frame.
Inputs:
restFrequency: floating point value in Hz or GHz, or a string with units
velocityLSRK: floating point value in km/s
datestring: "YYYY/MM/DD/HH:MM:SS"
ra: string "HH:MM:SS.SSSS"
dec: string "DD.MM.SS.SSSS" or "DD:MM:SS.SSSS" (colons will be replaced with .)
prec: only used to display the value when verbose=True
Returns: the TOPO frequency in Hz
-Todd Hunter
"""
topoFreqHz, diff1, diff2 = frames(velocityLSRK, datestring, ra, dec, equinox,
observatory, verbose=verbose,
restFreq=restFrequency, veltype=veltype)
return topoFreqHz
[docs]def frames(velocity=286.7, datestring="2005/11/01/00:00:00",
ra="05:35:28.105", dec="-069.16.10.99", equinox="J2000",
observatory="ALMA", prec=4, verbose=True,
restFreq=345.79599, veltype='optical'):
"""
Converts an optical velocity into barycentric, LSRK and TOPO frames.
Converts a radio LSRK velocity into TOPO frame.
Inputs:
velocity: in km/s
datestring: "YYYY/MM/DD/HH:MM:SS"
ra: "05:35:28.105"
dec: "-069.16.10.99"
equinox: "J2000"
observatory: "ALMA"
prec: precision to display (digits to the right of the decimal point)
veltype: 'radio' or 'optical'
restFreq: in Hz, GHz or a string with units
Returns:
* TOPO frequency in Hz
* difference between LSRK-TOPO in km/sec
* difference between LSRK-TOPO in Hz
- Todd Hunter
"""
lme = casa_tools.measures
lqa = casa_tools.quanta
if dec.find(':') >= 0:
dec = dec.replace(':', '.')
position = lme.direction(equinox, ra, dec)
obstime = lme.epoch('TAI', datestring)
if veltype.lower().find('opt') == 0:
velOpt = lqa.quantity(velocity, "km/s")
dopp = lme.doppler("OPTICAL", velOpt)
# CASA doesn't do Helio, but difference to Bary is hopefully small
rvelOpt = lme.toradialvelocity("BARY", dopp)
elif veltype.lower().find('rad') == 0:
rvelOpt = lme.radialvelocity('LSRK', str(velocity)+'km/s')
else:
print("veltype must be 'rad'io or 'opt'ical")
return
lme.doframe(position)
lme.doframe(lme.observatory(observatory))
lme.doframe(obstime)
lme.showframe()
rvelRad = lme.measure(rvelOpt, 'LSRK')
doppRad = lme.todoppler("RADIO", rvelRad)
restFreq = parseFrequencyArgumentToGHz(restFreq)
freqRad = lme.tofrequency('LSRK', doppRad, lme.frequency('rest', str(restFreq)+'GHz'))
lsrk = lqa.tos(rvelRad['m0'], prec=prec)
rvelTop = lme.measure(rvelOpt, 'TOPO')
doppTop = lme.todoppler("RADIO", rvelTop)
freqTop = lme.tofrequency('TOPO', doppTop, lme.frequency('rest', str(restFreq)+'GHz'))
topo = lqa.tos(rvelTop['m0'], prec=prec)
velocityDifference = 0.001*(rvelRad['m0']['value']-rvelTop['m0']['value'])
frequencyDifference = freqRad['m0']['value'] - freqTop['m0']['value']
lme.done()
lqa.done()
return freqTop['m0']['value'], velocityDifference, frequencyDifference
[docs]def RescaleTrans(trans, lim):
# Input: the array of transmission values and current y-axis limits
# Returns: arrays of the rescaled transmission values and the zero point
# values in units of the frame, and in amplitude.
debug = False
yrange = lim[1]-lim[0]
labelgap = 0.5 # Use this fraction of the margin to separate the top
# curve from the upper y-axis
TOP_MARGIN = 0.25
y2 = lim[1] - labelgap*yrange*TOP_MARGIN/(1.0+TOP_MARGIN)
y1 = lim[1] - yrange*TOP_MARGIN/(1.0+TOP_MARGIN)
transmissionRange = np.max(trans)-np.min(trans)
if (transmissionRange < 0.05):
# force there to be a minimum range of transmission display
# overemphasize tiny ozone lines
transmissionRange = 0.05
# convert transmission to amplitude
newtrans = y2 - (y2-y1)*(np.max(trans)-trans)/transmissionRange
# Use edge values
edgeValueTransmission = trans[-1]
otherEdgeValueTransmission = trans[0]
# Now convert the edge channels' transmission values into amplitude
edgeValueAmplitude = y2 - (y2-y1)*(np.max(trans)-trans[-1])/transmissionRange
otherEdgeValueAmplitude = y2 - (y2-y1)*(np.max(trans)-trans[0])/transmissionRange
# Now convert amplitude to frame units, offsetting downward by half
# the font size
fontoffset = 0.01
edgeValueFrame = (edgeValueAmplitude - lim[0])/yrange - fontoffset
otherEdgeValueFrame = (otherEdgeValueAmplitude - lim[0])/yrange - fontoffset
# scaleFactor is how large the plot is from the bottom x-axis
# up to the labelgap, in units of the transmissionRange
scaleFactor = (1+TOP_MARGIN*(1-labelgap)) / (TOP_MARGIN*(1-labelgap))
# compute the transmission at the bottom of the plot, and label it
y0transmission = np.max(trans) - transmissionRange*scaleFactor
y0transmissionFrame = 0
y0transmissionAmplitude = lim[0]
if (y0transmission <= 0):
# If the bottom of the plot is below zero transmission, then label
# the location of zero transmission instead.
y0transmissionAmplitude = y1-(y2-y1)*(np.min(trans)/transmissionRange)
y0transmissionFrame = (y0transmissionAmplitude-lim[0]) / (lim[1]-lim[0])
y0transmission = 0
return(newtrans, edgeValueFrame, y0transmission, y0transmissionFrame,
otherEdgeValueFrame, edgeValueTransmission,
otherEdgeValueTransmission, edgeValueAmplitude,
otherEdgeValueAmplitude, y0transmissionAmplitude)
[docs]def lsrkToRest(lsrkFrequency, velocityLSRK, datestring, ra, dec,
equinox='J2000', observatory='ALMA', prec=4, verbose=True):
"""
Converts an LSRK frequency, LSRK velocity, and observing date/direction
to the corresponding frequency in the rest frame.
Inputs:
lsrkFrequency: floating point value in Hz or GHz, or a string with units
velocityLSRK: floating point value in km/s
datestring: "YYYY/MM/DD/HH:MM:SS" (format = image header keyword 'date-obs')
ra: string "HH:MM:SS.SSSS"
dec: string "DD.MM.SS.SSSS" or "DD:MM:SS.SSSS" (colons will be replaced with .)
prec: only used to display the value when verbose=True
Returns: the Rest frequency in Hz
-Todd Hunter
"""
if dec.find(':') >= 0:
dec = dec.replace(':', '.')
if verbose:
print("Warning: replacing colons with decimals in the dec field.")
freqGHz = parseFrequencyArgumentToGHz(lsrkFrequency)
lqa = casa_tools.quanta
lme = casa_tools.measures
velocityRadio = lqa.quantity(velocityLSRK, "km/s")
position = lme.direction(equinox, ra, dec)
obstime = lme.epoch('TAI', datestring)
dopp = lme.doppler("RADIO", velocityRadio)
radialVelocityLSRK = lme.toradialvelocity("LSRK", dopp)
lme.doframe(position)
lme.doframe(lme.observatory(observatory))
lme.doframe(obstime)
rvelRad = lme.measure(radialVelocityLSRK, 'LSRK')
doppRad = lme.todoppler('RADIO', rvelRad)
freqRad = lme.torestfrequency(lme.frequency('LSRK', str(freqGHz)+'GHz'), dopp)
lqa.done()
lme.done()
return freqRad['m0']['value']
[docs]def parseFrequencyArgumentToGHz(bandwidth):
"""
Converts a frequency string into floating point in GHz, based on the units.
If the units are not present, then the value is assumed to be GHz if less
than 10000, otherwise it assumes Hz.
-Todd Hunter
"""
value = parseFrequencyArgument(bandwidth)
if (value > 10000 or str(bandwidth).lower().find('hz') >= 0):
value *= 1e-9
return(value)
[docs]def parseFrequencyArgument(bandwidth):
"""
Converts a string frequency (or dictionary) into floating point in Hz, based on
the units. If the units are not present, then the value is simply converted to float.
-Todd Hunter
"""
if (isinstance(bandwidth, dict)):
bandwidth = str(bandwidth['value']) + bandwidth['unit']
else:
bandwidth = str(bandwidth)
ghz = bandwidth.lower().find('ghz')
mhz = bandwidth.lower().find('mhz')
khz = bandwidth.lower().find('khz')
hz = bandwidth.lower().find('hz')
if (ghz>0):
bandwidth = 1e9*float(bandwidth[:ghz])
elif (mhz>0):
bandwidth = 1e6*float(bandwidth[:mhz])
elif (khz>0):
bandwidth = 1e3*float(bandwidth[:khz])
elif (hz>0):
bandwidth = float(bandwidth[:hz])
else:
bandwidth = float(bandwidth)
return(bandwidth)
# FIXME: function contains unresolved references, clean-up or remove if not needed by Pipeline.
[docs]def rad2radec(ra=0,dec=0,imfitdict=None, prec=5, verbose=True, component=0,
replaceDecDotsWithColons=True, hmsdms=False, delimiter=', ',
prependEquinox=False, hmdm=False):
"""
Convert a position in RA/Dec from radians to sexagesimal string which
is comma-delimited, e.g. '20:10:49.01, +057:17:44.806'.
The position can either be entered as scalars via the 'ra' and 'dec'
parameters, as a tuple via the 'ra' parameter, as an array of shape (2,1)
via the 'ra' parameter, or
as an imfit dictionary can be passed via the 'imfitdict' argument, and the
position of component 0 will be displayed in RA/Dec sexagesimal.
replaceDecDotsWithColons: replace dots with colons as the Declination d/m/s delimiter
hmsdms: produce output of format: '20h10m49.01s, +057d17m44.806s'
hmdm: produce output of format: '20h10m49.01, +057d17m44.806' (for simobserve)
delimiter: the character to use to delimit the RA and Dec strings output
prependEquinox: if True, insert "J2000" before coordinates (i.e. for clean or simobserve)
Todd Hunter
"""
if (isinstance(imfitdict, dict)):
comp = 'component%d' % (component)
ra = imfitdict['results'][comp]['shape']['direction']['m0']['value']
dec = imfitdict['results'][comp]['shape']['direction']['m1']['value']
if (isinstance(ra, tuple) or isinstance(ra, list) or isinstance(ra, np.ndarray)):
if (len(ra) == 2):
dec = ra[1] # must come first before ra is redefined
ra = ra[0]
else:
ra = ra[0]
dec = dec[0]
if np.shape(ra) == (2, 1):
dec = ra[1][0]
ra = ra[0][0]
lqa = casa_tools.quanta
myra = lqa.formxxx('%.12frad' % ra, format='hms', prec=prec+1)
mydec = lqa.formxxx('%.12frad' % dec, format='dms', prec=prec-1)
if replaceDecDotsWithColons:
mydec = mydec.replace('.', ':', 2)
if len(mydec.split(':')[0]) > 3:
mydec = mydec[0] + mydec[2:]
mystring = '%s, %s' % (myra, mydec)
lqa.done()
if (hmsdms):
# FIXME: fix unresolved reference and/or clean-up/remove function if not needed by Pipeline.
mystring = convertColonDelimitersToHMSDMS(mystring)
if (prependEquinox):
mystring = "J2000 " + mystring
elif (hmdm):
# FIXME: fix unresolved reference and/or clean-up/remove function if not needed by Pipeline.
mystring = convertColonDelimitersToHMSDMS(mystring, s=False)
if (prependEquinox):
mystring = "J2000 " + mystring
if (delimiter != ', '):
mystring = mystring.replace(', ', delimiter)
return(mystring)
[docs]def imheadlist(vis, omitBeam=False):
"""
Emulates imhead(mode='list') but leaves off the min/max/minpos/maxpos
keywords, the filling of which makes it take so long to run on large cubes.
-Todd Hunter
"""
if (not os.path.exists(vis)):
print("Could not find image.")
return
header = {}
keys = ['bunit', 'date-obs', 'equinox', 'imtype', 'masks',
'object', 'observer', 'projection', 'reffreqtype',
'restfreq', 'shape', 'telescope']
if not omitBeam:
singleBeam = imhead(vis, mode='get', hdkey='beammajor')
if (singleBeam == False):
header = imhead(vis, mode='list')
if (header is None):
print("No beam found. Re-run with omitBeam=True.")
return -1
if 'perplanebeams' not in header:
print("No beam found. Re-run with omitBeam=True.")
return -1
beammajor = []
beamminor = []
beampa = []
for beamchan in range(header['perplanebeams']['nChannels']):
beamdict = header['perplanebeams']['*'+str(beamchan)]
beammajor.append(beamdict['major']['value'])
beamminor.append(beamdict['minor']['value'])
beampa.append(beamdict['positionangle']['value'])
bmaj = np.median(beammajor)
bmin = np.median(beamminor)
sinbpa = np.sin(np.radians(np.array(beampa)))
cosbpa = np.cos(np.radians(np.array(beampa)))
bpa = degrees(np.median(np.arctan2(np.median(sinbpa), np.median(cosbpa))))
header['beammajor'] = bmaj
header['beamminor'] = bmin
header['beampa'] = bpa
else:
keys += ['beammajor', 'beamminor', 'beampa']
for key in keys:
try:
header[key] = imhead(vis, mode='get', hdkey=key)
except:
pass
for axis in range(len(header['shape'])):
for key in ['cdelt', 'crval']:
mykey = key+str(axis+1)
try:
result = imhead(vis, mode='get', hdkey=mykey)
if (isinstance(result, dict)):
header[mykey] = result['value']
else:
# crval3 (pol axis) will be an array (e.g. ['I']) not a dict
header[mykey] = result
except:
print("Failed to set header key: ", mykey)
pass
for key in ['crpix', 'ctype', 'cunit']:
mykey = key+str(axis+1)
try:
header[mykey] = imhead(vis, mode='get', hdkey=mykey)
except:
pass
return(header)
[docs]def CalcAtmTransmissionForImage(img, chanInfo='', airmass=1.5, pwv=-1,
spectralaxis=-1, value='transmission', P=-1, H=-1,
T=-1, altitude=-1):
"""
This function is called by atmosphereVariation.
Supported telescopes are VLA and ALMA (needed for default weather and PWV)
img: name of CASA image
value: 'transmission' or 'tsky'
chanInfo: a list containing nchan, firstFreqHz, lastFreqHz, channelWidthHz
pwv: in mm
P: in mbar
H: in percent
T: in Kelvin
Returns:
2 arrays: frequencies (in GHz) and values (Kelvin, or transmission: 0..1)
"""
with casa_tools.ImageReader(img) as image:
lcs = image.coordsys()
telescopeName = lcs.telescope()
lcs.done()
if chanInfo == '':
chanInfo = numberOfChannelsInCube(img, returnChannelWidth=True, returnFreqs=True)
freqs = np.linspace(chanInfo[1]*1e-9, chanInfo[2]*1e-9, chanInfo[0])
numchan = len(freqs)
lsrkwidth = (chanInfo[2] - chanInfo[1])/(numchan-1)
result = cubeLSRKToTopo(img, chanInfo[1:3])
if (result is None):
topofreqs = freqs
else:
topoWidth = (result[1]-result[0])/(numchan-1)
topofreqs = np.linspace(result[0], result[1], chanInfo[0]) * 1e-9
#print("Converted LSRK range,width (%f-%f,%f) to TOPO (%f-%f,%f) over %d channels" % (chanInfo[1]*1e-9, chanInfo[2]*1e-9,lsrkwidth,topofreqs[0],topofreqs[-1],topoWidth,numchan))
P0 = 1000.0 # mbar
H0 = 20.0 # percent
T0 = 273.0 # Kelvin
if (telescopeName.find('ALMA') >= 0 or telescopeName.find('ACA') >= 0):
pwv0 = 1.0
P0 = 563.0
H0 = 20.0
T0 = 273.0
altitude0 = 5059
elif (telescopeName.find('VLA') >= 0):
P0 = 786.0
pwv0 = 5.0
altitude0 = 2124
else:
pwv0 = 10.0
altitude0 = 0
if (pwv < 0):
pwv = pwv0
if (T < 0):
T = T0
if (H < 0):
H = H0
if (P < 0):
P = P0
if (altitude < 0):
altitude = altitude0
tropical = 1
midLatitudeSummer = 2
midLatitudeWinter = 3
reffreq = 0.5*(topofreqs[numchan//2-1]+topofreqs[numchan//2])
# reffreq = np.mean(topofreqs)
numchanModel = numchan*1
chansepModel = (topofreqs[-1]-topofreqs[0])/(numchanModel-1)
nbands = 1
lqa = casa_tools.quanta
fCenter = lqa.quantity(reffreq, 'GHz')
fResolution = lqa.quantity(chansepModel, 'GHz')
fWidth = lqa.quantity(numchanModel*chansepModel, 'GHz')
myat = casa_tools.atmosphere
myat.initAtmProfile(humidity=H, temperature=lqa.quantity(T, "K"),
altitude=lqa.quantity(altitude, "m"),
pressure=lqa.quantity(P, 'mbar'), atmType=midLatitudeWinter)
myat.initSpectralWindow(nbands, fCenter, fWidth, fResolution)
myat.setUserWH2O(lqa.quantity(pwv, 'mm'))
# myat.setAirMass() # This does not affect the opacity, but it does effect TebbSky, so do it manually.
lqa.done()
dry = np.array(myat.getDryOpacitySpec(0)[1])
wet = np.array(myat.getWetOpacitySpec(0)[1]['value'])
TebbSky = myat.getTebbSkySpec(spwid=0)[1]['value']
# readback the values to be sure they got set
if (myat.getRefFreq()['unit'] != 'GHz'):
LOG.warn("There is a unit mismatch for refFreq in the atm code.")
if (myat.getChanSep()['unit'] != 'MHz'):
LOG.warn("There is a unit mismatch for chanSep in the atm code.")
numchanModel = myat.getNumChan()
freq0 = myat.getChanFreq(0)['value']
freq1 = myat.getChanFreq(numchanModel-1)['value']
# We keep the original LSRK freqs for overlay on the LSRK spectrum, but associate
# the transmission values from the equivalent TOPO freqs
newfreqs = np.linspace(freqs[0], freqs[-1], numchanModel) # fix for SCOPS-4815
transmission = np.exp(-airmass*(wet+dry))
TebbSky *= (1-np.exp(-airmass*(wet+dry)))/(1-np.exp(-wet-dry))
if value=='transmission':
values = transmission
else:
values = TebbSky
del myat
return(newfreqs, values)
[docs]def plot_spectra(image_robust_rms_and_spectra, rec_info, plotfile):
"""
Takes a pipeline-produced cube and plots the spectrum within the clean
mask (pixels with value=1 in the mask), and the noise spectrum from outside
the mask and within the 0.2-0.3 level (auto adjusted upward if image size
has been mitigated).
image_robust_rms_and_spectra: dictionary of spectra and metadata
rec_info: dictionary of receiver information (type (DSB/TSB), LO1 frequency)
"""
qaTool = casa_tools.quanta
cube = os.path.basename(image_robust_rms_and_spectra['nonpbcor_imagename'])
# Get spectral frame
with casa_tools.ImageReader(cube) as image:
lcs = image.coordsys()
frame = lcs.referencecode('spectral')[0]
lcs.done()
unmaskedPixels = image_robust_rms_and_spectra['nonpbcor_image_cleanmask_npoints']
if unmaskedPixels is None:
unmaskedPixels = 0
plt.clf()
desc = plt.subplot(111)
units = 'mJy'
fontsize = 9
# x axes
nchan = len(image_robust_rms_and_spectra['nonpbcor_image_cleanmask_spectrum'])
channels = np.arange(1, nchan + 1)
freq_ch1 = qaTool.getvalue(qaTool.convert(image_robust_rms_and_spectra['nonpbcor_image_non_cleanmask_freq_ch1'], 'Hz'))
freq_chN = qaTool.getvalue(qaTool.convert(image_robust_rms_and_spectra['nonpbcor_image_non_cleanmask_freq_chN'], 'Hz'))
freqs = np.linspace(freq_ch1, freq_chN, nchan)
# Flux density spectrum
intensity = image_robust_rms_and_spectra['nonpbcor_image_cleanmask_spectrum'] * 1000.
if unmaskedPixels > 0:
mycolor = 'r'
message = 'Red spectrum from %d pixels in flattened clean mask' % (unmaskedPixels)
plt.text(0.025, 0.96, message, transform=desc.transAxes, ha='left', fontsize=fontsize + 1, color='r')
else:
mycolor = 'b'
pblimit = image_robust_rms_and_spectra['nonpbcor_image_cleanmask_spectrum_pblimit']
plt.text(0.025, 0.96,
'Blue spectrum is mean from pixels above the pb=%.2f level (no clean mask was found)' % pblimit,
transform=desc.transAxes, ha='left', fontsize=fontsize+1, color='b')
# Noise spectrum
noise = image_robust_rms_and_spectra['nonpbcor_image_non_cleanmask_robust_rms'] * 1000.
plt.text(0.025, 0.93, 'Black spectrum is per-channel scaled MAD from imstat annulus and outside clean mask',
transform=desc.transAxes, ha='left', fontsize=fontsize+1)
# Turn off rightside y-axis ticks to make way for second y-axis
desc.yaxis.set_ticks_position('left')
plt.plot(channels, intensity, '-', color=mycolor)
plt.xlabel('%d Channels' % len(channels))
plt.ylabel('Flux density (%s)' % units, color=mycolor)
plt.tick_params(axis='y', labelcolor=mycolor, color=mycolor)
# Give a buffer between final data point and y-axes in order
# to be able to see high edge channel values
plt.xlim([channels[0] - len(channels) // 100, channels[-1] + len(channels) // 100])
freqDelta = (freqs[1]-freqs[0])*(len(channels)//100) # in Hz
freqLimits = np.array([freqs[0]-freqDelta, freqs[-1]+freqDelta]) # Hz
# Ignore edge channels (otherwise the first channel in ephemeris
# cubes may go way off scale
ylimits = np.array([np.min(intensity[1:-1]), np.max(intensity[1:-1])])
ylimits[1] += ylimits[1]-ylimits[0]
# Be sure the mean intensity spectrum is separated from lower y-axis by
# reducing the y lower limit by 5 percent.
yrange = ylimits[1]-ylimits[0]
ylimits[0] -= 0.05*yrange
ylimits[1] += 0.05*yrange
noiseLimits = [2*np.min(noise[1:-1])-np.max(noise[1:-1]),
np.max(noise[1:-1])]
yrange = noiseLimits[1]-noiseLimits[0]
noiseLimits[0] -= 0.1*yrange
noiseLimits[1] += 0.1*yrange
plt.ylim(ylimits)
# Plot horizontal dotted line at zero intensity
plt.plot(plt.xlim(), [0, 0], 'k:')
plt.text(0.5, 1.085, cube, transform=desc.transAxes, fontsize=fontsize, ha='center')
addFrequencyAxisAbove(plt.gca(), freqs[0], freqs[-1], frame,
twinx=True, ylimits=noiseLimits,
xlimits=freqLimits)
plt.plot(freqs * 1e-9, noise, 'k-')
# Plot continuum frequency ranges
fpattern = re.compile('([\d.]*)(~)([\d.]*)(\D*)')
cont_freq_ranges = fpattern.findall(image_robust_rms_and_spectra['cont_freq_ranges'].replace(';', ''))
for cont_freq_range in cont_freq_ranges:
fLowGHz = qaTool.getvalue(qaTool.convert(qaTool.quantity(float(cont_freq_range[0]), cont_freq_range[3]), 'GHz'))
fHighGHz = qaTool.getvalue(qaTool.convert(qaTool.quantity(float(cont_freq_range[2]), cont_freq_range[3]), 'GHz'))
fcLevel = plt.ylim()[0] + yrange * 0.025
plt.plot([fLowGHz, fHighGHz], [fcLevel] * 2, 'c-', lw=2)
# Overlay atmosphere transmission
freq, transmission = CalcAtmTransmissionForImage(cube)
rescaledY, edgeYvalue, zeroValue, zeroYValue, otherEdgeYvalue, edgeT, otherEdgeT, edgeValueAmplitude, otherEdgeValueAmplitude, zeroValueAmplitude = RescaleTrans(transmission, plt.ylim())
plt.plot(freq, rescaledY, 'm-')
if rec_info['type'] == 'DSB':
LO1 = float(qaTool.getvalue(qaTool.convert(rec_info['LO1'], 'GHz')))
imageFreq0 = 2.0 * LO1 - freq[0]
imageFreq1 = 2.0 * LO1 - freq[-1]
chanInfo = [len(freq), imageFreq0, imageFreq1, -(freqs[1]-freqs[0])]
imageFreq, imageTransmission = CalcAtmTransmissionForImage(cube, chanInfo)
results = RescaleTrans(imageTransmission, plt.ylim())
rescaledImage = results[0]
# You need to keep the signal sideband frequency range so that the overlay works!
plt.plot(freq, rescaledImage, 'm--')
plt.draw()
plt.savefig(plotfile)
plt.clf()
plt.close()