WSL/SLF GitLab Repository

Commit eb5f9d5e authored by Mathias Bavay's avatar Mathias Bavay
Browse files

The STD_PRESS spatial interpolation algorithm has been extended and can now...

The STD_PRESS spatial interpolation algorithm has been extended and can now deal with multiple input stations (with 2 different methods for handling this case). The Imis plugin now reads the air pressure and convert it to station air pressure (instead of sea level). Otherwise, improved error messages in libsmet and documentation / minor code cleanup.
parent 898a35cb
......@@ -88,7 +88,7 @@ size_t InterpolationAlgorithm::getData(const Date& i_date, const MeteoData::Para
for (size_t ii=0; ii<vecMeteo.size(); ii++){
const double& val = vecMeteo[ii](i_param);
if (val != IOUtils::nodata) {
o_vecData.push_back(val);
o_vecData.push_back( val );
}
}
......@@ -104,8 +104,8 @@ size_t InterpolationAlgorithm::getData(const Date& i_date, const MeteoData::Para
for (size_t ii=0; ii<vecMeteo.size(); ii++){
const double& val = vecMeteo[ii](i_param);
if (val != IOUtils::nodata){
o_vecData.push_back(val);
o_vecMeta.push_back(vecMeteo[ii].meta);
o_vecData.push_back( val );
o_vecMeta.push_back( vecMeteo[ii].meta );
}
}
......@@ -119,7 +119,7 @@ size_t InterpolationAlgorithm::getStationAltitudes(const std::vector<StationData
for (size_t ii=0; ii<i_vecMeta.size(); ii++){
const double& alt = i_vecMeta[ii].position.getAltitude();
if (alt != IOUtils::nodata) {
o_vecData.push_back(alt);
o_vecData.push_back( alt );
}
}
......@@ -301,11 +301,22 @@ double StandardPressureAlgorithm::getQualityRating(const Date& i_date, const Met
date = i_date;
param = in_param;
nrOfMeasurments = getData(date, param, vecData, vecMeta);
const size_t nr_args = vecArgs.size();
if (nr_args>1)
throw InvalidArgumentException("Wrong number of arguments supplied for the "+algo+" algorithm", AT);
if (nr_args==1) {
if (IOUtils::strToUpper(vecArgs[0])=="USE_RESIDUALS") {
use_residuals = true;
} else
throw InvalidArgumentException("Unknown argument \""+vecArgs[0]+"\" supplied to the "+algo+" interpolation", AT);
}
if (param != MeteoData::P)
return 0.0;
if (nrOfMeasurments == 0)
if (nrOfMeasurments <=1 || use_residuals)
return 1.0;
return 0.1;
......@@ -314,6 +325,35 @@ double StandardPressureAlgorithm::getQualityRating(const Date& i_date, const Met
void StandardPressureAlgorithm::calculate(const DEMObject& dem, Grid2DObject& grid)
{
Interpol2D::stdPressure(dem, grid);
if (nrOfMeasurments==1 || !use_residuals) { //correct the locally measured offset to std pressure
double offset = 0.;
size_t count = 0;
for (size_t ii=0; ii<nrOfMeasurments; ii++) {
const double altitude = vecMeta[ii].position.getAltitude();
if (altitude!=IOUtils::nodata) {
offset += vecData[ii] - Atmosphere::stdAirPressure( altitude );
count++;
}
}
if (count>0) {
offset /= static_cast<double>( count );
grid += offset;
}
} else if (use_residuals && nrOfMeasurments>1) { //spatially distribute the residuals
std::vector<double> residuals;
for (size_t ii=0; ii<nrOfMeasurments; ii++) {
const double altitude = vecMeta[ii].position.getAltitude();
if (altitude!=IOUtils::nodata)
residuals.push_back( vecData[ii] - Atmosphere::stdAirPressure( altitude ) );
}
if (residuals.empty())
throw IOException("Not enough data for spatially interpolating parameter " + MeteoData::getParameterName(param), AT);
Grid2DObject offset;
Interpol2D::IDW(residuals, vecMeta, dem, offset);
grid += offset;
}
}
......@@ -331,7 +371,8 @@ double ConstAlgorithm::getQualityRating(const Date& i_date, const MeteoData::Par
return 0.01;
}
void ConstAlgorithm::calculate(const DEMObject& dem, Grid2DObject& grid) {
void ConstAlgorithm::calculate(const DEMObject& dem, Grid2DObject& grid)
{
Interpol2D::constant(user_cst, dem, grid);
}
......@@ -351,7 +392,8 @@ double AvgAlgorithm::getQualityRating(const Date& i_date, const MeteoData::Param
return 0.0;
}
void AvgAlgorithm::calculate(const DEMObject& dem, Grid2DObject& grid) {
void AvgAlgorithm::calculate(const DEMObject& dem, Grid2DObject& grid)
{
Interpol2D::constant(Interpol1D::arithmeticMean(vecData), dem, grid);
}
......@@ -1363,10 +1405,10 @@ double SWRadInterpolation::getQualityRating(const Date& i_date, const MeteoData:
date = i_date;
param = in_param;
const size_t nrArgs = vecArgs.size();
if ( nrArgs>1 ) {
const size_t nr_args = vecArgs.size();
if ( nr_args>1 ) {
throw InvalidArgumentException("Wrong arguments supplied to the "+algo+" interpolation.", AT);
} else if (nrArgs==1) {
} else if (nr_args==1) {
if (IOUtils::strToUpper(vecArgs[0])=="NO_SHADING") {
shading = false;
} else
......
......@@ -216,16 +216,29 @@ class NoneAlgorithm : public InterpolationAlgorithm {
/**
* @class StandardPressureAlgorithm
* @brief Standard atmospheric pressure interpolation algorithm.
* Fill the grid with the standard atmosphere's pressure, depending on the local elevation.
* This first fills the grid with the standard atmosphere's pressure, depending on the local elevation. Then, depending on the available data:
* - if there are no measured atmospheric pressure, nothing else happens;
* - if one station has measured local atmospheric pressure, its offset to the standard atmospheric pressure is computed and applied to
* the computed grid;
* - if multiple stations have measured local atmospheric pressure:
* - default: the average offset will be applied to the computed grid;
* - USE_RESIDUALS option: the residuals are computed at each station, spatially distributed (with IDW) and applied to the computed grid;
*
* @code
* P::algorithms = STD_PRESS
* P::Std_Press = USE_RESIDUALS
* @endcode
*/
class StandardPressureAlgorithm : public InterpolationAlgorithm {
public:
StandardPressureAlgorithm(Meteo2DInterpolator& i_mi,
const std::vector<std::string>& i_vecArgs,
const std::string& i_algo, TimeSeriesManager& i_tsmanager, GridsManager& i_gridsmanager)
: InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, i_tsmanager, i_gridsmanager) {}
: InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, i_tsmanager, i_gridsmanager), use_residuals(false) {}
virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param);
virtual void calculate(const DEMObject& dem, Grid2DObject& grid);
private:
bool use_residuals; ///< should we compute residuals ate each station and distribute them spatially?
};
/**
......
......@@ -242,7 +242,7 @@ double Atmosphere::WBGT_index(const double& TA, const double& RH, const double&
}
/**
* @brief Standard water vapor saturation
* @brief Standard water vapor saturation pressure
* @param T air temperature (K)
* @return standard water vapor saturation pressure (Pa)
*/
......
......@@ -17,6 +17,7 @@
*/
#include <meteoio/plugins/ImisIO.h>
#include <meteoio/meteoLaws/Meteoconst.h>
#include <meteoio/meteoLaws/Atmosphere.h> //for p_sea to p_local conversion
#include <meteoio/MathOptim.h>
using namespace std;
......@@ -67,9 +68,9 @@ const string ImisIO::sqlQueryStationMetaData = "SELECT stao_name, stao_x, stao_y
const string ImisIO::sqlQuerySensorDepths = "SELECT hts1_1, hts1_2, hts1_3 FROM station2.v_station_standort WHERE stat_abk LIKE :1 AND stao_nr=:2"; ///< Sensor depths at station
const string ImisIO::sqlQueryMeteoDataDrift = "SELECT TO_CHAR(a.datum, 'YYYY-MM-DD HH24:MI') AS thedate, a.ta, a.iswr, a.vw, a.dw, a.vw_max, a.rh, a.ilwr, a.hnw, a.tsg, a.tss, a.hs, a.rswr, b.vw AS vw_drift, b.dw AS dw_drift, a.ts1, a.ts2, a.ts3 FROM (SELECT * FROM ams.v_ams_raw WHERE stat_abk=:1 AND stao_nr=:2 AND datum>=:3 AND datum<=:4) a LEFT OUTER JOIN (SELECT case when to_char(datum,'MI')=40 then trunc(datum,'HH24')+0.5/24 else datum end as datum, vw, dw FROM ams.v_ams_raw WHERE stat_abk=:5 AND stao_nr=:6 AND datum>=:3 AND datum<=:4) b ON a.datum=b.datum ORDER BY thedate"; ///< C. Marty's Data query with wind drift station; gets wind from enet stations for imis snow station too! [2010-02-24]
const string ImisIO::sqlQueryMeteoDataDrift = "SELECT TO_CHAR(a.datum, 'YYYY-MM-DD HH24:MI') AS thedate, a.ta, a.iswr, a.vw, a.dw, a.vw_max, a.rh, a.ilwr, a.hnw, a.tsg, a.tss, a.hs, a.rswr, a.ap, b.vw AS vw_drift, b.dw AS dw_drift, a.ts1, a.ts2, a.ts3 FROM (SELECT * FROM ams.v_ams_raw WHERE stat_abk=:1 AND stao_nr=:2 AND datum>=:3 AND datum<=:4) a LEFT OUTER JOIN (SELECT case when to_char(datum,'MI')=40 then trunc(datum,'HH24')+0.5/24 else datum end as datum, vw, dw FROM ams.v_ams_raw WHERE stat_abk=:5 AND stao_nr=:6 AND datum>=:3 AND datum<=:4) b ON a.datum=b.datum ORDER BY thedate"; ///< C. Marty's Data query with wind drift station; gets wind from enet stations for imis snow station too! [2010-02-24]
const string ImisIO::sqlQueryMeteoData = "SELECT TO_CHAR(datum, 'YYYY-MM-DD HH24:MI') AS thedate, ta, iswr, vw, dw, vw_max, rh, ilwr, hnw, tsg, tss, hs, rswr, ts1, ts2, ts3 FROM ams.v_ams_raw WHERE stat_abk=:1 AND stao_nr=:2 AND datum>=:3 AND datum<=:4 ORDER BY thedate ASC"; ///< Data query without wind drift station
const string ImisIO::sqlQueryMeteoData = "SELECT TO_CHAR(datum, 'YYYY-MM-DD HH24:MI') AS thedate, ta, iswr, vw, dw, vw_max, rh, ilwr, hnw, tsg, tss, hs, rswr, ap, ts1, ts2, ts3 FROM ams.v_ams_raw WHERE stat_abk=:1 AND stao_nr=:2 AND datum>=:3 AND datum<=:4 ORDER BY thedate ASC"; ///< Data query without wind drift station
const string ImisIO::sqlQuerySWEData = "SELECT TO_CHAR(datum, 'YYYY-MM-DD HH24:MI') AS thedate, swe FROM snowpack.ams_pmod WHERE stat_abk=:1 AND stao_nr=:2 AND datum>=:3 AND datum<=:4 ORDER BY thedate ASC"; ///< Query SWE as calculated by SNOWPACK to feed into PSUM
......@@ -748,8 +749,9 @@ void ImisIO::parseDataSet(const std::vector<std::string>& i_meteo, MeteoData& md
IOUtils::convertString(md(MeteoData::TSS), i_meteo.at(10), std::dec);
IOUtils::convertString(md(MeteoData::HS), i_meteo.at(11), std::dec);
IOUtils::convertString(md(MeteoData::RSWR), i_meteo.at(12), std::dec);
IOUtils::convertString(md(MeteoData::P), i_meteo.at(13), std::dec);
unsigned int ii = 13;
unsigned int ii = 14;
if (fullStation) {
if (!md.param_exists("VW_DRIFT")) md.addParameter("VW_DRIFT");
IOUtils::convertString(md("VW_DRIFT"), i_meteo.at(ii++), std::dec);
......@@ -895,7 +897,7 @@ size_t ImisIO::getStationMetaData(const std::string& stat_abk, const std::string
/**
* @brief Gets data from ams.v_ams_raw which is a table of SDB and
* retrieves the temperature sensor depths from station2.standort \n
* retrieves the temperature sensor depths from station2.standort.
* Each record returned are vector of strings which are pushed back in vecMeteoData.
* @param stat_abk : a string key of ams.v_ams_raw
* @param stao_nr : a string key of ams.v_ams_raw
......@@ -984,7 +986,7 @@ void ImisIO::convertUnits(MeteoData& meteo)
{
meteo.standardizeNodata(plugin_nodata);
//converts C to Kelvin, converts RH to [0,1]
//converts C to Kelvin, converts RH to [0,1], HS to m and P to Pa
double& ta = meteo(MeteoData::TA);
ta = IOUtils::C_TO_K(ta);
......@@ -1001,6 +1003,13 @@ void ImisIO::convertUnits(MeteoData& meteo)
double& hs = meteo(MeteoData::HS);
if (hs != IOUtils::nodata)
hs /= 100.0;
double& p = meteo(MeteoData::P);
if (p != IOUtils::nodata) {
p *= 100.0;
if (ta==IOUtils::nodata) p=IOUtils::nodata;
else p *= Atmosphere::stdAirPressure(meteo.meta.position.getAltitude()) / Cst::std_press;
}
//convert extra parameters (if present) //HACK TODO: find a dynamic way...
convertSnowTemperature(meteo, "TS1");
......
......@@ -548,7 +548,7 @@ void PNGIO::write2DGrid(const Grid2DObject& grid_in, const MeteoGrids::Parameter
gradient.set(Gradient::bg_isomorphic, min, max, autoscale);
} else if (parameter==MeteoGrids::P) {
if (!autoscale) {
//lowest and highest pressures ever recorded on Earth: 87000 and 108570
//lowest and highest sea level pressures ever recorded on Earth: 87000 and 108570
min = 87000.; max = 115650.; //centered around 1 atm
gradient.set(Gradient::bluewhitered, min, max, autoscale);
} else {
......
......@@ -346,13 +346,13 @@ SMETWriter::SMETWriter(const std::string& in_filename, const std::string& in_fie
if (nr_of_fields!=vecFields.size()) {
std::ostringstream ss;
ss << "Trying to write " << vecFields.size() << " fields in file '" << filename << "' that has " << nr_of_fields << " fields";
ss << "Trying to append " << vecFields.size() << " fields in existing file '" << filename << "' that has " << nr_of_fields << " fields";
throw SMETException(ss.str(), AT);
}
for (size_t ii=0; ii<nr_of_fields; ii++) {
if (vecFields[ii] != reader.get_field_name(ii)) {
std::ostringstream ss;
ss << "Trying to write field '" << vecFields[ii] << "' at position " << ii << " in file '" << filename << "' when it should be '" << reader.get_field_name(ii);
ss << "Trying to append field '" << vecFields[ii] << "' at position " << ii << " in existing file '" << filename << "' when it should be '" << reader.get_field_name(ii);
throw SMETException(ss.str(), AT);
}
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment