WSL/SLF GitLab Repository

Commit 7a61d513 authored by Mathias Bavay's avatar Mathias Bavay
Browse files

The estimation of the average sampling rate was quite wrong: it was not...

The estimation of the average sampling rate was quite wrong: it was not correct for multiple stations or for partially empty buffers (ie giving an end date way beyond the real end of the data). It was also claiming to return a value in 1/s but instead it was 1/days. This now only takes into account the time interval when there is some data and properly handles multiple stations and returns a value in Hz (or 1/s). The documentation has been expanded. Also some code cleanup (constification).
parent 8708ea04
......@@ -4,7 +4,7 @@
using namespace mio; //The MeteoIO namespace is called mio
//This is the most basic example. 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
//provide date as ISO formatted, for example 2008-12-01T15:35:00 and
//it will retrieve the data for this date according to the io.ini configuration file
int main(int /*argc*/, char** argv) {
Date d1;
......@@ -18,11 +18,12 @@ int main(int /*argc*/, char** argv) {
//io.setProcessingLevel(IOManager::raw); //set the processing level: raw, filtered or resampled
io.getMeteoData(d1, vecMeteo);
std::cout << vecMeteo.size() << " stations with an average sampling rate of " << io.getAvgSamplingRate() << "\n";
//writing some data out in order to prove that it really worked!
for (unsigned int ii=0; ii < vecMeteo.size(); ii++) {
std::cout << "---------- Station: " << (ii+1) << " / " << vecMeteo.size() << std::endl;
std::cout << vecMeteo[ii] << std::endl;
}
//for (unsigned int ii=0; ii < vecMeteo.size(); ii++) {
// std::cout << "---------- Station: " << (ii+1) << " / " << vecMeteo.size() << std::endl;
// std::cout << vecMeteo[ii] << std::endl;
//}
return 0;
}
......@@ -65,7 +65,7 @@ void BufferedIOHandler::bufferGrid(const Grid2DObject& in_grid2Dobj, const std::
void BufferedIOHandler::read2DGrid(Grid2DObject& in_grid2Dobj, const std::string& in_filename)
{
if(max_grids>0) {
std::map<std::string, Grid2DObject>::iterator it = mapBufferedGrids.find(in_filename);
const std::map<std::string, Grid2DObject>::iterator it = mapBufferedGrids.find(in_filename);
if (it != mapBufferedGrids.end()) { //already in map
in_grid2Dobj = (*it).second;
return;
......@@ -84,7 +84,7 @@ void BufferedIOHandler::read2DGrid(Grid2DObject& grid_out, const MeteoGrids::Par
{
if(max_grids>0) {
const string buffer_name = date.toString(Date::ISO)+"::"+MeteoGrids::getParameterName(parameter);
std::map<std::string, Grid2DObject>::iterator it = mapBufferedGrids.find(buffer_name);
const std::map<std::string, Grid2DObject>::iterator it = mapBufferedGrids.find(buffer_name);
if (it != mapBufferedGrids.end()) { //already in map
grid_out = (*it).second;
return;
......@@ -102,7 +102,7 @@ void BufferedIOHandler::read2DGrid(Grid2DObject& grid_out, const MeteoGrids::Par
void BufferedIOHandler::readDEM(DEMObject& in_grid2Dobj)
{
if(max_grids>0) {
std::map<std::string, Grid2DObject>::iterator it = mapBufferedGrids.find("/:DEM");
const std::map<std::string, Grid2DObject>::iterator it = mapBufferedGrids.find("/:DEM");
if (it != mapBufferedGrids.end()) {
//already in map. If the update properties have changed,
//we copy the ones given in input and force the update of the object
......@@ -130,7 +130,7 @@ void BufferedIOHandler::readDEM(DEMObject& in_grid2Dobj)
void BufferedIOHandler::readLanduse(Grid2DObject& in_grid2Dobj)
{
if(max_grids>0) {
std::map<std::string, Grid2DObject>::iterator it = mapBufferedGrids.find("/:LANDUSE");
const std::map<std::string, Grid2DObject>::iterator it = mapBufferedGrids.find("/:LANDUSE");
if (it != mapBufferedGrids.end()) { //already in map
in_grid2Dobj = (*it).second;
return;
......@@ -149,7 +149,7 @@ void BufferedIOHandler::readLanduse(Grid2DObject& in_grid2Dobj)
void BufferedIOHandler::readAssimilationData(const Date& in_date, Grid2DObject& in_grid2Dobj)
{
if(max_grids>0) {
std::map<std::string, Grid2DObject>::iterator it = mapBufferedGrids.find("/:ASSIMILATIONDATA" + in_date.toString(Date::FULL));
const std::map<std::string, Grid2DObject>::iterator it = mapBufferedGrids.find("/:ASSIMILATIONDATA" + in_date.toString(Date::FULL));
if (it != mapBufferedGrids.end()) { //already in map
in_grid2Dobj = (*it).second;
return;
......@@ -235,15 +235,21 @@ void BufferedIOHandler::setMinBufferRequirements(const double& i_chunk_size, con
double BufferedIOHandler::getAvgSamplingRate()
{
if (vec_buffer_meteo.size() > 0){
unsigned long sum = 0;
for (size_t ii=0; ii<vec_buffer_meteo.size(); ii++){
//count all data
sum += (unsigned long)vec_buffer_meteo[ii].size();
const size_t nr_stations = vec_buffer_meteo.size();
if (nr_stations > 0){
double sum = 0;
for (size_t ii=0; ii<nr_stations; ii++){ //loop over all stations
const size_t nr_data_pts = vec_buffer_meteo[ii].size();
if(nr_data_pts>1) {
const std::vector<MeteoData>& curr_station = vec_buffer_meteo[ii];
const double days = curr_station[nr_data_pts-1].date.getJulian() - curr_station[0].date.getJulian();
//add the average sampling rate for this station
if(days>0.) sum += (double)(nr_data_pts-1) / days; //the interval story: 2 points define 1 interval!
}
}
if (sum > 0){
double days = buffer_end.getJulian() - buffer_start.getJulian();
return ((double)sum / days);
if (sum > 0.){
return ((double)sum / ((double)nr_stations*24*3600)); //in points per seconds, ie Hz
}
}
......
......@@ -117,8 +117,12 @@ class BufferedIOHandler : public IOInterface {
/**
* @brief Returns the average sampling rate in the data.
* This computes the average sampling rate from the data that is contained in the buffer.
* @return average sampling rate in seconds, nodata if the buffer is empty
* This computes the average sampling rate of the data that is contained in the buffer. This is a quick
* estimate, centered on how often a station measures "something" (ie, how many timestamps do we have
* for this station in the buffer). if the station measures TA at h+0 and h+30 and
* RH at h+15 and h+45, it would return 4 measurements per hour. If the station measures TA and RH at h+0 and h+30,
* it would return 2 measurements per hour.
* @return average sampling rate in Hz, nodata if the buffer is empty
*/
double getAvgSamplingRate();
......
......@@ -922,17 +922,13 @@ double DEMObject::fillMissingGradient(const double& delta1, const double& delta2
void DEMObject::surfaceGradient(double& dx_sum, double& dy_sum, double A[4][4]) {
//Compute the gradient for a given cell (i,j) accross its eight surrounding cells (Horn, 1981)
double dx1, dx2, dx3;
double dy1, dy2, dy3;
double dx1 = lineGradient(A[3][1], A[3][2], A[3][3]);
double dx2 = lineGradient(A[2][1], A[2][2], A[2][3]);
double dx3 = lineGradient(A[1][1], A[1][2], A[1][3]);
//general calculation
dx1 = lineGradient(A[3][1], A[3][2], A[3][3]);
dx2 = lineGradient(A[2][1], A[2][2], A[2][3]);
dx3 = lineGradient(A[1][1], A[1][2], A[1][3]);
dy1 = lineGradient(A[3][1], A[2][1], A[1][1]);
dy2 = lineGradient(A[3][2], A[2][2], A[1][2]);
dy3 = lineGradient(A[3][3], A[2][3], A[1][3]);
double dy1 = lineGradient(A[3][1], A[2][1], A[1][1]);
double dy2 = lineGradient(A[3][2], A[2][2], A[1][2]);
double dy3 = lineGradient(A[3][3], A[2][3], A[1][3]);
//now trying to fill whatever could not be filled...
if(dx1==IOUtils::nodata) dx1 = fillMissingGradient(dx2, dx3);
......
......@@ -210,8 +210,12 @@ class IOManager {
/**
* @brief Returns the average sampling rate in the data.
* This computes the average sampling rate from the data that is contained in the buffer.
* @return average sampling rate in seconds, nodata if the buffer is empty
* This computes the average sampling rate of the data that is contained in the buffer. This is a quick
* estimate, centered on how often a station measures "something" (ie, how many timestamps do we have
* for this station in the buffer). if the station measures TA at h+0 and h+30 and
* RH at h+15 and h+45, it would return 4 measurements per hour. If the station measures TA and RH at h+0 and h+30,
* it would return 2 measurements per hour.
* @return average sampling rate in Hz, nodata if the buffer is empty
*/
double getAvgSamplingRate();
......
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