WSL/SLF GitLab Repository

Commit 9b9ba4c3 authored by Mathias Bavay's avatar Mathias Bavay
Browse files

A bug has been found in the processing of the variance: we were using a naive...

A bug has been found in the processing of the variance: we were using a naive formula that is ill-fitted to signals showing small variations around a large mean. This lead to catastrophic cancellation and negative values... This has been replaced by a (slower) compensated formula (that works in two passes). The standard deviation filter has been fixed. 

The Tukey filter has been brought back into the new infrastructure and should now be usable. However it is not suitable to snow heights: in such a case, the standard deviation is too small (or even zero) and the signal might show some sudden variations (like in the case of a snow fall) that would therefore be falsely rejected. This filter seems much better suited to wind speeds that exhibit some constant background noise, thus a non-null standard deviation.

The rate filter can now take 2 arguments: in such a case, the first one is the minimum rate and the second the maximum rate. This would allow to define a descent rate different from the climb rate.

Several error messages have been improved in A3DIO and general plugin loading.

Finally, two small tools for working with smet files have been contributed, one allowing transparent loading of smet files in xmgrace (smet2agr). And the example io.ini has been fixed (it had been broken and was not working with the examples anymore)
parent f643892b
......@@ -142,11 +142,16 @@ INSTALL(FILES ${hfilterfiles} DESTINATION include/meteoio/meteofilters COMPONENT
IF(WIN32) #handle the special case of WIN32 needing import libraries
FILE(GLOB importlibs "${CMAKE_CURRENT_BINARY_DIR}/libmeteoio*.lib*")
FILE(GLOB solibs "${CMAKE_BINARY_DIR}/meteoio/libmeteoio*.${DLL_EXT}*")
FILE(GLOB plugins "${CMAKE_BINARY_DIR}/meteoio/plugins/lib*.${DLL_EXT}*")
FILE(GLOB alibs "${CMAKE_BINARY_DIR}/meteoio/*.${STAT_EXT}*")
ELSE(WIN32)
FILE(GLOB solibs "${CMAKE_BINARY_DIR}/lib/libmeteoio*.${DLL_EXT}*")
FILE(GLOB plugins "${CMAKE_BINARY_DIR}/lib/plugins/lib*.${DLL_EXT}*")
FILE(GLOB alibs "${CMAKE_BINARY_DIR}/lib/*.${STAT_EXT}*")
ENDIF(WIN32)
FILE(GLOB solibs "${CMAKE_BINARY_DIR}/lib/libmeteoio*.${DLL_EXT}*")
FILE(GLOB plugins "${CMAKE_BINARY_DIR}/lib/plugins/lib*.${DLL_EXT}*")
FILE(GLOB alibs "${CMAKE_BINARY_DIR}/lib/*.${STAT_EXT}*")
IF(DEST MATCHES "par")
FILE(GLOB modules "${CMAKE_BINARY_DIR}/lib/*.module")
INSTALL(FILES ${modules} DESTINATION lib/meteoio COMPONENT libraries)
......
......@@ -28,6 +28,8 @@ int main(int /*argc*/, char** argv) {
io.write2DGrid(param,"hnw.asc");
io.interpolate(d1, dem, MeteoData::RH, param);
io.write2DGrid(param,"rh.asc");
io.interpolate(d1, dem, MeteoData::ILWR, param);
io.write2DGrid(param,"ilwr.asc");
return 0;
}
[General]
PLUGINPATH = ../../lib/plugins
BUFF_CHUNK_SIZE = 15
BUFF_BEFORE = 1.5
PLUGINPATH = ../../lib/plugins
BUFF_CHUNK_SIZE = 30
BUFF_BEFORE = 1.5
[Input]
COORDSYS = CH1903
COORDPARAM = NULL
TIME_ZONE = 1
; COORDSYS = PROJ4
; COORDPARAM = 21781
COORDSYS = CH1903
TIME_ZONE = 1
; COORDSYS = PROJ4
; COORDPARAM = 21781
#reading ARC dem
; DEM = ARC
; DEMFILE = ./input/surface-grids/Switzerland_1000m.asc
;
; GRID2D = ARC
DEM = ARC
DEMFILE = ./input/surface-grids/Switzerland_1000m.asc
GRID2D = ARC
#reading ARPS dem
; DEM = ARPS
; DEMFILE = ./wgrt10r2_vw4.asc
; ARPS_X = 653400
; ARPS_Y = 112204
; DEM = ARPS
; DEMFILE = ./wgrt10r2_vw4.asc
; ARPS_X = 653400
; ARPS_Y = 112204
#reading PGM image as dem
; DEM = PGM
; DEMFILE = ./Switzerland.pgm
; PGM_XCOORD = 479500.
; PGM_YCOORD = 73500.
; PGM_CELLSIZE = 1000.
; PGM_MIN = 193.
; PGM_MAX = 4204.
; DEM = PGM
; DEMFILE = ./Switzerland.pgm
; PGM_XCOORD = 479500.
; PGM_YCOORD = 73500.
; PGM_CELLSIZE = 1000.
; PGM_MIN = 193.
; PGM_MAX = 4204.
; LANDUSE = ARC
; LANDUSEFILE = ch_as97.asc
; LANDUSE = ARC
; LANDUSEFILE = ch_as97.asc
#Alpine3D traditional inputs -> A3D plugin
; METEO = A3D
; METEOPATH = ./input/meteo
; METEO = A3D
; METEOPATH = ./input/meteo
#Borma
; METEO = BORMA
; METEOPATH = ./input/xml
; NROFSTATIONS = 1
; STATION1 = 00.00.00.2
; METEO = BORMA
; METEOPATH = ./input/xml
; NROFSTATIONS = 1
; STATION1 = 00.00.00.2
#Snowpack input
; METEO = SNOWPACK
; METEOPATH = input
; NROFSTATIONS = 1
; METAFILE = IMIS_Extracted_Info.txt ;metadata for all stations
; STATION1 = DAV1 ;this is used as a stationID to get meta data from METAFILE
; METEOFILE1 = MST96_RR.inp
; METEO = SNOWPACK
; METEOPATH = input
; NROFSTATIONS = 1
; METAFILE = IMIS_Extracted_Info.txt ;metadata for all stations
; STATION1 = DAV1 ;this is used as a stationID to get meta data from METAFILE
; METEOFILE1 = MST96_RR.inp
#SMET meteorological file format
METEO = SMET ;SNOWPACK ;
METEOPATH = ./input/meteo
NROFSTATIONS = 7
STATION1 = FLU2.smet
STATION2 = FIR2.smet
STATION3 = FRA2.smet
STATION4 = GLA2.smet
STATION5 = ILI2.smet
STATION6 = OTT2.smet
STATION7 = TUJ3.smet
METEO = SMET
METEOPATH = ./input/meteo
STATION1 = FLU2.smet
STATION2 = FIR2.smet
STATION3 = FRA2.smet
STATION4 = GLA2.smet
STATION5 = ILI2.smet
STATION6 = OTT2.smet
STATION7 = TUJ3.smet
#IMIS network database input -> IMIS plugin
; METEO = IMIS
; DBNAME = sdbo
; DBUSER = XXX
; DBPASS = XXX
; NROFSTATIONS = 6
; STATION1 = DAV2
; STATION2 = KLO3
; STATION3 = BED1
; STATION4 = WFJ1
; STATION5 = *SAM0
; STATION6 = MORN2
; METEO = IMIS
; DBNAME = sdbo
; DBUSER = XXX
; DBPASS = XXX
; NROFSTATIONS = 4
; STATION1 = MORN2
; STATION2 = DAV3
; STATION3 = KLO2
; STATION4 = *SAM0
#GEOtop traditional inputs -> GEOTOP plugin
; METEO = GEOTOP
; METEOPATH = meteo/
; METEOPREFIX = _meteo
; METEO = GEOTOP
; METEOPATH = meteo/
; METEOPREFIX = _meteo
#GSN direct input -> GSN plugin
; METEO = GSN
; ENDPOINT = http://montblanc.slf.ch:22001/services/A3DWebService
; NROFSTATIONS = 4
; STATION1 = wan_sen14_2008
; STATION2 = lafouly_st_1043
; STATION3 = lafouly_st_1042
; STATION4 = lafouly_st_1040
; METEO = GSN
; ENDPOINT = http://montblanc.slf.ch:22001/services/A3DWebService
; NROFSTATIONS = 4
; STATION1 = wan_sen14_2008
; STATION2 = lafouly_st_1043
; STATION3 = lafouly_st_1042
; STATION4 = lafouly_st_1040
[Output]
COORDSYS = CH1903
COORDPARAM = NULL
TIME_ZONE = 1
COORDSYS = CH1903
TIME_ZONE = 1
GRID2D = PGM
GRID2D = PGM
METEO = SMET ;SNOWPACK ;
METEOPATH = ./output
METEO = SMET
METEOPATH = ./
[Filters]
HNW::filter1 = min
HNW::arg1 = soft 0.0
TA::filter1 = min_max
TA::arg1 = 240 320
TA::filter1 = min_max
TA::arg1 = 240 320
RH::filter1 = min_max
RH::arg1 = 0.01 1.2
RH::filter2 = min_max
RH::arg2 = soft 0.05 1.0
RH::filter1 = min_max
RH::arg1 = 0.01 1.2
RH::filter2 = min_max
RH::arg2 = soft 0.05 1.0
HNW::filter1 = min
HNW::arg1 = -0.1
HNW::filter2 = min
HNW::arg2 = soft 0.
ISWR::filter1 = min_max
ISWR::arg1 = -10 1500
ISWR::filter2 = min_max
ISWR::arg2 = soft 0 1500
ISWR::filter1 = min_max
ISWR::arg1 = -10. 1500.
ISWR::filter2 = min
ISWR::arg2 = soft 0.
RSWR::filter1 = min_max
RSWR::arg1 = -10 1500
RSWR::filter2 = min_max
RSWR::arg2 = soft 0 1500
RSWR::filter1 = min_max
RSWR::arg1 = -10 1500
RSWR::filter2 = min
RSWR::arg2 = soft 0
#for TA between 240 and 320 K
ILWR::filter1 = min_max
ILWR::arg1 = 188 600
ILWR::filter2 = min_max
ILWR::arg2 = soft 200 400
ILWR::filter1 = min_max
ILWR::arg1 = 188 600
ILWR::filter2 = min_max
ILWR::arg2 = soft 200 400
#we need to consider time with no snow -> TSS>0
#min(TSS) in db since 1998: -50C
TSS::filter1 = min_max
TSS::arg1 = 200 320
TSS::filter1 = min_max
TSS::arg1 = 200 320
#idem
TSG::filter1 = min_max
TSG::arg1 = 200 320
TSG::filter1 = min_max
TSG::arg1 = 200 320
HS::filter1 = min
HS::arg1 = soft 0.0
HS::filter2 = rate
HS::arg2 = 5.55e-5 ;0.20 m/h
HS::filter1 = min
HS::arg1 = soft 0.0
HS::filter2 = rate
HS::arg2 = 5.55e-5 ;0.20 m/h
VW::filter1 = min_max
VW::arg1 = -2 70
VW::filter2 = min_max
VW::arg2 = soft 0.2 50.0
VW::filter1 = min_max
VW::arg1 = -2 70
VW::filter2 = min_max
VW::arg2 = soft 0.2 50.0
[Interpolations1D]
WINDOW_SIZE = 86400
TA::resample = linear
TA::resample = linear
RH::resample = linear
RH::resample = linear
HS::resample = no
HS::resample = no
VW::resample = nearest_neighbour
VW::args = extrapolate
VW::resample = nearest_neighbour
VW::args = extrapolate
HNW::resample = linear
HNW::resample = linear
[Interpolations2D]
TA::algorithms = USER LIDW_LAPSE IDW_LAPSE CST_LAPSE
TA::cst_lapse = -0.008 soft
TA::lidw_lapse = 7
TA::user = ./test
TA::algorithms = USER LIDW_LAPSE IDW_LAPSE CST_LAPSE
TA::cst_lapse = -0.008 soft
TA::lidw_lapse = 7
TA::user = ./test
RH::algorithms = RH IDW_LAPSE CST
RH::algorithms = RH IDW_LAPSE CST
HNW::algorithms = IDW_LAPSE CST_LAPSE CST
HNW::cst_lapse = 0.0005 frac
HNW::algorithms = IDW_LAPSE CST_LAPSE CST
HNW::cst_lapse = 0.0005 frac
VW::algorithms = IDW_LAPSE CST
VW::algorithms = IDW_LAPSE CST
P::algorithms = STD_PRESS
P::algorithms = STD_PRESS
......@@ -242,15 +242,17 @@ void A3DIO::read1DStation(std::string& file_1d, StationData& sd)
try {
location.check();
} catch(...) {
std::cerr << "[E] Error in geographic coordinates in file " << file_1d << " trapped at " << AT << std::endl;
throw;
std::stringstream msg;
msg << "[E] Inconsistent geographic coordinates in file " << file_1d;
throw InvalidArgumentException(msg.str(), AT);
}
sd.setStationData(location, "meteo1d", "Meteo1D station");
} catch(...) {
std::cout << "[E] " << AT << ": "<< std::endl;
cleanup();
throw;
std::stringstream msg;
msg << "[E] Error processing header of file " << file_1d;
throw InvalidFormatException(msg.str(), AT);
}
}
......@@ -313,9 +315,10 @@ void A3DIO::read1DMeteo(const Date& dateStart, const Date& dateEnd, std::vector<
convertUnits(tmpdata);
}
} catch(...) {
std::cout << "[E] " << AT << ": "<< std::endl;
cleanup();
throw;
std::stringstream msg;
msg << "[E] Error processing data section of file " << file_1d << " possibly at line \"" << line << "\"";
throw InvalidFormatException(msg.str(), AT);
}
cleanup();
......@@ -552,8 +555,8 @@ unsigned int A3DIO::getNrOfStations(std::vector<std::string>& filenames, std::ma
}
void A3DIO::read2DMeteoData(const std::string& filename, const std::string& parameter,
std::map<std::string,unsigned int>& hashStations,
std::vector< std::vector<MeteoData> >& vecM, unsigned int& bufferindex)
std::map<std::string,unsigned int>& hashStations,
std::vector< std::vector<MeteoData> >& vecM, unsigned int& bufferindex)
{
std::string line_in = "";
......@@ -642,7 +645,7 @@ void A3DIO::read2DMeteoData(const std::string& filename, const std::string& para
}
void A3DIO::read2DMeteoHeader(const std::string& filename, std::map<std::string,unsigned int>& hashStations,
std::vector<StationData>& vecS)
std::vector<StationData>& vecS)
{
std::string line_in = "";
unsigned int columns = 0;
......
......@@ -146,7 +146,6 @@ void IOHandler::deletePlugin(DynamicLibrary*& dynLibrary, IOInterface*& io) thro
void IOHandler::loadPlugin(const std::string& libname, const std::string& classname, DynamicLibrary*& dynLibrary, IOInterface*& io)
{
std::cout << "[i] " << AT << ": Loading dynamic plugin: " << libname << std::endl;
std::string pluginpath = "";
try {
......@@ -155,20 +154,21 @@ void IOHandler::loadPlugin(const std::string& libname, const std::string& classn
pluginpath += "/";
//Which dynamic library needs to be loaded
std::cout << "\t" << "Trying to load " << libname << " ... ";
std::string filename = pluginpath + libname;
dynLibrary = DynamicLoader::loadObjectFile(filename);
if(dynLibrary == NULL) {
std::cout << "failed\n\tCouldn't load the dynamic library " << filename << "\n\t" << DynamicLoader::getErrorMessage() << std::endl;
std::cout << AT << ": [E] Failed loading dynamic plugin " << classname << " from " << filename << std::endl;
std::cout << "\t" << DynamicLoader::getErrorMessage() << std::endl;
std::cout << "Please check your PLUGINPATH in your configuration file!" << std::endl;
} else {
io = dynamic_cast<IOInterface*>((dynLibrary)->newObject(classname, cfg));
if(io == NULL) {
std::cout << "failed" << std::endl;
std::cout << AT << ": [E] Failed loading dynamic plugin " << classname << " from " << filename << "(NULL pointer to plugin's class)" << std::endl;
//delete dynLibrary; This causes a segfault !!
} else {
std::cout << "success" << std::endl;
std::cout << "[i] Success loading dynamic plugin " << classname << " from " << filename << std::endl;
}
}
} catch (std::exception& e) {
......
......@@ -131,7 +131,9 @@ double Interpol1D::getMedianAverageDeviation(const std::vector<double>& vecData)
}
double Interpol1D::variance(const std::vector<double>& X)
{
{//The variance is computed using a compensated variance algorithm,
//(see https://secure.wikimedia.org/wikipedia/en/wiki/Algorithms_for_calculating_variance)
//in order to be more robust to small variations around the mean.
const unsigned int n = X.size();
unsigned int count=0;
......@@ -147,8 +149,18 @@ double Interpol1D::variance(const std::vector<double>& X)
}
if(count<=1) return IOUtils::nodata;
const double mean = sum/(double)count;
return (sum_sq - sum*mean)/((double)count-1.);
double sum2=0., sum3=0.;
for(unsigned int i=0; i<n; i++) {
const double value=X[i];
if(value!=IOUtils::nodata) {
sum2 = sum2 + (value - mean)*(value - mean);
sum3 = sum3 + (value - mean);
}
}
const double variance = (sum2 - sum3*sum3/count) / (count - 1);
return variance;
}
double Interpol1D::std_dev(const std::vector<double>& X)
......@@ -157,7 +169,7 @@ double Interpol1D::std_dev(const std::vector<double>& X)
}
double Interpol1D::covariance(const std::vector<double>& X, const std::vector<double>& Y)
{
{//HACK: we should use a compensated formula here, similarly to the variance computation!
if(X.size()!=Y.size())
throw IOException("Vectors should have the same size for covariance!", AT);
const unsigned int n = X.size();
......
......@@ -9,6 +9,7 @@ SET(meteofilters_sources
meteofilters/FilterMeanAvg.cc
meteofilters/FilterWindAvg.cc
meteofilters/FilterStdDev.cc
meteofilters/FilterTukey.cc
meteofilters/RateFilter.cc
meteofilters/FilterBlock.cc
meteofilters/WindowedFilter.cc
......
......@@ -83,16 +83,15 @@ void FilterStdDev::process(const unsigned int& index, const std::vector<MeteoDat
}
}
void FilterStdDev::getStat(const std::vector<const MeteoData*>& vec_window, const unsigned int& index, double& stddev, double& mean) {
void FilterStdDev::getStat(const std::vector<const MeteoData*>& vec_window, const unsigned int& paramindex, double& stddev, double& mean) {
unsigned int count=0;
double sum=0., sum_sq=0.;
double sum=0.;
unsigned int n = vec_window.size();
for(unsigned int ii=0; ii<n; ii++) {
const double& value = (*vec_window[ii]).param(index);
const double& value = (*vec_window[ii]).param(paramindex);
if(value!=IOUtils::nodata) {
sum += value;
sum_sq += value*value;
count++;
}
}
......@@ -101,8 +100,17 @@ void FilterStdDev::getStat(const std::vector<const MeteoData*>& vec_window, cons
mean = IOUtils::nodata;
stddev = IOUtils::nodata;
} else {
//compensated variance algorithm, see https://secure.wikimedia.org/wikipedia/en/wiki/Algorithms_for_calculating_variance
mean = sum/(double)count;
const double variance = (sum_sq - sum*mean)/((double)count-1.);
double sum2=0., sum3=0.;
for(unsigned int ii=0; ii<n; ii++) {
const double& value = (*vec_window[ii]).param(paramindex);
if(value!=IOUtils::nodata) {
sum2 = sum2 + (value - mean)*(value - mean);
sum3 = sum3 + (value - mean);
}
}
const double variance = (sum2 - sum3*sum3/count) / (count - 1);
stddev = sqrt(variance);
}
}
......
......@@ -51,7 +51,7 @@ class FilterStdDev : public WindowedFilter {
private:
void parse_args(std::vector<std::string> vec_args);
void getStat(const std::vector<const MeteoData*>& vec_window, const unsigned int& index, double& stddev, double& mean);
void getStat(const std::vector<const MeteoData*>& vec_window, const unsigned int& paramindex, double& stddev, double& mean);
static const double sigma; ///<How many times the stddev allowed for valid points
};
......
/***********************************************************************************/
/* Copyright 2011 WSL Institute for Snow and Avalanche Research SLF-DAVOS */
/***********************************************************************************/
/* This file is part of MeteoIO.
MeteoIO is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
MeteoIO is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with MeteoIO. If not, see <http://www.gnu.org/licenses/>.
*/
#include <meteoio/meteofilters/FilterTukey.h>
#include <meteoio/IOUtils.h>
#include <cmath>
using namespace std;
namespace mio {
const double FilterTukey::k = 1.5; ///<How many times the stddev allowed as deviation to the smooth signal for valid points
FilterTukey::FilterTukey(const std::vector<std::string>& vec_args) : WindowedFilter("TUKEY")
{
parse_args(vec_args);
//This is safe, but maybe too imprecise:
properties.time_before = min_time_span;
properties.time_after = min_time_span;
properties.points_before = min_data_points;
properties.points_after = min_data_points;
}
void FilterTukey::process(const unsigned int& index, const std::vector<MeteoData>& ivec,
std::vector<MeteoData>& ovec)
{
ovec.clear();
for (unsigned int ii=0; ii<ivec.size(); ii++){ //for every element in ivec, get a window
ovec.push_back(ivec[ii]);
double& value = ovec[ii].param(index);
//Calculate std deviation
const double std_dev = getStdDev(ivec, ii, index); //HACK: we need to start at ii=0 for getWindow
const double u3 = getU3(ivec, ii, index);
if(std_dev!=IOUtils::nodata && u3!=IOUtils::nodata) {
if( abs(value-u3) > k*std_dev ) {
value = IOUtils::nodata;
}
}
}
}
double FilterTukey::getStdDev(const std::vector<MeteoData>& ivec, const unsigned int& ii, const unsigned int& paramindex) {
unsigned int count=0;
double sum=0.;
const vector<const MeteoData*>& vec_window = get_window(ii, ivec);
unsigned int n = vec_window.size();
for(unsigned int ii=0; ii<n; ii++) {
const double& value = (*vec_window[ii]).param(paramindex);
if(value!=IOUtils::nodata) {
sum += value;
count++;
}
}
if(count<=1) {
return IOUtils::nodata;
}
//compensated variance algorithm, see https://secure.wikimedia.org/wikipedia/en/wiki/Algorithms_for_calculating_variance
const double mean = sum/(double)count;
double sum2=0., sum3=0.;
for(unsigned int ii=0; ii<n; ii++) {
const double& value = (*vec_window[ii]).param(paramindex);
if(value!=IOUtils::nodata) {
sum2 = sum2 + (value - mean)*(value - mean);
sum3 = sum3 + (value - mean);
}
}
const double variance = (sum2 - sum3*sum3/count) / (count - 1);
return sqrt(variance);
}
double FilterTukey::getU3(const std::vector<MeteoData>& ivec, const unsigned int& ii, const unsigned int& paramindex) {