casa / casa6 / 03ba04cb248

Commits

Jorge Lopez authored and Ville Suoranta committed 03ba04cb248 Merge 12 Sep 2022

Pull request #631: CAS-13910: use np.median instead of np.mean in getWeather

Merge in CASA/casa6 from CAS-13910 to master

* commit '613335d42381d1b8409d5c7381f496728ac9ccee':
  CAS-13910: use np.median instead of np.mean in getWeather

casatasks/src/private/task_plotbandpass.py

Modified

-6776
+    mysep = np.abs(mytimes[0]-mytime)
-6777
+    for m in range(1,len(mytimes)):
-6778
+        if (np.abs(mytimes[m] - mytime) < mysep):
-6779
+            mysep = np.abs(mytimes[m] - mytime)
-6780
+            myindex = m
-6781
+    return(myindex)
-6782
-6783
+def getWeather(vis='', scan='', antenna='0',verbose=False, mymsmd=None):
-6784
+    """
-6785
+    Queries the WEATHER and ANTENNA tables of an .ms by scan number or
--
+    list of scan numbers in order to return mean values of: angleToSun,
-+
+    list of scan numbers in order to return median values of: angleToSun,
-6787
+      pressure, temperature, humidity, dew point, wind speed, wind direction,
-6788
+      azimuth, elevation, solarangle, solarelev, solarazim.
-6789
+    If the sun is below the horizon, the solarangle returned is negated.
-6790
+    -- Todd Hunter
-6791
+    """
-6792
+    if (verbose):
-6793
+        print("Entered getWeather with vis,scan,antenna = %s,%s,%s" % (str(vis), str(scan), str(antenna)))
-6794
+    try:
-6795
+        if str(antenna).isdigit():
-6796
+            antennaName = mymsmd.antennanames(antenna)[0]
-6871
+        listfield = ""
-6872
+        for field in listfields:
-6873
+            listfield += "%s" % field
-6874
+            if (field != listfields[-1]):
-6875
+                listfield += ","
-6876
+    else:
-6877
+        listscan = str(scan)
-6878
+        listfield = mymsmd.fieldsforscan(scan)
-6879
+    [az,el] = ComputeSolarAzElForObservatory(myTimes[0], mymsmd)
-6880
+    [az2,el2] = ComputeSolarAzElForObservatory(myTimes[-1], mymsmd)
--
+    azsun = np.mean([az,az2])
--
+    elsun = np.mean([el,el2])
-+
+    azsun = np.median([az,az2])
-+
+    elsun = np.median([el,el2])
-6883
+    direction = subtable.getcol("DIRECTION")
-6884
+    azeltime = subtable.getcol("TIME")
-6885
+    subtable.close()
-6886
+    telescopeName = mymsmd.observatorynames()[0]
-6887
+    if (len(direction) > 0 and telescopeName.find('VLA') < 0 and telescopeName.find('NRO') < 0):
-6888
+      azimuth = direction[0][0]*180.0/math.pi  # a list of values
-6889
+      elevation = direction[1][0]*180.0/math.pi # a list of values
-6890
+      npat = np.array(azeltime)
-6891
+      matches = np.where(npat>myTimes[0])[0]
-6892
+      matches2 = np.where(npat<myTimes[-1])[0]
-6893
+      if (len(matches2) > 0 and len(matches) > 0):
-6894
+          if verbose: print("matches[0]=%d, matches2[-1]=%d" % (matches[0],matches[-1]))
-6895
+          matchingIndices = list(range(matches[0],matches2[-1]+1))
-6896
+      else:
-6897
+          matchingIndices = []
-6898
+      if (len(matchingIndices) > 0):  # CAS-8440
--
+          conditions['azimuth'] = np.mean(azimuth[matches[0]:matches2[-1]+1])
--
+          conditions['elevation'] = np.mean(elevation[matches[0]:matches2[-1]+1])
-+
+          conditions['azimuth'] = np.median(azimuth[matches[0]:matches2[-1]+1])
-+
+          conditions['elevation'] = np.median(elevation[matches[0]:matches2[-1]+1])
-6901
+      elif (len(matches) > 0):        # CAS-8440
--
+          if verbose: print("using mean of all az/el values after time 0")
--
+          conditions['azimuth'] = np.mean(azimuth[matches[0]])
--
+          conditions['elevation'] = np.mean(elevation[matches[0]])
-+
+          if verbose: print("using median of all az/el values after time 0")
-+
+          conditions['azimuth'] = np.median(azimuth[matches[0]])
-+
+          conditions['elevation'] = np.median(elevation[matches[0]])
-6905
+      else:                           # CAS-8440
--
+          if verbose: print("using mean of all az/el values")
--
+          conditions['azimuth'] = np.mean(azimuth)
--
+          conditions['elevation'] = np.mean(elevation)
-+
+          if verbose: print("using median of all az/el values")
-+
+          conditions['azimuth'] = np.median(azimuth)
-+
+          conditions['elevation'] = np.median(elevation)
-6909
+      conditions['solarangle'] = angularSeparation(azsun,elsun,conditions['azimuth'],conditions['elevation'])
-6910
+      conditions['solarelev'] = elsun
-6911
+      conditions['solarazim'] = azsun
-6912
+      if (verbose):
--
+          print("Using antenna = %s to retrieve mean azimuth and elevation" % (antennaName))
-+
+          print("Using antenna = %s to retrieve median azimuth and elevation" % (antennaName))
-6914
+          print("Separation from sun = %f deg" % (abs(conditions['solarangle'])))
-6915
+      if (elsun<0):
-6916
+        conditions['solarangle'] = -conditions['solarangle']
-6917
+        if (verbose):
-6918
+            print("Sun is below horizon (elev=%.1f deg)" % (elsun))
-6919
+      else:
-6920
+        if (verbose):
-6921
+            print("Sun is above horizon (elev=%.1f deg)" % (elsun))
-6922
+      if (verbose):
-6923
+          print("Average azimuth = %.2f, elevation = %.2f degrees" % (conditions['azimuth'],conditions['elevation']))
-6926
+      if (type(scan) == int or type(scan)==np.int32):
-6927
+          # compute Az/El for this scan
-6928
+        myfieldId = mymsmd.fieldsforscan(scan)
-6929
+        if (type(myfieldId) == list or type(myfieldId) == type(np.ndarray(0))):
-6930
+            myfieldId = myfieldId[0]
-6931
+        fieldName = mymsmd.namesforfields(myfieldId)
-6932
+        if (type(fieldName) == list or type(fieldName) == type(np.ndarray(0))):
-6933
+            fieldName = fieldName[0]
-6934
+#        print("A) fieldname = ", fieldName)
-6935
+#        print("myfieldId = ", myfieldId)
--
+        myscantime = np.mean(mymsmd.timesforscan(scan))
-+
+        myscantime = np.median(mymsmd.timesforscan(scan))
-6937
+#        print("Calling getRADecForField")
-6938
+        mydirection = getRADecForField(vis, myfieldId, verbose)
-6939
+        if (verbose): print("mydirection= %s" % (str(mydirection)))
-6940
+        if (len(telescopeName) < 1):
-6941
+            telescopeName = 'ALMA'
-6942
+        myazel = computeAzElFromRADecMJD(mydirection, myscantime/86400., telescopeName)
-6943
+        conditions['elevation'] = myazel[1] * 180/math.pi
-6944
+        conditions['azimuth'] = myazel[0] * 180/math.pi
-6945
+        conditions['solarangle'] = angularSeparation(azsun,elsun,conditions['azimuth'],conditions['elevation'])
-6946
+        conditions['solarelev'] = elsun
-6963
+              print("Scans to loop over = %s" % (str(scan)))
-6964
+          for s in scan:
-6965
+              fieldName = mymsmd.fieldsforscan(s)
-6966
+              if (type(fieldName) == list):
-6967
+                  # take only the first pointing in the mosaic
-6968
+                  fieldName = fieldName[0]
-6969
+              myfieldId = mymsmd.fieldsforname(fieldName)
-6970
+              if (type(myfieldId) == list or type(myfieldId)==type(np.ndarray(0))):
-6971
+                  # If the same field name has two IDs (this happens in EVLA data)
-6972
+                  myfieldId = myfieldId[0]
--
+              myscantime = np.mean(mymsmd.timesforscan(s))
-+
+              myscantime = np.median(mymsmd.timesforscan(s))
-6974
+              mydirection = getRADecForField(vis, myfieldId, verbose)
-6975
+              telescopeName = mymsmd.observatorynames()[0]
-6976
+              if (len(telescopeName) < 1):
-6977
+                  telescopeName = 'ALMA'
-6978
+              myazel = computeAzElFromRADecMJD(mydirection, myscantime/86400., telescopeName)
-6979
+              myaz.append(myazel[0]*180/math.pi)
-6980
+              myel.append(myazel[1]*180/math.pi)
--
+          conditions['azimuth'] = np.mean(myaz)
--
+          conditions['elevation'] = np.mean(myel)
-+
+          conditions['azimuth'] = np.median(myaz)
-+
+          conditions['elevation'] = np.median(myel)
-6983
+          conditions['solarangle'] = angularSeparation(azsun,elsun,conditions['azimuth'],conditions['elevation'])
-6984
+          conditions['solarelev'] = elsun
-6985
+          conditions['solarazim'] = azsun
-6986
+          if (verbose):
--
+              print("Using antenna = %s to retrieve mean azimuth and elevation" % (antennaName))
-+
+              print("Using antenna = %s to retrieve median azimuth and elevation" % (antennaName))
-6988
+              print("Separation from sun = %f deg" % (abs(conditions['solarangle'])))
-6989
+          if (elsun<0):
-6990
+              conditions['solarangle'] = -conditions['solarangle']
-6991
+              if (verbose):
-6992
+                  print("Sun is below horizon (elev=%.1f deg)" % (elsun))
-6993
+          else:
-6994
+              if (verbose):
-6995
+                  print("Sun is above horizon (elev=%.1f deg)" % (elsun))
-6996
+          if (verbose):
-6997
+              print("Average azimuth = %.2f, elevation = %.2f degrees" % (conditions['azimuth'],conditions['elevation']))
-7008
+        print("Could not open the WEATHER table for this ms.")
-7009
+        mytb.done()
-7010
+        return([conditions,myTimes])
-7011
+    if (True):
-7012
+        mjdsec = mytb.getcol('TIME')
-7013
+        indices = np.argsort(mjdsec)
-7014
+        mjd = mjdsec/86400.
-7015
+        pressure = mytb.getcol('PRESSURE')
-7016
+        relativeHumidity = mytb.getcol('REL_HUMIDITY')
-7017
+        temperature = mytb.getcol('TEMPERATURE')
--
+        if (np.mean(temperature) > 100):
-+
+        if (np.median(temperature) > 100):
-7019
+            # must be in units of Kelvin, so convert to C
-7020
+            temperature -= 273.15
-7021
+        if 'DEW_POINT' in mytb.colnames():
-7022
+            dewPoint = mytb.getcol('DEW_POINT')
--
+            if (np.mean(dewPoint) > 100):
-+
+            if (np.median(dewPoint) > 100):
-7024
+                # must be in units of Kelvin, so convert to C
-7025
+                dewPoint -= 273.15
--
+            if (np.mean(dewPoint) == 0):
-+
+            if (np.median(dewPoint) == 0):
-7027
+                # assume it is not measured and use NOAA formula to compute from humidity:
-7028
+                dewPoint = ComputeDewPointCFromRHAndTempC(relativeHumidity, temperature)
-7029
+        else:
-7030
+            dewPoint = None  # Nobeyama measurement sets do not have a dewpoint column
-7031
+        sinWindDirection = np.sin(mytb.getcol('WIND_DIRECTION'))
-7032
+        cosWindDirection = np.cos(mytb.getcol('WIND_DIRECTION'))
-7033
+        windSpeed = mytb.getcol('WIND_SPEED')
-7034
+        mytb.done()
-7035
-7036
+        # put values into time order (they mostly are, but there can be small differences)
-7080
+              conditions['readings'] = 1
-7081
+              if (verbose):
-7082
+                  print("selectedValues=%d, myTimes[0]=%.0f, len(matches)=%d, len(matches2)=%d" % (selectedValues,
-7083
+                     myTimes[0], len(matches), len(matches2)))
-7084
+                  if (len(matches) > 0):
-7085
+                      print("matches[0]=%f, matches[-1]=%f" % (matches[0], matches[-1]))
-7086
+                  if (len(matches2) > 0):
-7087
+                      print("matches2[0]=%f, matches2[-1]=%d" % (matches2[0], matches2[-1]))
-7088
+          else:
-7089
+              conditions['readings'] = len(selectedValues)
--
+          conditions['pressure'] = np.mean(pressure[selectedValues])
-+
+          conditions['pressure'] = np.median(pressure[selectedValues])
-7091
+          if (conditions['pressure'] != conditions['pressure']):
-7092
+              # A nan value got through, due to no selected values (should be impossible)"
-7093
+              if (verbose):
-7094
+                  print(">>>>>>>>>>>>>>>>>>>>>>>>  selectedValues = %s" % (str(selectedValues)))
-7095
+                  print("len(matches)=%d, len(matches2)=%d" % (len(matches), len(matches2)))
-7096
+                  print("matches[0]=%f, matches[-1]=%f, matches2[0]=%f, matches2[-1]=%d" % (matches[0], matches[-1], matches2[0], matches2[-1]))
--
+          conditions['temperature'] = np.mean(temperature[selectedValues])
--
+          conditions['humidity'] = np.mean(relativeHumidity[selectedValues])
-+
+          conditions['temperature'] = np.median(temperature[selectedValues])
-+
+          conditions['humidity'] = np.median(relativeHumidity[selectedValues])
-7099
+          if dewPoint is not None:
--
+              conditions['dewpoint'] = np.nanmean(dewPoint[selectedValues])
--
+          conditions['windspeed'] = np.mean(windSpeed[selectedValues])
--
+          conditions['winddirection'] = (180./math.pi)*np.arctan2(np.mean(sinWindDirection[selectedValues]),np.mean(cosWindDirection[selectedValues]))
-+
+              conditions['dewpoint'] = np.nanmedian(dewPoint[selectedValues])
-+
+          conditions['windspeed'] = np.median(windSpeed[selectedValues])
-+
+          conditions['winddirection'] = (180./math.pi)*np.arctan2(np.median(sinWindDirection[selectedValues]),np.median(cosWindDirection[selectedValues]))
-7103
+          if (conditions['winddirection'] < 0):
-7104
+              conditions['winddirection'] += 360
-7105
+          if (verbose):
--
+              print("Mean weather values for scan %s (field %s)" % (listscan,listfield))
-+
+              print("Median weather values for scan %s (field %s)" % (listscan,listfield))
-7107
+              print("  Pressure = %.2f mb" % (conditions['pressure']))
-7108
+              print("  Temperature = %.2f C" % (conditions['temperature']))
-7109
+              if dewPoint is not None:
-7110
+                  print("  Dew point = %.2f C" % (conditions['dewpoint']))
-7111
+              print("  Relative Humidity = %.2f %%" % (conditions['humidity']))
-7112
+              print("  Wind speed = %.2f m/s" % (conditions['windspeed']))
-7113
+              print("  Wind direction = %.2f deg" % (conditions['winddirection']))
-7114
-7115
+    return([conditions,myTimes])
-7116
+    # end of getWeather

casa6

Commits

Jira Issues

Add shortcut