Commits
Kazuhiko Shimada authored and Ville Suoranta committed 22928258c5e Merge
18 18 | import math |
19 19 | import os |
20 20 | import re # used for testing if a string is a float |
21 21 | import time |
22 22 | import copy |
23 23 | |
24 24 | import matplotlib.transforms |
25 25 | import numpy as np |
26 26 | import pylab as pb |
27 27 | from casatasks import casalog |
28 + | from casatasks.private import sdutil, simutil |
28 29 | from casatools import (atmosphere, ctsys, measures, ms, msmetadata, quanta, |
29 30 | table) |
30 31 | from matplotlib.ticker import (FormatStrFormatter, MultipleLocator, |
31 32 | ScalarFormatter) |
32 33 | from six.moves import input, range |
33 34 | |
34 35 | # CAS-13722, CAS-13385 |
35 36 | import warnings |
36 37 | import matplotlib.cbook |
37 38 | #warnings.filterwarnings("ignore",category=matplotlib.cbook.MatplotlibDeprecationWarning) |
5663 5664 | else: |
5664 5665 | try: |
5665 5666 | pwvmean = float(pwv) |
5666 5667 | except: |
5667 5668 | pwvmean = defaultPWV |
5668 5669 | |
5669 5670 | if (verbose): |
5670 5671 | print("Using PWV = %.2f mm" % pwvmean) |
5671 5672 | |
5672 5673 | # default values in case we can't find them below |
5674 + | myqa = quanta() |
5673 5675 | airmass = 1.5 |
5674 5676 | P = 563.0 |
5675 5677 | H = 20.0 |
5676 5678 | T = 273.0 |
5677 5679 | roundedScanTimes = [] |
5678 5680 | if (type(mymsmd) != str): |
5679 5681 | if (verbose): |
5680 5682 | print("Looking for scans for field integer = %d" % (field)) |
5681 5683 | scans = mymsmd.scansforfield(field) |
5682 5684 | if (verbose): |
5698 5700 | # if (verbose): print("tdiff = %s" % (str(tdiff))) |
5699 5701 | if (tdiff < mindiff): |
5700 5702 | bestscan = scans[i] |
5701 5703 | if (verbose): print("bestscan = %s" % (str(bestscan))) |
5702 5704 | mindiff = tdiff |
5703 5705 | if (verbose): |
5704 5706 | print("For timestamp=%.1f, got closest scan = %d, %.0f sec away" %(timestamp, bestscan,mindiff)) |
5705 5707 | if (verbose): print("Calling getWeather()") |
5706 5708 | [conditions,myTimes] = getWeather(vis,bestscan,antenna,verbose,mymsmd) |
5707 5709 | if (verbose): print("Done getWeather()") |
5708 - | P = conditions['pressure'] |
5710 + | |
5711 + | # convert pressure with unit to the value in mbar |
5712 + | P = myqa.convert(myqa.quantity(conditions['pressure'], conditions['pressure_unit']), 'mbar')['value'] |
5713 + | |
5709 5714 | H = conditions['humidity'] |
5710 5715 | T = conditions['temperature']+273.15 |
5711 5716 | if (P <= 0.0): |
5712 5717 | P = 563 |
5713 5718 | if (H <= 0.0): |
5714 5719 | H = 20 |
5715 5720 | else: |
5716 5721 | conditions = {} |
5717 5722 | if (type(mymsmd) != str): |
5718 5723 | if ('elevation' not in list(conditions.keys())): |
5739 5744 | bestscan = -1 |
5740 5745 | if (verbose): |
5741 5746 | print("CalcAtm: found elevation=%f (airmass=%.3f) for scan: %s" % (conditions['elevation'],1/np.sin(conditions['elevation']*np.pi/180.), str(bestscan))) |
5742 5747 | print("P,H,T = %f,%f,%f" % (P,H,T)) |
5743 5748 | if (conditions['elevation'] <= 3): |
5744 5749 | print("Using 45 deg elevation instead") |
5745 5750 | airmass = 1.0/math.cos(45*math.pi/180.) |
5746 5751 | else: |
5747 5752 | airmass = 1.0/math.cos((90-conditions['elevation'])*math.pi/180.) |
5748 5753 | |
5754 + | # get the elevation of the antenna |
5755 + | geodetic_elevation = 5059 |
5756 + | if os.path.exists(os.path.join(vis, 'ANTENNA')): |
5757 + | with sdutil.table_manager(os.path.join(vis, 'ANTENNA')) as tb: |
5758 + | _X, _Y, _Z = (float(i) for i in tb.getcell('POSITION', antenna)) |
5759 + | geodetic_elevation = simutil.simutil().xyz2long(_X, _Y, _Z, 'WGS84')[2] |
5760 + | |
5749 5761 | tropical = 1 |
5750 5762 | midLatitudeSummer = 2 |
5751 5763 | midLatitudeWinter = 3 |
5752 5764 | numchan = len(freqs) |
5753 5765 | # Set the reference freq to be the middle of the middle two channels |
5754 5766 | reffreq=0.5*(freqs[numchan//2-1]+freqs[numchan//2]) |
5755 5767 | originalnumchan = numchan |
5756 5768 | while (numchan > MAX_ATM_CALC_CHANNELS): |
5757 5769 | numchan //= 2 |
5758 5770 | # print("Reducing numchan to ", numchan) |
5759 5771 | chans = list(range(0,originalnumchan,(originalnumchan//numchan))) |
5760 5772 | |
5761 5773 | chansep = (freqs[-1]-freqs[0])/(numchan-1) |
5762 5774 | nbands = 1 |
5763 5775 | if (verbose): print("Opening casac.atmosphere()") |
5764 5776 | myat = atmosphere() |
5765 5777 | if (verbose): print("Opened") |
5766 - | myqa = quanta() |
5767 5778 | fCenter = myqa.quantity(reffreq,'GHz') |
5768 5779 | fResolution = myqa.quantity(chansep,'GHz') |
5769 5780 | fWidth = myqa.quantity(numchan*chansep,'GHz') |
5770 - | myat.initAtmProfile(humidity=H,temperature=myqa.quantity(T,"K"),altitude=myqa.quantity(5059,"m"),pressure=myqa.quantity(P,'mbar'),atmType=midLatitudeWinter) |
5781 + | myat.initAtmProfile(humidity=H, temperature=myqa.quantity(T, "K"), |
5782 + | altitude=myqa.quantity(geodetic_elevation, "m"), |
5783 + | pressure=myqa.quantity(P, 'mbar'), atmType=midLatitudeWinter) |
5771 5784 | myat.initSpectralWindow(nbands,fCenter,fWidth,fResolution) |
5772 5785 | myat.setUserWH2O(myqa.quantity(pwvmean,'mm')) |
5773 5786 | |
5774 5787 | # myat.setAirMass() # This does not affect the opacity, but it does effect TebbSky, so do it manually. |
5775 5788 | |
5776 5789 | n = myat.getNumChan() |
5777 5790 | if (verbose): print("numchan = %s" % (str(n))) |
5778 5791 | if ctsys.compare_version('<',[4,0,0]): |
5779 5792 | dry = np.array(myat.getDryOpacitySpec(0)['dryOpacity']) |
5780 5793 | wet = np.array(myat.getWetOpacitySpec(0)['wetOpacity'].value) |
7186 7199 | mytb.open("%s/WEATHER" % vis) |
7187 7200 | except: |
7188 7201 | print("Could not open the WEATHER table for this ms.") |
7189 7202 | mytb.done() |
7190 7203 | return([conditions,myTimes]) |
7191 7204 | if (True): |
7192 7205 | mjdsec = mytb.getcol('TIME') |
7193 7206 | indices = np.argsort(mjdsec) |
7194 7207 | mjd = mjdsec/86400. |
7195 7208 | pressure = mytb.getcol('PRESSURE') |
7209 + | conditions['pressure_unit'] = mytb.getcolkeywords('PRESSURE').get('QuantumUnits', ['mbar'])[0] |
7196 7210 | relativeHumidity = mytb.getcol('REL_HUMIDITY') |
7197 7211 | temperature = mytb.getcol('TEMPERATURE') |
7198 7212 | if (np.median(temperature) > 100): |
7199 7213 | # must be in units of Kelvin, so convert to C |
7200 7214 | temperature -= 273.15 |
7201 7215 | if 'DEW_POINT' in mytb.colnames(): |
7202 7216 | dewPoint = mytb.getcol('DEW_POINT') |
7203 7217 | if (np.median(dewPoint) > 100): |
7204 7218 | # must be in units of Kelvin, so convert to C |
7205 7219 | dewPoint -= 273.15 |