casa6
Public
Commits
Bob Garwood authored b2b17e3f2a8 Merge
| 23 23 | import matplotlib.transforms |
| 24 24 | import numpy as np |
| 25 25 | import pylab as pb |
| 26 26 | from casatasks import casalog |
| 27 27 | from casatools import (atmosphere, ctsys, measures, ms, msmetadata, quanta, |
| 28 28 | table) |
| 29 29 | from matplotlib.ticker import (FormatStrFormatter, MultipleLocator, |
| 30 30 | ScalarFormatter) |
| 31 31 | from six.moves import input, range |
| 32 32 | ​ |
| 33 + | # CAS-13722, CAS-13385 |
| 34 + | import warnings |
| 35 + | import matplotlib.cbook |
| 36 + | warnings.filterwarnings("ignore",category=matplotlib.cbook.MatplotlibDeprecationWarning) |
| 33 37 | ​ |
| 34 38 | PLOTBANDPASS_REVISION_STRING = "$Id: task_plotbandpass.py,v 1.102 2018/01/21 14:45:41 thunter Exp $" |
| 35 39 | TOP_MARGIN = 0.25 # Used if showatm=T or showtksy=T |
| 36 40 | BOTTOM_MARGIN = 0.25 # Used if showfdm=T |
| 37 41 | MAX_ATM_CALC_CHANNELS = 512 |
| 38 42 | ​ |
| 39 43 | markeredgewidth = 0.0 |
| 40 44 | ​ |
| 41 45 | # This is a color sequence found online which has distinguishable colors |
| 42 46 | overlayColors = [ |
| 715 719 | """ |
| 716 720 | This is a task to plot bandpass and Tsys calibration tables faster and more |
| 717 721 | flexibly than plotcal, including the ability to overlay the atmospheric |
| 718 722 | transmission, and to create multi-page plots as pngs and combine them |
| 719 723 | into single PDF documents. |
| 720 724 | It works with both the old style and new style cal tables. The source code |
| 721 725 | is in task_plotbandpass.py. For more detailed help, see examples at: |
| 722 726 | http://casaguides.nrao.edu/index.php?title=Plotbandpass |
| 723 727 | -- Todd Hunter |
| 724 728 | """ |
| 725 - | axes = dict() # keep track of already created axes |
| 726 729 | def safe_pb_subplot(xframe): |
| 727 730 | """ |
| 728 731 | CAS-12786: old pyplots (up to CASA 5.6.1 used to accept "220" in the pos parameter |
| 729 732 | Newer pyplots won't. Assume the index effectively used was 1 ("221") |
| 730 - | |
| 731 - | CAS-13276: pb.subplot will return a different instance in future matplotlib versions rather than |
| 732 - | the same instance. We need to keep track of the axes that have been already created; otherwise, we will lose plots. |
| 733 733 | """ |
| 734 - | if (axes.get(xframe) != None): |
| 735 - | return axes.get(xframe) |
| 736 - | ​ |
| 737 - | if str(xframe).endswith('0'): |
| 738 - | adesc = pb.subplot(xframe + 1) |
| 739 - | else: |
| 740 - | adesc = pb.subplot(xframe) |
| 741 - | |
| 742 - | axes[xframe] = adesc |
| 743 - | return adesc |
| 734 + | xf = (xframe + 1) if str(xframe).endswith('0') else xframe |
| 735 + | return pb.subplot(xf) |
| 736 + | |
| 737 + | def safe_pb_clf(): |
| 738 + | pb.clf() |
| 744 739 | ​ |
| 745 740 | casalog.origin('plotbandpass') |
| 746 741 | casalogPost(debug,"%s" % (PLOTBANDPASS_REVISION_STRING)) |
| 747 742 | DEBUG = debug |
| 748 743 | help = False |
| 749 744 | vm = '' # unused variable, now that msmd is available in casa |
| 750 745 | if (help): |
| 751 746 | print("Usage: plotbandpass(caltable='', antenna='', field='', spw='', yaxis='amp',") |
| 752 747 | print(" xaxis='chan', figfile='', plotrange=[0,0,0,0], caltable2='',") |
| 753 748 | print(" overlay='', showflagged=False, timeranges='', buildpdf=False, caltable3='',") |
| 2158 2153 | # I added pb.ion() because Remy suggested it. |
| 2159 2154 | if (interactive): |
| 2160 2155 | pb.ion() # This will open a new window if not present. |
| 2161 2156 | else: |
| 2162 2157 | pb.ioff() # This will not destroy an existing window or prevent new plots from appearing there. |
| 2163 2158 | # # The call to pb.figure() causes an additional new window everytime. |
| 2164 2159 | # # pb.figure() |
| 2165 2160 | |
| 2166 2161 | newylimits = [LARGE_POSITIVE, LARGE_NEGATIVE] |
| 2167 2162 | |
| 2168 - | pb.clf() |
| 2163 + | safe_pb_clf() # pb.clf() |
| 2169 2164 | if (bpoly): |
| 2170 2165 | # The number of polarizations cannot be reliably inferred from the shape of |
| 2171 2166 | # the GAIN column in the caltable. Must use the shape of the DATA column |
| 2172 2167 | # in the ms. |
| 2173 2168 | if (debug): print("in bpoly") |
| 2174 2169 | if (msFound): |
| 2175 2170 | (corr_type, corr_type_string, nPolarizations) = getCorrType(msName, spwsToPlot, mymsmd, debug) |
| 2176 2171 | casalogPost(debug,"nPolarizations in first spw to plot = %s" % (str(nPolarizations))) |
| 2177 2172 | else: |
| 2178 2173 | print("With no ms available, I will assume ALMA data: XX, YY, and refFreq=first channel.") |
| 2229 2224 | fieldString = str(field) |
| 2230 2225 | timeString = ', t%d/%d %s' % (mytime,nUniqueTimes-1,utstring(uniqueTimes[mytime],3)) |
| 2231 2226 | if (scansForUniqueTimes != []): |
| 2232 2227 | if (scansForUniqueTimes[mytime]>=0): |
| 2233 2228 | timeString = ', scan%d %s' % (scansForUniqueTimes[mytime],utstring(uniqueTimes[mytime],3)) |
| 2234 2229 | if ((yaxis.find('amp')>=0 or amplitudeWithPhase) and myap==0): |
| 2235 2230 | xframe += 1 |
| 2236 2231 | myUniqueColor = [] |
| 2237 2232 | if (debug): |
| 2238 2233 | print("v) incrementing xframe to %d" % xframe) |
| 2239 - | adesc = pb.subplot(xframe) |
| 2234 + | adesc = pb.safe_pb_subplot(xframe) |
| 2240 2235 | previousSubplot = xframe |
| 2241 2236 | if (ispw==originalSpw[ispw]): |
| 2242 2237 | # all this was added mistakenly here. If it causes a bug, remove it. |
| 2243 2238 | if (overlayTimes and len(fieldsToPlot) > 1): |
| 2244 2239 | indices = fstring = '' |
| 2245 2240 | for f in fieldIndicesToPlot: |
| 2246 2241 | if (f != fieldIndicesToPlot[0]): |
| 2247 2242 | indices += ',' |
| 2248 2243 | fstring += ',' |
| 2249 2244 | indices += str(uniqueFields[f]) |
| 2344 2339 | xant = antennasToPlot[xctr] |
| 2345 2340 | antstring, Antstring = buildAntString(xant,msFound,msAnt) |
| 2346 2341 | ispw = spwsToPlot[spwctr] |
| 2347 2342 | redisplay = True |
| 2348 2343 | else: |
| 2349 2344 | pagectr += 1 |
| 2350 2345 | if (pagectr >= len(pages)): |
| 2351 2346 | pages.append([xctr,spwctr,mytime,1]) |
| 2352 2347 | # print("appending [%d,%d,%d,%d]" % (xctr,spwctr,mytime,1)) |
| 2353 2348 | newpage = 0 |
| 2354 - | pb.clf() |
| 2349 + | safe_pb_clf() |
| 2355 2350 | |
| 2356 2351 | if (yaxis.find('phase')>=0 or amplitudeWithPhase): |
| 2357 2352 | xframe += 1 |
| 2358 2353 | myUniqueColor = [] |
| 2359 2354 | # # print("w) incrementing xframe to %d" % xframe) |
| 2360 - | adesc = pb.subplot(xframe) |
| 2355 + | adesc = pb.safe_pb_subplot(xframe) |
| 2361 2356 | previousSubplot = xframe |
| 2362 2357 | if (ispw==originalSpw[ispw]): |
| 2363 2358 | pb.title("%sspw%2d, field %d: %s%s" % (antennaString,ispw, |
| 2364 2359 | uniqueFields[fieldIndex],fieldString,timeString), size=titlesize) |
| 2365 2360 | else: |
| 2366 2361 | pb.title("%sspw%2d (%d), field %d: %s%s" % (antennaString,ispw,originalSpw[ispw], |
| 2367 2362 | uniqueFields[fieldIndex],fieldString,timeString), size=titlesize) |
| 2368 2363 | phaseSolutionX = calcChebyshev(polynomialPhase[index][0:nPolyPhase[index]], validDomain, frequenciesGHz[index]*1e+9) * 180/math.pi |
| 2369 2364 | phaseSolutionY = calcChebyshev(polynomialPhase[index][nPolyPhase[index]:2*nPolyPhase[index]], validDomain, frequenciesGHz[index]*1e+9) * 180/math.pi |
| 2370 2365 | if (nPolarizations == 1): |
| 2436 2431 | pagectr += 1 |
| 2437 2432 | if (pagectr >= len(pages)): |
| 2438 2433 | newpage = 1 |
| 2439 2434 | else: |
| 2440 2435 | newpage = 0 |
| 2441 2436 | if (debug): |
| 2442 2437 | print("1)Setting xframe to %d" % xframeStart) |
| 2443 2438 | xframe = xframeStart |
| 2444 2439 | if (xctr+1 < len(antennasToPlot)): |
| 2445 2440 | # don't clear the final plot when finished |
| 2446 - | pb.clf() |
| 2441 + | safe_pb_clf() |
| 2447 2442 | if (spwctr+1<len(spwsToPlot) or mytime+1<nUniqueTimes): |
| 2448 2443 | # don't clear the final plot when finished |
| 2449 - | pb.clf() |
| 2444 + | safe_pb_clf() |
| 2450 2445 | pb.subplots_adjust(hspace=myhspace, wspace=mywspace) |
| 2451 2446 | if (redisplay == False): |
| 2452 2447 | mytime += 1 |
| 2453 2448 | if (debug): |
| 2454 2449 | print("Incrementing mytime to %d" % (mytime)) |
| 2455 2450 | # end while(mytime) |
| 2456 2451 | if (redisplay == False): |
| 2457 2452 | spwctr +=1 |
| 2458 2453 | if (debug): |
| 2459 2454 | print("Incrementing spwctr to %d" % (spwctr)) |
| 3229 3224 | print("1)setting myUniqueTime to %d" % (mytime)) |
| 3230 3225 | myUniqueTime = mytime |
| 3231 3226 | ctr += 1 |
| 3232 3227 | if (ctr > len(fieldIndicesToPlot) and bOverlay==False): |
| 3233 3228 | print("multi-field time overlay *************** why are there 2 matches?") |
| 3234 3229 | # # # # if (ctr == 0): |
| 3235 3230 | # # # # print("No match for %.1f in "%(t), uTPFPS) |
| 3236 3231 | |
| 3237 3232 | # # # # print("Overlay antenna %d, myUniqueTime=%d" % (xctr, myUniqueTime)) |
| 3238 3233 | if (xframe == xframeStart): |
| 3239 - | pb.clf() |
| 3234 + | safe_pb_clf() |
| 3240 3235 | xflag = [item for sublist in xflag for item in sublist] |
| 3241 3236 | yflag = [item for sublist in yflag for item in sublist] |
| 3242 3237 | # # pflag = [xflag, yflag] |
| 3243 3238 | # # flagfrequencies = [frequencies, frequencies2] |
| 3244 3239 | antstring, Antstring = buildAntString(xant,msFound,msAnt) |
| 3245 3240 | if (msFound): |
| 3246 3241 | fieldString = msFields[uniqueFields[fieldIndex]] |
| 3247 3242 | else: |
| 3248 3243 | fieldString = str(field) |
| 3249 3244 | if (overlayTimes): |
| 4256 4251 | #redisplay the current page by setting ctrs back to the value they had at start of that page |
| 4257 4252 | xctr = pages[pagectr][PAGE_ANT] |
| 4258 4253 | spwctr = pages[pagectr][PAGE_SPW] |
| 4259 4254 | mytime = pages[pagectr][PAGE_TIME] |
| 4260 4255 | myap = pages[pagectr][PAGE_AP] |
| 4261 4256 | xant = antennasToPlot[xctr] |
| 4262 4257 | antstring, Antstring = buildAntString(xant,msFound,msAnt) |
| 4263 4258 | ispw = spwsToPlot[spwctr] |
| 4264 4259 | # # # # print("Returning to [%d,%d,%d,%d]" % (xctr,spwctr,mytime,myap)) |
| 4265 4260 | if (xctr==pages[0][PAGE_ANT] and spwctr==pages[0][PAGE_SPW] and mytime==pages[0][PAGE_TIME] and pages[0][PAGE_AP]==myap): |
| 4266 - | pb.clf() |
| 4261 + | safe_pb_clf() |
| 4267 4262 | if (debug): |
| 4268 4263 | print("2)Setting xframe to %d" % xframeStart) |
| 4269 4264 | xframe = xframeStart |
| 4270 4265 | myUniqueColor = [] |
| 4271 4266 | continue |
| 4272 4267 | else: |
| 4273 4268 | pagectr += 1 |
| 4274 4269 | if (pagectr >= len(pages)): |
| 4275 4270 | if (xframe == lastFrame and overlayTimes and overlayAntennas and xctr+1==len(antennasToPlot) and |
| 4276 4271 | yaxis=='amp'): |
| 4277 4272 | # I'm not sure why this works, but is needed to fix CAS-7154 |
| 4278 4273 | myspwctr = spwctr+1 |
| 4279 4274 | else: |
| 4280 4275 | myspwctr = spwctr |
| 4281 4276 | pages.append([xctr,myspwctr,mytime,1]) |
| 4282 4277 | if (debug): |
| 4283 4278 | print("amp: appending [%d,%d,%d,%d]" % (xctr,myspwctr,mytime,1)) |
| 4284 4279 | newpage = 0 |
| 4285 - | pb.clf() |
| 4280 + | safe_pb_clf() |
| 4286 4281 | if (debug): |
| 4287 4282 | print("3)Setting xframe to %d" % xframeStart) |
| 4288 4283 | xframe = xframeStart |
| 4289 4284 | myUniqueColor = [] |
| 4290 4285 | else: |
| 4291 4286 | if (debug): |
| 4292 4287 | print("::: Not done page: Not checking whether we need to set xframe=xframeStart") |
| 4293 4288 | print("::: xframe=%d ?= lastFrame=%d, amplitudeWithPhase=" % (xframe, lastFrame), amplitudeWithPhase) |
| 4294 4289 | print("::: xctr+1=%d ?= len(antennasToPlot)=%d" % (xctr+1,len(antennasToPlot))) |
| 4295 4290 | print(":::: mytimeTest = %s" % (mytimeTest)) |