Commits

Kana Sugimoto authored 6dad6fce682 
A fix by Todd in how channel frequencies in atmospheric model are obtained.

Frequencies of the edge channels are used instead of reference channel in the band.
No tags
a7e91aae486
master

gcwrap/python/scripts/task_plotbandpass.py

Modified
6 6 # bandpass solution tables with options to overlay them in various
7 7 # combinations, and/or with an atmospheric transmission or sky temperature
8 8 # model. It works with both the 'new' (casa 3.4) and 'old' calibration
9 9 # table formats, and allows for mixed mode spws (e.g. TDM and FDM for ALMA).
10 10 # It uses the new msmd tool to access the information about an ms.
11 11 #
12 12 # Todd R. Hunter February 2013
13 13 #
14 14 # To test: see plotbandpass_regression.py
15 15 #
16 -PLOTBANDPASS_REVISION_STRING = "$Id: task_plotbandpass.py,v 1.96 2017/05/19 15:28:20 thunter Exp $"
16 +PLOTBANDPASS_REVISION_STRING = "$Id: task_plotbandpass.py,v 1.100 2017/06/09 15:57:51 thunter Exp $"
17 17 import pylab as pb
18 18 import math, os, sys, re
19 19 import time as timeUtilities
20 20 import numpy as np
21 21 import re # used for testing if a string is a float
22 22 # import casadef # necessary to read the casa version strings, but only in pre-CASA 5.0
23 23 from taskinit import * # necessary for tb.open() to work
24 24 from matplotlib.ticker import MultipleLocator, FormatStrFormatter, ScalarFormatter
25 25 import matplotlib.transforms
26 26 import inspect
83 83 developerEmail = "thunter@nrao.edu"
84 84 ​
85 85 #class Polarization:
86 86 # taken from Stokes.h in casa, for reference only
87 87 # (Undefined, I,Q,U,V,RR,RL,LR,LL,XX,XY,YX,YY) = range(13)
88 88 ​
89 89 def version(showfile=True):
90 90 """
91 91 Returns the CVS revision number.
92 92 """
93 - myversion = "$Id: task_plotbandpass.py,v 1.96 2017/05/19 15:28:20 thunter Exp $"
93 + myversion = "$Id: task_plotbandpass.py,v 1.100 2017/06/09 15:57:51 thunter Exp $"
94 94 if (showfile):
95 95 print "Loaded from %s" % (__file__)
96 96 return myversion
97 97 ​
98 98 def refTypeToString(rtype):
99 99 rtypes = ['REST','LSRK','LSRD','BARY','GEO','TOPO','GALACTO','LGROUP','CMB']
100 100 return(rtypes[rtype])
101 101 ​
102 102 def corrTypeToString(ptype):
103 103 ptypes = ['Undefined','I','Q','U','V','RR','RL','LR','LL','XX','XY','YX','YY']
660 660 casalog.post(mystring)
661 661 if (debug): print mystring
662 662
663 663 def computeHighestSpwIndexInSpwsToPlotThatHasCurrentScan(spwsToPlot, scansToPlotPerSpw, scan):
664 664 highestSpwIndex = -1
665 665 for i,spw in enumerate(spwsToPlot):
666 666 if (scan in scansToPlotPerSpw[spw]):
667 667 highestSpwIndex = i
668 668 return(highestSpwIndex)
669 669 ​
670 -​
671 670 DEFAULT_PLATFORMING_THRESHOLD = 10.0 # unused if platformingSigma != 0
672 671 def plotbandpass(caltable='', antenna='', field='', spw='', yaxis='amp',
673 672 xaxis='chan', figfile='', plotrange=[0,0,0,0],
674 673 caltable2='', overlay='', showflagged=False, timeranges='',
675 674 buildpdf=False, caltable3='', markersize=3, density=108,
676 675 interactive=True, showpoints='auto', showlines='auto',
677 676 subplot='22', zoom='', poln='', showatm=False, pwv='auto',
678 677 gs='gs', convert='convert', chanrange='',
679 678 solutionTimeThresholdSeconds=30.0, debug=False,
680 679 phase='', vis='',showtsky=False, showfdm=False,showatmfield='',
3076 3075 if (type(uFFI) == list or type(uFFI) == type(np.ndarray(0))):
3077 3076 uFFI = uFFI[0]
3078 3077 if (debug):
3079 3078 print "uFFI = %s" % (str(uFFI))
3080 3079 (atmfreqImage,atmchanImage,transmissionImage,pwvmean,atmairmass,TebbSkyImage,missingCalWVRErrorPrinted) = \
3081 3080 CalcAtmTransmission(channels, frequenciesImage, xaxis,
3082 3081 pwv, vm, mymsmd, msName, asdm, xant, uniqueTimes[mytime],
3083 3082 interval, uFFI, refFreq[originalSpw[ispw]],
3084 3083 net_sideband[originalSpw[ispw]], mytime,
3085 3084 missingCalWVRErrorPrinted, caltable, verbose=DEBUG)
3086 - # difference between the requested Atm freq and actual Atm calcuation
3087 - chanDifference = atmfreq[0]-frequencies[0]
3088 - chanImageDifference = atmfreqImage[0]-frequenciesImage[-1]
3089 - atmfreqImage = list(2*LO1 - np.array(atmfreqImage) + 2*chanDifference + chanImageDifference) # CAS-7715 adds final 2 terms
3090 - atmfreqImage.reverse()
3085 + # difference between the requested Atm freq and actual Atm calcuation CAS-7715. No longer necessary after CAS-10228.
3086 + #chanDifference = atmfreq[0]-frequencies[0]
3087 + #chanImageDifference = atmfreqImage[0]-frequenciesImage[-1]
3088 + #print "signal SB difference = %f, image SB difference = %f" % (chanDifference, chanImageDifference)
3089 + atmfreqImage = list(2*LO1 - np.array(atmfreqImage)) # + 2*chanDifference + chanImageDifference) # CAS-7715 adds final 2 terms
3091 3090 atmchanImage.reverse()
3092 3091
3093 3092 if (overlayTimes):
3094 3093 atmString = 'PWV %.2fmm, airmass %.2f (field %d)' % (pwvmean,atmairmass,showatmfield)
3095 3094 else:
3096 3095 atmString = 'PWV %.2fmm, airmass %.3f' % (pwvmean,atmairmass)
3097 3096 if (bOverlay):
3098 3097 for i in range(nRows2):
3099 3098 if (overlayTimes or overlayAntennas or len(fieldsToPlot)>1 or
3100 3099 (nFields>1 and len(fieldlist)<nFields)):
5539 5538 if (conditions['elevation'] <= 3):
5540 5539 print "Using 45 deg elevation instead"
5541 5540 airmass = 1.0/math.cos(45*math.pi/180.)
5542 5541 else:
5543 5542 airmass = 1.0/math.cos((90-conditions['elevation'])*math.pi/180.)
5544 5543 ​
5545 5544 tropical = 1
5546 5545 midLatitudeSummer = 2
5547 5546 midLatitudeWinter = 3
5548 5547 numchan = len(freqs)
5548 + # Set the reference freq to be the middle of the middle two channels
5549 5549 reffreq=0.5*(freqs[numchan/2-1]+freqs[numchan/2])
5550 5550 originalnumchan = numchan
5551 5551 while (numchan > MAX_ATM_CALC_CHANNELS):
5552 5552 numchan /= 2
5553 5553 # print "Reducing numchan to ", numchan
5554 5554 chans = range(0,originalnumchan,(originalnumchan/numchan))
5555 5555 ​
5556 5556 chansep = (freqs[-1]-freqs[0])/(numchan-1)
5557 5557 nbands = 1
5558 5558 if (verbose): print "Opening casac.atmosphere()"
5582 5582 n = myat.getNumChan()
5583 5583 if (verbose): print "numchan = %s" % (str(n))
5584 5584 if cu.compare_version('<',[4,0,0]):
5585 5585 dry = np.array(myat.getDryOpacitySpec(0)['dryOpacity'])
5586 5586 wet = np.array(myat.getWetOpacitySpec(0)['wetOpacity'].value)
5587 5587 TebbSky = []
5588 5588 for chan in range(n): # do NOT use numchan here, use n
5589 5589 TebbSky.append(myat.getTebbSky(nc=chan, spwid=0).value)
5590 5590 TebbSky = np.array(TebbSky)
5591 5591 # readback the values to be sure they got set
5592 - rf = myat.getRefFreq().value
5593 - cs = myat.getChanSep().value
5592 + #rf = myat.getRefFreq().value
5593 + #cs = myat.getChanSep().value
5594 5594 else: # casa >=4.0
5595 5595 dry = np.array(myat.getDryOpacitySpec(0)[1])
5596 5596 wet = np.array(myat.getWetOpacitySpec(0)[1]['value'])
5597 5597 TebbSky = myat.getTebbSkySpec(spwid=0)[1]['value']
5598 5598 # readback the values to be sure they got set
5599 - rf = myat.getRefFreq()['value']
5600 - cs = myat.getChanSep()['value']
5601 - if (myat.getRefFreq()['unit'] != 'GHz'):
5602 - print "There is a unit mismatch for refFreq in the code."
5603 - if (myat.getChanSep()['unit'] != 'MHz'):
5604 - print "There is a unit mismatch for chanSep in the code."
5599 + #rf = myqa.convert(myat.getRefFreq(),'GHz')['value']
5600 + #cs = myqa.convert(myat.getChanSep(),'GHz')['value']
5605 5601 ​
5606 5602 transmission = np.exp(-airmass*(wet+dry))
5607 5603 TebbSky *= (1-np.exp(-airmass*(wet+dry)))/(1-np.exp(-wet-dry))
5608 5604 ​
5609 5605 if (refFreqInTable*1e-9>np.mean(freqs)):
5610 5606 if ((net_sideband % 2) == 0):
5611 5607 sense = 1
5612 5608 else:
5613 5609 sense = 2
5614 5610 else:
5615 5611 if ((net_sideband % 2) == 0):
5616 5612 sense = 2
5617 5613 else:
5618 5614 sense = 1
5619 5615 ​
5620 5616 if (sense == 1):
5621 5617 # The following looks right for LSB sense=1
5622 -# freq = rf.value + cs.value*0.001*(0.5*n-1-np.array(range(n)))
5623 -# freq = rf + cs*0.001*(0.5*n-1-np.array(range(n)))
5624 5618 if (xaxis.find('chan')>=0):
5625 5619 trans = np.zeros(len(transmission))
5626 5620 Tebb = np.zeros(len(TebbSky))
5627 5621 for i in range(len(transmission)):
5628 5622 trans[i] = transmission[len(transmission)-1-i]
5629 5623 Tebb[i] = TebbSky[len(TebbSky)-1-i]
5630 5624 transmission = trans
5631 5625 TebbSky = Tebb
5632 -# else:
5633 - # Using numchan can cause an inconsistency for small number of channels
5634 -# freq = rf.value+cs.value*0.001*(np.array(range(numchan))-0.5*numchan+1)
5635 - # The following looks right for USB sense=2
5636 -# freq = rf+cs*0.001*(np.array(range(n))-0.5*n+1)
5637 5626
5638 5627 # Be sure that number of frequencies matched number of transmission values - CAS-10123
5639 5628 numchan = len(transmission)
5640 - chansep = cs*0.001
5641 5629 chans = range(len(transmission))
5642 - if sense == 2:
5643 - freq = np.linspace(reffreq-numchan*chansep/2, reffreq+numchan*chansep/2, numchan)
5644 - else:
5645 - freq = np.linspace(reffreq+numchan*chansep/2, reffreq-numchan*chansep/2, numchan) + chansep
5646 - if cu.compare_version('>=',[5,0,0]):
5647 - freq += chansep*0.5
5630 + # Note that getChanFreq returns units of GHz, but use convert to be sure.
5631 + startFreq = myqa.convert(myat.getChanFreq(0),'GHz')['value']
5632 + endFreq = myqa.convert(myat.getChanFreq(numchan-1),'GHz')['value']
5633 + # print "startFreq=%f endFreq=%f " % (startFreq, endFreq)
5634 + freq = np.linspace(startFreq, endFreq, numchan)
5635 +# old method that fails on spws with an even number of channels, i.e. when the integer refchan
5636 +# is half a channel from the center of the span
5637 +# if sense == 2:
5638 +# freq = np.linspace(rf-((numchan-1)/2.)*chansepGHz, rf+((numchan-1)/2.)*chansepGHz, numchan)
5639 +# else:
5640 +# freq = np.linspace(rf+((numchan-1)/2.)*chansepGHz,
5641 +# rf-((numchan-1)/2.)*chansepGHz, numchan)
5642 + # Therewas a 1-channel offset in CASA 5.0.x (CAS-10228), but it was fixed.
5643 +# if (cu.compare_version('<',[5,1,0])):
5644 +# freq += chansepGHz
5645 +​
5648 5646 if (verbose): print "Done CalcAtmTransmission"
5649 5647 return(freq, chans, transmission, pwvmean, airmass, TebbSky, missingCalWVRErrorPrinted)
5650 5648 ​
5651 5649 def RescaleTrans(trans, lim, subplotRows, lo1='', xframe=0):
5652 5650 # Input: the array of transmission or TebbSky values and current limits
5653 5651 # Returns: arrays of the rescaled transmission values and the zero point
5654 5652 # values in units of the frame, and in amplitude.
5655 5653 debug = False
5656 5654 yrange = lim[1]-lim[0]
5657 5655 if (lo1 == ''):
7255 7253 print myline
7256 7254 ​
7257 7255 # Loop over all rows in the ASDM_RECEIVER table, unless we've split, in
7258 7256 # which case this will loop over the N spws in the table.
7259 7257 needToClose = False
7260 7258 if mymsmd is None or mymsmd == '':
7261 7259 mymsmd = createCasaTool(msmdtool)
7262 7260 mymsmd.open(vis)
7263 7261 needToClose = True
7264 7262 if (spwsForIntent == None):
7265 - scienceSpws = np.setdiff1d(mymsmd.spwsforintent(intent),mymsmd.wvrspws())
7263 + if intent in mymsmd.intents(): # prevent a warning of OBSERVE_TARGET does not exist
7264 + scienceSpws = np.setdiff1d(mymsmd.spwsforintent(intent),mymsmd.wvrspws())
7265 + else:
7266 + scienceSpws = []
7266 7267 else:
7267 7268 scienceSpws = spwsForIntent
7268 7269 birdieIF = 0
7269 7270 if (birdieFreq is not None):
7270 7271 birdieFreq = parseFrequencyArgumentToHz(birdieFreq)
7271 7272 birdieFreqGHz = parseFrequencyArgumentToGHz(birdieFreq)
7272 7273 fdmSpws = mymsmd.almaspws(fdm=True)
7273 7274 for i in range(len(index)):
7274 7275 if (verbose):
7275 7276 print "index[%d]=%d" % (i,index[i])

Everything looks good. We'll let you know here if there's anything you should know about.

Add shortcut