WSL/SLF GitLab Repository

Commit 18ec0b40 authored by Mathias Bavay's avatar Mathias Bavay
Browse files

Cleaner definition of epsilon in MeteoData. Finally, working 2D_interpolations...

Cleaner definition of epsilon in MeteoData. Finally, working 2D_interpolations test (with reference files and relying on numdiff).
parent b9c322ed
......@@ -24,12 +24,12 @@ int main(int /*argc*/, char** argv) {
Grid2DObject param;
io.interpolate(d1, dem, MeteoData::TA, param);
io.write2DGrid(param, MeteoGrids::TA, d1);
io.interpolate(d1, dem, MeteoData::HNW, param);
io.write2DGrid(param, MeteoGrids::HNW, d1);
//io.interpolate(d1, dem, MeteoData::HNW, param);
//io.write2DGrid(param, MeteoGrids::HNW, d1);
io.interpolate(d1, dem, MeteoData::RH, param);
io.write2DGrid(param, MeteoGrids::RH, d1);
io.interpolate(d1, dem, MeteoData::RSWR, param);
io.write2DGrid(param, MeteoGrids::RSWR, d1);
//io.interpolate(d1, dem, MeteoData::ILWR, param);
//io.write2DGrid(param, MeteoGrids::ILWR, d1);
return 0;
}
......@@ -3,8 +3,7 @@ CFLAGS = -Wall -Wextra
DEBUG = -g #-O3 # -DDEBUG -ggdb
METEOIODIR = ../../
LIBS_LINUX_EXTRA = -rdynamic -ldl
LIBS = -lstdc++ -L$(METEOIODIR)/lib -lmeteoio $(LIBS_LINUX_EXTRA)
LIBS = -rdynamic -lstdc++ -ldl -L../../lib -lmeteoio -ldl
INCLUDE=-I. -I$(METEOIODIR) -I$(METEOIODIR)/include
INCLUDE_POPC=-I. -I$(METEOIODIR) -I$(METEOIODIR)/include
......
......@@ -27,6 +27,8 @@ int main(int argc, char** argv) {
//Very basic conversion: get the whole data set at once, with its original sampling rate
//io.getMeteoData(d1, d2, vecMeteo);
Timer timer;
timer.start();
//More elaborate conversion: sample the data to a specific rate
//by looping over the time and calling readMeteoData for each timestep
std::vector<MeteoData> Meteo; //we need some intermediate storage, for storing data sets for 1 timestep
......@@ -39,10 +41,11 @@ int main(int argc, char** argv) {
}
}
timer.stop();
//In both case, we write the data out
std::cout << "Writing output data" << std::endl;
io.writeMeteoData(vecMeteo);
std::cout << "Done!!" << std::endl;
std::cout << "Done!! in " << timer.getElapsed() << " s" << std::endl;
return 0;
}
......@@ -5,43 +5,30 @@ using namespace mio; //The MeteoIO namespace is called mio
//This is a basic example of using as dem: the dem is read, the grid coordinates of a point given by its (lat,long) are retrieved
//and a sub-dem is extracted starting at these coordinates and extending dist_x and dist_y and written out.
int main(void) {
DEMObject dem;
//DEMObject dem;
Grid2DObject dem;
Config cfg("io.ini");
IOManager io(cfg);
//reading dem
dem.setUpdatePpt(DEMObject::SLOPE);
io.readDEM(dem);
//dem.setUpdatePpt(DEMObject::SLOPE);
//io.readDEM(dem);
io.read2DGrid(dem, "./input/surface-grids/Switzerland_1000m.asc");
//writing some statistics about this dem
//dem.grid2D.getMin() scans the DEM grid to get the min, while dem.min_altitude is cached and therefore very cheap
//The raw content of the 2D grids can also be accessed, for example dem.grid2D.getMin(IOUtils::RAW_NODATA). In this case, there would be no interpretation of some values as nodata.
std::cout << "DEM information: \n";
/*std::cout << "DEM information: \n";
std::cout << "\tmin=" << dem.grid2D.getMin() << " max=" << dem.grid2D.getMax() << " mean=" << dem.grid2D.getMean() << "\n";
std::cout << "\tmin slope=" << dem.min_slope << " max slope=" << dem.max_slope << std::endl;
std::cout << "\tmin slope=" << dem.min_slope << " max slope=" << dem.max_slope << std::endl;*/
Grid2DObject slope(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner, dem.slope);
/*Grid2DObject slope(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner, dem.slope);
io.write2DGrid(slope, MeteoGrids::SLOPE, Date(0.));
Grid2DObject azi(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner, dem.azi);
io.write2DGrid(azi, MeteoGrids::AZI, Date(0.));
io.write2DGrid(azi,"azi.asc");*/
//retrieving grid coordinates of a real world point
Coords point;
point.copyProj(dem.llcorner); //we use the same projection parameters as the DEM
point.setLatLon(46.232103, 7.362185, IOUtils::nodata);
dem.gridify(point);
//computing grid distances from real world distances
const double dist_x=70000., dist_y=120000.;
const unsigned int ncols = (unsigned int)ceil(dist_x/dem.cellsize);
const unsigned int nrows = (unsigned int)ceil(dist_y/dem.cellsize);
//extracting a sub-dem starting at the given coordinates and extending a given distance along x and along y
DEMObject sub_dem(dem, point.getGridI(), point.getGridJ(), ncols, nrows);
//writing the sub-dem out
io.write2DGrid(sub_dem, MeteoGrids::DEM, Date(0.));
io.write2DGrid(dem,MeteoGrids::DEM, Date(0.));
return 0;
}
[General]
PLUGINPATH = ../../lib/plugins
BUFF_CHUNK_SIZE = 30
BUFF_CHUNK_SIZE = 370
BUFF_BEFORE = 1.5
[Input]
......@@ -50,13 +50,14 @@ DEMFILE = ./input/surface-grids/Switzerland_1000m.asc
#SMET meteorological file format
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
STATION1 = *WFJ.smet
; 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
......@@ -84,15 +85,12 @@ STATION7 = TUJ3.smet
COORDSYS = CH1903
TIME_ZONE = 1
GRID2D = ARC
GRID2DPATH = ./
GRID2D = PNG
;png_legend = false
png_autoscale = false
;png_min_size = 400x400
png_max_size = 1366*768
; GRID2D = PNG
; GRID2DPATH = ./
; PNG_MIN_SIZE = 800x600
; PNG_WORLD_FILE = true
; PNG_AUTOSCALE = false
; PNG_SCALING = nearest
METEO = SMET
METEOPATH = ./
......
......@@ -634,27 +634,38 @@ size_t IOUtils::seek(const Date& soughtdate, const std::vector<MeteoData>& vecM,
{
//returns index of element, if element does not exist it returns closest index after soughtdate
//the element needs to be an exact hit or embedded between two measurments
if (vecM.size() <= 0) {//no elements in buffer
return IOUtils::npos;
}
size_t first = 0, last = vecM.size()-1;
const double start_val=vecM[first].date.getJulianDate(true);
const double end_val=vecM[last].date.getJulianDate(true);
const double curr_val = soughtdate.getJulianDate(true);
//if we reach this point: at least one element in buffer
if (vecM[0].date > soughtdate) {
if (vecM[first].date > soughtdate) {
return IOUtils::npos;
}
if (vecM[vecM.size()-1].date < soughtdate) {//last element is earlier, return npos
if (vecM[last].date < soughtdate) {//last element is earlier, return npos
return IOUtils::npos;
}
if (vecM[0].date == soughtdate) {//closest element
if (vecM[first].date == soughtdate) {//closest element
return 0;
}
const double raw_pos = (curr_val-start_val) / (end_val-start_val);
const size_t start = MAX( (size_t)(floor(raw_pos*last*.8)), first);
const size_t end = MIN( (size_t)ceil(raw_pos*last*1.2), last);
if(vecM[start].date.getJulianDate(true)<curr_val) first=start;
if(vecM[end].date.getJulianDate(true)>=curr_val) last=end;
//if we reach this point: the date is spanned by the buffer and there are at least two elements
if (exactmatch){
size_t first = 1, last = vecM.size()-1;
//first = 1; last = vecM.size()-1;
//size_t first = 1, last = vecM.size()-1;
//perform binary search
while (first <= last) {
......@@ -667,8 +678,7 @@ size_t IOUtils::seek(const Date& soughtdate, const std::vector<MeteoData>& vecM,
return mid; // found it. return position
}
} else {
size_t first = 0, last = vecM.size()-1;
//first = 0; last = vecM.size()-1;
//perform binary search
while (first <= last) {
const size_t mid = (first + last) / 2; // compute mid point
......
......@@ -71,7 +71,7 @@ const std::string& MeteoGrids::getParameterName(const size_t& parindex)
/************************************************************
* static section *
************************************************************/
const double MeteoData::epsilon = 1e-6;
const double MeteoData::epsilon = 1e-5;
const size_t MeteoData::nrOfParameters = MeteoData::lastparam - MeteoData::firstparam + 1;
map<size_t, string> MeteoData::static_meteoparamname;
std::vector<std::string> MeteoData::s_default_paramname;
......@@ -211,20 +211,18 @@ void MeteoData::setResampled(const bool& in_resampled)
bool MeteoData::operator==(const MeteoData& in) const
{
//An object is equal if the date is equal and all meteo parameters are equal
if (date != in.date){
if (date != in.date) {
return false;
}
if (nrOfAllParameters != in.nrOfAllParameters){ //the number of meteo parameters has to be consistent
if (nrOfAllParameters != in.nrOfAllParameters) { //the number of meteo parameters has to be consistent
return false;
}
for (size_t ii=0; ii<nrOfAllParameters; ii++) {
//const double epsilon = (fabs(data[ii]) < fabs(in.data[ii]) ? fabs(in.data[ii]) : fabs(data[ii])) * std::numeric_limits<double>::epsilon(); // Hack not working...
//const double epsilon = (fabs(data[ii]) < fabs(in.data[ii]) ? fabs(in.data[ii]) : fabs(data[ii])) * 0.0000001; // Hack not working with 0 == 0 ....
const double epsilon = 1.0e-5; // HACK bad is not relativ....
//const double epsilon_rel = (fabs(data[ii]) < fabs(in.data[ii]) ? fabs(in.data[ii]) : fabs(data[ii])) * std::numeric_limits<double>::epsilon(); // Hack not working...
//const double epsilon_rel = (fabs(data[ii]) < fabs(in.data[ii]) ? fabs(in.data[ii]) : fabs(data[ii])) * 0.0000001; // Hack not working with 0 == 0 ....
if( !IOUtils::checkEpsilonEquality(data[ii], in.data[ii], epsilon) ){
//cerr<<getParameterName(ii) << " : " <<data[ii] << " != " << in.data[ii] << "( epsilon : " << epsilon << ")" << endl; Hack still here for debug
return false;
}
}
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -4,12 +4,6 @@
using namespace mio; //The MeteoIO namespace is called mio
using namespace std;
// HACK SPEACK WITH Mathias if use ndiff-2 !!! (Problem no-data -999, relativ error, do it here ???)
//This is the a basic example of spatial interpolations. It does not check any exceptions, it only tries to be as c-like as possible
//provide date as ISO formatted, for example 2008-12-01T15:35:00 and
//it will retrieve and interpolate the data for this date according to the io.ini configuration file
int main(int /*argc*/, char** argv) {
Date d1;
......@@ -35,37 +29,5 @@ int main(int /*argc*/, char** argv) {
io.interpolate(d1, dem, MeteoData::RSWR, param);
io.write2DGrid(param, MeteoGrids::RSWR, d1);
string search_string = "-999.000";
string replace_string = "0";
string inbuf;
ifstream input_file("2009-01-19T12.00_HNW.asc");
ofstream output_file("first.txt");
while (!input_file.eof()){
size_t spot=0;
size_t hiho = 0;
getline(input_file, inbuf);
do {
spot = inbuf.find(search_string);
if(spot >= 0){
string tmpstring = inbuf.substr(0,spot);
tmpstring += replace_string;
tmpstring += inbuf.substr(spot+search_string.length(), inbuf.length());
inbuf = tmpstring;
}
cout << hiho++<< endl;
} while (spot != string::npos);
output_file << inbuf << endl;
}
//ifstream ifout("");
return 0;
}
#!/bin/sh
./2D_interpolations 2009-01-19T12:00
DATE="2009-01-19T12.00"
numdiff -r 1e-4 ${DATE}_HNW_ref.asc ${DATE}_HNW.asc
numdiff -r 1e-4 ${DATE}_RH_ref.asc ${DATE}_RH.asc
numdiff -r 1e-4 ${DATE}_RSWR_ref.asc ${DATE}_RSWR.asc
numdiff -r 1e-4 ${DATE}_TA_ref.asc ${DATE}_TA.asc
......@@ -5,9 +5,11 @@ ADD_EXECUTABLE(2D_interpolations 2D_interpolations.cc)
TARGET_LINK_LIBRARIES(2D_interpolations ${LIBRARIES})
# add the tests
ADD_TEST(2D_interpolations.smoke 2D_interpolations 2009-01-19T12:00)
#ADD_TEST(2D_interpolations.smoke 2D_interpolations 2009-01-19T12:00)
ADD_TEST(2D_interpolations.smoke 2D_interpolations.sh)
SET_TESTS_PROPERTIES(2D_interpolations.smoke
PROPERTIES LABELS smoke)
PROPERTIES LABELS smoke
FAIL_REGULAR_EXPRESSION "error")
......
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