casa / casa5 / d62d2389590

Commits

Pam Ford authored d62d2389590 12 Apr 2018

Select weather table to remove invalid values and ALMA station

code/plotms/Data/PlotMSAtm.cc

Modified

-357
+    if (pwv == 0.0) {
-358
+        if (telescopeName_=="ALMA") pwv = 1.0;
-359
+        else pwv = 5.0;
-360
+        parent_->logmesg("load_cache",
-361
+            "Using default pwv " + casacore::String::toString(pwv) + " for telescope " + telescopeName_);
-362
+    }
-363
+    pwv_ = pwv;
-364
+}
-365
-366
+void PlotMSAtm::getMeanWeather() {
--
+    // Info from MS WEATHER table in weather_ Record
--
+    // set defaults
--
+    casacore::Float pressure, humidity(20.0), temperature(273.15);
--
+    // NB: plotbandpass uses default pressure 563 in all cases;
--
+    // see CAS-9053 algorithm #2
--
+    pressure = (telescopeName_=="ALMA" ? 563.0 : 786.0);  // mb
--
--
+    // values from WEATHER table if it exists
--
+    bool noWeather(false);
-+
+    // Fill weather_ Record with mean values from MS WEATHER table
-+
+    // Set defaults; CAS-9053 default pressure depends on telescope
-+
+    casacore::Float humidity(20.0), temperature(273.15);
-+
+    casacore::Float pressure = (telescopeName_=="ALMA" ? 563.0 : 786.0);  // mb
-+
-+
+    // Use values from WEATHER table if it exists
-+
+    bool noWeather(true);
-374
+    if (!tableName_.empty()) {
--
+        casacore::Table mstab(tableName_);
--
+        if (mstab.keywordSet().fieldNumber("WEATHER") > -1) {
--
+            casacore::String subname = tableName_ + "::WEATHER";
--
+            try {
--
+                casacore::Table subtable = Table::openTable(subname);
--
+                // get columns and temp. units
--
+                casacore::Vector<casacore::Float> pressureCol, humidityCol,
--
+                    temperatureCol;
--
+                casacore::Vector<casacore::Double> timeCol;
--
+                pressureCol = ScalarColumn<casacore::Float>(subtable, "PRESSURE").getColumn();
--
+                humidityCol = ScalarColumn<casacore::Float>(subtable, "REL_HUMIDITY").getColumn();
--
+                ScalarColumn<casacore::Float> tempColumn =
--
+                    ScalarColumn<casacore::Float>(subtable, "TEMPERATURE");
--
+                temperatureCol = tempColumn.getColumn();
--
+                casacore::String units =
--
+                    tempColumn.keywordSet().asArrayString("QuantumUnits").tovector()[0];
--
+                timeCol = ScalarColumn<casacore::Double>(subtable, "TIME").getColumn();
--
+                mstab.closeSubTables();
--
--
+                // now select based on cal times
--
+                casacore::Vector<casacore::Float> selpressure, selhumidity,
--
+                    seltemperature;
--
+                // pressure
--
+                if (!pressureCol.empty()) {
--
+                    selpressure = getValuesNearTimes(pressureCol, timeCol);
--
+                    if (!selpressure.empty()) pressure = mean(selpressure);
--
+                }
--
+                // humidity
--
+                if (!humidityCol.empty()) {
--
+                    selhumidity = getValuesNearTimes(humidityCol, timeCol);
--
+                    if (!selhumidity.empty()) humidity = mean(selhumidity);
--
+                }
--
+                // temperature
--
+                if (!temperatureCol.empty()) {
--
+                    seltemperature = getValuesNearTimes(temperatureCol, timeCol);
--
+                    if (!seltemperature.empty()) {
--
+                        temperature = mean(seltemperature);
--
+                        if (units.compare("K")!=0)
--
+                            temperature += (float)273.15;  // convert C to K
--
+                    }
-+
+        try {
-+
+            // openTable throws exception if table doesn't exist
-+
+            casacore::Table wtable = Table::openTable(tableName_ + "::WEATHER");
-+
+            TableColumn tempCol = TableColumn(wtable, "TEMPERATURE");
-+
+            String tempUnits = tempCol.keywordSet().asArrayString("QuantumUnits").tovector()[0];
-+
+            // select valid stations and values in range
-+
+            casacore::Table selwtable = selectWeatherTable(wtable, tempUnits);
-+
-+
+            // get columns
-+
+            casacore::Vector<casacore::Float>
-+
+                pressureCol(ScalarColumn<casacore::Float>(selwtable, "PRESSURE").getColumn()),
-+
+                humidityCol(ScalarColumn<casacore::Float>(selwtable, "REL_HUMIDITY").getColumn()),
-+
+                temperatureCol(ScalarColumn<casacore::Float>(selwtable, "TEMPERATURE").getColumn());
-+
+            casacore::Vector<casacore::Double>
-+
+                timeCol(ScalarColumn<casacore::Double>(selwtable, "TIME").getColumn());
-+
-+
+            // select values based on cal times
-+
+            casacore::Vector<casacore::Float> selpressure, selhumidity, seltemperature;
-+
+            casacore::Float meanp(0.0), meanh(0.0), meant(0.0);
-+
+            // pressure
-+
+            if (!pressureCol.empty()) {
-+
+                selpressure = getValuesNearTimes(pressureCol, timeCol);
-+
+                if (!selpressure.empty())
-+
+                    meanp = mean(selpressure);
-+
+                if (meanp==0.0)
-+
+                    parent_->logmesg("load_cache", "WEATHER pressure is zero, using default value instead.");
-+
+                else
-+
+                    pressure = meanp;
-+
+            }
-+
+            // humidity
-+
+            if (!humidityCol.empty()) {
-+
+                selhumidity = getValuesNearTimes(humidityCol, timeCol);
-+
+                if (!selhumidity.empty())
-+
+                    meanh = mean(selhumidity);
-+
+                if (meanh==0.0)
-+
+                    parent_->logmesg("load_cache", "WEATHER humidity is zero, using default value instead.");
-+
+                else
-+
+                    humidity = meanh;
-+
+            }
-+
-+
+            // temperature
-+
+            if (!temperatureCol.empty()) {
-+
+                seltemperature = getValuesNearTimes(temperatureCol, timeCol);
-+
+                if (!seltemperature.empty())
-+
+                    meant = mean(seltemperature);
-+
+                if (meant==0.0) {
-+
+                    parent_->logmesg("load_cache", "WEATHER temperature is zero, using default value instead.");
-+
+                } else {
-+
+                    if (tempUnits=="C") meant += (float)273.15;  // convert C to K
-+
+                    temperature = meant;
-425
+                }
--
+            } catch (AipsError & err) {
--
+                noWeather = true;
-426
+            }
--
+        } else {
--
+            noWeather = true;
-+
+            noWeather = false;
-+
+        } catch (AipsError & err) {
-+
+            cout << err.getMesg() << endl;
-430
+        }
--
+    } else {
--
+        noWeather = true;
--
+    }
--
+    if (noWeather) {
--
+        parent_->logmesg("load_cache", "Could not open WEATHER table, using default values: humidity " + casacore::String::toString(humidity) + ", pressure " + casacore::String::toString(pressure) + ", temperature " + casacore::String::toString(temperature) + " for Atm/Tsky");
--
+    }
--
+    // reset to defaults
--
+    if (humidity==0.0) {
--
+        humidity = 20.0;
--
+        parent_->logmesg("load_cache", "WEATHER table humidity is 0.0, using default value " + casacore::String::toString(humidity));
--
+    }
--
+    if (pressure==0.0) {
--
+        pressure = (telescopeName_=="ALMA" ? 563.0 : 786.0);  // mb
--
+        parent_->logmesg("load_cache", "WEATHER table pressure is 0.0, using default value " + casacore::String::toString(pressure));
--
+    }
--
+    if (temperature==0.0) {
--
+        temperature = 273.15;
--
+        parent_->logmesg("load_cache", "WEATHER table temperature is 0.0, using default value " + casacore::String::toString(temperature));
-431
+    }
-+
-+
+	String msg = "Using";
-+
+    msg += (noWeather ? " default " : " ");
-+
+	msg += "weather values for Atm/Tsky:\n";
-+
+    msg += "  humidity " + casacore::String::toString(humidity) +
-+
+        ", pressure " + casacore::String::toString(pressure) +
-+
+        ", temperature " + casacore::String::toString(temperature);
-+
+    parent_->logmesg("load_cache", msg);
-440
+    // to use in atmosphere.initAtmProfile (tool)
-441
+    weather_.define("humidity", humidity);       // %
-442
+    weather_.define("pressure", pressure);       // mb
-443
+    weather_.define("temperature", temperature); // K
-444
+}
-445
-+
+Table PlotMSAtm::selectWeatherTable(Table& wtable, String tempUnits) {
-+
+    // Remove invalid station, values from weather table
-+
+    TableExprNode ten;
-+
+    if (telescopeName_=="ALMA") {
-+
+        try {
-+
+            // do not use weather station "MeteoItinerant"
-+
+            casacore::Table staTable = Table::openTable(tableName_ + "::ASDM_STATION");
-+
+            casacore::Table result = staTable(staTable.col("name")=="MeteoItinerant");
-+
+            Vector<uInt> antIds = result.rowNumbers(staTable);
-+
+            if (!antIds.empty()) ten = (wtable.col("NS_WX_STATION_ID") != antIds(0));
-+
+        } catch (AipsError & err) {
-+
+            // table does not exist
-+
+        }
-+
+    }
-+
+    // SCIREQ-1241
-+
+    // temperature in the range -30 to +30C
-+
+    // pressure in the range 450 to 800mb
-+
+    // Humidity in the range -10 to +110%
-+
+    Float mintemp = (tempUnits=="C" ? -30.0 : 243.15);
-+
+    Float maxtemp = (tempUnits=="C" ? 30.0 : 303.15);
-+
+    ten = ten && (wtable.col("TEMPERATURE")>=mintemp && wtable.col("TEMPERATURE")<=maxtemp);
-+
+    ten = ten && (wtable.col("PRESSURE")>=450.0 && wtable.col("PRESSURE")<=800.0);
-+
+    ten = ten && (wtable.col("REL_HUMIDITY")>=-10.0 && wtable.col("REL_HUMIDITY")<=110.0);
-+
+    casacore::Table selwtable(wtable(ten));
-+
+    return selwtable;
-+
+}
-+
-473
+casacore::Double PlotMSAtm::computeMeanAirmass() {
-474
+    // Calculate airmass from elevation
-475
+    casacore::Vector<casacore::Double> airmasses;
-476
+    airmasses.resize(fields_.size());
-477
+    casacore::Double elevation;
-478
+    for (uInt i=0; i<fields_.size(); ++i) {
-479
+        elevation = getElevation(fields_(i));
-480
+        if (elevation <= 3.0) elevation = 45.0;
-481
+        airmasses(i) = 1.0 / std::cos((90.0 - elevation) * C::pi / 180.0);
-482
+    }

casa5

Commits

Jira Issues

Add shortcut