Commits

Kana Sugimoto authored a9c362e230c Merge 
Merge pull request #101 in CASA/casa from bugfix/CAS-10228 to master

* commit '6dad6fce682f6b95194c107438e0abfa617a1825':
  A fix by Todd in how channel frequencies in atmospheric model are obtained.
  A workaround fix to 1-based refChan in getRefChan().
  Fixed initSpectralWindow and addSpectralWindow to workaround ICT-9490.
6541319b7b5
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