WSL/SLF GitLab Repository

Commit 328a1236 authored by Mathias Bavay's avatar Mathias Bavay
Browse files

A bug in the MeteoData::mergeTimeSeries() method has been fixed that was...

A bug in the MeteoData::mergeTimeSeries() method has been fixed that was related to the handling of extra parameters. The SHADE filter can now compute masks directly from the DEM, but this has required to propagate the Config object all the way down to the filter. Some minor code cleanup and documentation fixes.

Please notice that the SHADE filter is currently only appropriate for a single station! The most probable use case is to merge different stations into a new location (or to downscale to a new location) and properly mask the radiation by the surrounding terrain.
parent 859b9169
......@@ -70,10 +70,11 @@ class IOHandler : public IOInterface {
void create_exclude_map();
void create_keep_map();
void create_merge_map();
void checkTimestamps(const std::vector<METEO_SET>& vecVecMeteo) const;
void copy_parameters(const size_t& stationindex, std::vector< METEO_SET >& vecMeteo) const;
void exclude_params(std::vector<METEO_SET>& vecVecMeteo) const;
void keep_params(std::vector<METEO_SET>& vecVecMeteo) const;
void copy_parameters(const size_t& stationindex, std::vector< METEO_SET >& vecMeteo) const;
void merge_stations(std::vector<METEO_SET>& vecVecMeteo) const;
void merge_stations(STATIONS_SET& vecStation) const;
......
......@@ -204,8 +204,8 @@ size_t TimeSeriesManager::getMeteoData(const Date& i_date, METEO_SET& vecMeteo)
}
if ((IOUtils::resampled & processing_level) == IOUtils::resampled) { //resampling required
MeteoData md;
for (size_t ii=0; ii<(*data).size(); ii++) { //for every station
MeteoData md;
const bool success = meteoprocessor.resample(i_date, (*data)[ii], md);
if (success) vecMeteo.push_back(md);
}
......
......@@ -700,8 +700,8 @@ double DEMObject::getHorizon(const Coords& point, const double& bearing) const
void DEMObject::getHorizon(const Coords& point, const double& increment, std::vector< std::pair<double,double> >& horizon) const
{
for (double bearing=0.0; bearing <360.; bearing += increment) {
const double alpha = getHorizon(point, bearing);
horizon.push_back( make_pair(bearing, tan(alpha)*Cst::to_deg) );
const double tan_alpha = getHorizon(point, bearing);
horizon.push_back( make_pair(bearing, atan(tan_alpha)*Cst::to_deg) );
}
}
......
......@@ -253,7 +253,7 @@ bool MeteoData::operator!=(const MeteoData& in) const
double& MeteoData::operator()(const size_t& parindex)
{
#ifndef NOSAFECHECKS
if (parindex >= nrOfAllParameters)//getNrOfParameters())
if (parindex >= nrOfAllParameters)
throw IndexOutOfBoundsException("Trying to access meteo parameter that does not exist", AT);
#endif
return data[parindex];
......@@ -262,7 +262,7 @@ double& MeteoData::operator()(const size_t& parindex)
const double& MeteoData::operator()(const size_t& parindex) const
{
#ifndef NOSAFECHECKS
if (parindex >= nrOfAllParameters)//getNrOfParameters())
if (parindex >= nrOfAllParameters)
throw IndexOutOfBoundsException("Trying to access meteo parameter that does not exist", AT);
#endif
return data[parindex];
......@@ -365,13 +365,15 @@ void MeteoData::mergeTimeSeries(std::vector<MeteoData>& vec1, const std::vector<
{
if (vec1.empty() || vec2.empty()) return;
//if some new fields need to be created, do it on the front element
//so they can be found by other modules by testing the front element
//add any extra parameter as found in the first element of vec2
const size_t nrExtra2 = vec2.front().nrOfAllParameters - nrOfParameters;
for (size_t ii=0; ii<nrExtra2; ii++) {
const string extra_name = vec2.front().extra_param_name[ii];
if (vec1.front().getParameterIndex(extra_name)==IOUtils::npos)
vec1.front().addParameter( extra_name );
if (vec1.front().getParameterIndex(extra_name)==IOUtils::npos) {
for (size_t jj=0; jj<vec1.size(); jj++) {
vec1[jj].addParameter( extra_name );
}
}
}
//optimizations for special cases
......
......@@ -22,6 +22,7 @@
#include <algorithm>
#include <meteoio/meteoFilters/ProcShade.h>
#include <meteoio/IOHandler.h>
#include <meteoio/FileUtils.h>
using namespace std;
......@@ -39,8 +40,8 @@ const double ProcShade::soil_albedo = .23; //grass
const double ProcShade::snow_albedo = .85; //snow
const double ProcShade::snow_thresh = .1; //if snow height greater than this threshold -> snow albedo
ProcShade::ProcShade(const std::vector<std::string>& vec_args, const std::string& name, const std::string& i_root_path)
: ProcessingBlock(name), Suns(), mask(), root_path(i_root_path)
ProcShade::ProcShade(const std::vector<std::string>& vec_args, const std::string& name, const Config &i_cfg)
: ProcessingBlock(name), cfg(i_cfg), dem(), Suns(), mask()
{
parse_args(vec_args);
properties.stage = ProcessingProperties::first; //for the rest: default values
......@@ -62,6 +63,14 @@ void ProcShade::process(const unsigned int& param, const std::vector<MeteoData>&
it = Suns.find( stationHash );
}
//compute the mask from DEM if a DEM is available
if (mask.empty()) {
Coords position( ovec[0].meta.position );
if (!dem.gridify(position))
throw NoDataException("In filter '"+block_name+"', station '"+ovec[0].meta.stationID+"' is not included in the DEM", AT);
dem.getHorizon(position, 10., mask);
}
for (size_t ii=0; ii<ovec.size(); ii++){
double& tmp = ovec[ii](param);
if (tmp == IOUtils::nodata) continue; //preserve nodata values
......@@ -183,10 +192,11 @@ void ProcShade::readMask(const std::string& filter, const std::string& filename,
void ProcShade::parse_args(const std::vector<std::string>& vec_args)
{
const size_t nrArgs = vec_args.size();
if (nrArgs==0) {
//compute from DEM
throw InvalidArgumentException("Shading from DEM not implemented yet in filter " + getName(), AT);
if (nrArgs==0) { //compute from DEM
IOHandler io(cfg);
io.readDEM(dem);
} else if (nrArgs==1) {
const std::string root_path( cfg.getConfigRootDir() );
//if this is a relative path, prefix the path with the current path
const std::string in_filename = vec_args[0];
const std::string prefix = ( IOUtils::isAbsolutePath(in_filename) )? "" : root_path+"/";
......
......@@ -20,6 +20,8 @@
#include <meteoio/meteoFilters/FilterBlock.h>
#include <meteoio/meteoLaws/Sun.h>
#include <meteoio/Config.h>
#include <meteoio/dataClasses/DEMObject.h>
#include <vector>
#include <string>
......@@ -47,11 +49,14 @@ namespace mio {
* ISWR::arg1 = ../input/iswr_mask.dat
* @endcode
*
* If no arguments are provided, then it will compute the mask from the Digital Elevation Model. In such as case,
* a DEM must be declared in the [Input] section and must contain the stations of interest. Please make sure that the extend of the
* DEM is appropriate to correctly compute the shading effects!
*/
class ProcShade : public ProcessingBlock {
public:
ProcShade(const std::vector<std::string>& vec_args, const std::string& name, const std::string& i_root_path);
ProcShade(const std::vector<std::string>& vec_args, const std::string& name, const Config &i_cfg);
virtual void process(const unsigned int& param, const std::vector<MeteoData>& ivec,
std::vector<MeteoData>& ovec);
......@@ -62,9 +67,10 @@ class ProcShade : public ProcessingBlock {
void parse_args(const std::vector<std::string>& vec_args);
double getMaskElevation(const double& azimuth) const;
const Config &cfg;
DEMObject dem;
std::map<std::string, SunObject> Suns;
std::vector< std::pair<double,double> > mask;
std::string root_path;
static const double soil_albedo, snow_albedo, snow_thresh; ///< parametrize the albedo from HS
};
......
......@@ -117,7 +117,7 @@ namespace mio {
*
*/
ProcessingBlock* BlockFactory::getBlock(const std::string& blockname, const std::vector<std::string>& vec_args, const std::string& root_path)
ProcessingBlock* BlockFactory::getBlock(const std::string& blockname, const std::vector<std::string>& vec_args, const Config& cfg)
{
if (blockname == "SUPPR"){
return new FilterSuppr(vec_args, blockname);
......@@ -152,13 +152,13 @@ ProcessingBlock* BlockFactory::getBlock(const std::string& blockname, const std:
} else if (blockname == "UNVENTILATED_T"){
return new ProcUnventilatedT(vec_args, blockname);
} else if (blockname == "SHADE"){
return new ProcShade(vec_args, blockname, root_path);
return new ProcShade(vec_args, blockname, cfg);
} else if (blockname == "UNSHADE"){
return new ProcUnshade(vec_args, blockname);
} else if (blockname == "MULT"){
return new ProcMult(vec_args, blockname, root_path);
return new ProcMult(vec_args, blockname, cfg.getConfigRootDir());
} else if (blockname == "ADD"){
return new ProcAdd(vec_args, blockname, root_path);
return new ProcAdd(vec_args, blockname, cfg.getConfigRootDir());
} else if (blockname == "EXP_SMOOTHING"){
return new ProcExpSmoothing(vec_args, blockname);
} else if (blockname == "WMA_SMOOTHING"){
......
......@@ -19,6 +19,7 @@
#define PROCESSINGBLOCK_H
#include <meteoio/dataClasses/MeteoData.h>
#include <meteoio/Config.h>
#include <vector>
#include <string>
#include <set>
......@@ -86,7 +87,7 @@ class ProcessingBlock {
class BlockFactory {
public:
static ProcessingBlock* getBlock(const std::string& blockname, const std::vector<std::string>& vec_args, const std::string& root_path);
static ProcessingBlock* getBlock(const std::string& blockname, const std::vector<std::string>& vec_args, const Config& cfg);
};
} //end namespace
......
......@@ -25,8 +25,6 @@ namespace mio {
ProcessingStack::ProcessingStack(const Config& cfg, const std::string& parname) : filter_stack(), param_name(parname)
{
//this is required by filters that need to read some parameters in extra files
const string root_path = cfg.getConfigRootDir();
vector<string> vecFilters;
cfg.getValues(parname+"::filter", "Filters", vecFilters);
......@@ -40,7 +38,7 @@ ProcessingStack::ProcessingStack(const Config& cfg, const std::string& parname)
//Read arguments
cfg.getValue(tmp.str(), "Filters", vec_args, IOUtils::nothrow);
filter_stack.push_back( BlockFactory::getBlock(block_name, vec_args, root_path) );
filter_stack.push_back( BlockFactory::getBlock(block_name, vec_args, cfg) );
}
}
......
......@@ -280,7 +280,7 @@ void SunObject::getSlopeRadiation(const double& slope_azi, const double& slope_e
* D. G. Erbs, S.A. Klein, J.A. Duffie, <i>"Estimation of the diffuse radiation fraction for hourly, daily and monthly-average global radiation"</i>, Solar Energy, <b>28</b>, 4, 1982, Pages 293-302 and
* M. Iqbal, <i>"An introduction to solar radiation"</i>, 1983, Academic Press, ISBN: 0-12-373750-8 and
* D. T. Reindl, W. A. Beckman, J. A. Duffle, <i>"Diffuse fraction correlations</i>, Solar Energy, <b>45</b>, 1990, pp 1-7.
* @param iswr_modeled modelled radiation, it should be horizontal Top Of Atmosphere Radiation (W/m²)
* @param iswr_modeled modelled radiation, it should be horizontal Top Of %Atmosphere Radiation (W/m²)
* @param iswr_measured measured Incoming Short Wave Radiation on the ground (W/m²)
* @return splitting coefficient (between 0 and 1, 1 being 100% diffuse radiation)
*/
......@@ -325,7 +325,7 @@ double SunObject::getSplitting(const double& iswr_modeled, const double& iswr_me
* measured incoming global radiation to top of the atmosphere radiation toa_h.
* This is based on Boland, John, Lynne Scott, and Mark Luther, <i>"Modelling the diffuse fraction
* of global solar radiation on a horizontal surface"</i>, Environmetrics <b>12.2</b>, 2001, pp103-116.
* @param iswr_modeled modelled radiation, it should be horizontal Top Of Atmosphere Radiation (W/m²)
* @param iswr_modeled modelled radiation, it should be horizontal Top Of %Atmosphere Radiation (W/m²)
* @param iswr_measured measured Incoming Short Wave Radiation on the ground (W/m²)
* @param t solar time of day, ie solar time between 0 and 24
* @return splitting coefficient (between 0 and 1, 1 being 100% diffuse radiation)
......
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