WSL/SLF GitLab Repository

Commit 66541e39 authored by Mathias Bavay's avatar Mathias Bavay
Browse files

A new concept has been introduced: virtual stations. These are points where...

A new concept has been introduced: virtual stations. These are points where the meteo data will be transparently spatially interpolated and returned to getMeteoData calls. This allows a point model like Snowpack to work on virtual stations without any internal changes. A few new keys have been introduced in order to provide the necessary information: virtual stations location, meteo parameters to interpolate. For example:
VIRTUAL_STATIONS = true
VIRTUAL_PARAMETERS = TA RH HNW HS VW RSWR ILWR
VSTATION1 = 46.8 9.833333
VSTATION2 = 46.7 9.9

A few drawbacks still exist: the time range getMeteoData call is not supported yet. A DEM must be provided (ie DEM key in ini file). Obviously, all interpolated parameters must have been associated with interpolation algorithms. The virtual stations have IDs such as VIR1, VIR2, etc For performance reasons, the interpolation is then performed on a one cell dem, which means that algorithms using dem characteristics would (currently) be messed up. Basically, this is still a work in progress before everything is robust!

Otherwise, some POPC ifdef have been removed and the "iomanager.interpolate(... Grid2DObject)" call has been removed since the "iomanager.getMeteoData(... Grid2DObject)" call should be used instead (for consistency).
parent e9cfeddb
......@@ -22,11 +22,11 @@ int main(int /*argc*/, char** argv) {
//performing spatial interpolations
Grid2DObject param;
io.interpolate(d1, dem, MeteoData::TA, param);
io.getMeteoData(d1, dem, MeteoData::TA, param);
io.write2DGrid(param, MeteoGrids::TA, d1);
io.interpolate(d1, dem, MeteoData::HNW, param);
io.getMeteoData(d1, dem, MeteoData::HNW, param);
io.write2DGrid(param, MeteoGrids::HNW, d1);
io.interpolate(d1, dem, MeteoData::RH, param);
io.getMeteoData(d1, dem, MeteoData::RH, param);
io.write2DGrid(param, MeteoGrids::RH, d1);
return 0;
......
......@@ -338,6 +338,40 @@ Coords::Coords(const Coords& c) : ref_latitude(c.ref_latitude), ref_longitude(c.
grid_i(c.grid_i), grid_j(c.grid_j), grid_k(c.grid_k),
coordsystem(c.coordsystem), coordparam(c.coordparam), distance_algo(c.distance_algo) {}
/**
* @brief Local projection onstructor: this constructor is only suitable for building a local projection.
* Such a projection defines easting and northing as the distance (in meters) to a reference point
* which coordinates have to be provided here.
* @param[in] in_coordinatesystem string identifying the coordinate system to use
* @param[in] in_parameters string giving some additional parameters for the projection (empty string if not applicable)
* @param[in] coord_spec coordinate specification
*
* The coordinate specification is given as either: "easting northing epsg" or "lat lon".
*/
Coords::Coords(const std::string& in_coordinatesystem, const std::string& in_parameters, const std::string& coord_spec)
: ref_latitude(IOUtils::nodata), ref_longitude(IOUtils::nodata),
altitude(IOUtils::nodata), latitude(IOUtils::nodata), longitude(IOUtils::nodata),
easting(IOUtils::nodata), northing(IOUtils::nodata),
grid_i(IOUtils::inodata), grid_j(IOUtils::inodata), grid_k(IOUtils::inodata),
coordsystem(in_coordinatesystem), coordparam(in_parameters), distance_algo(GEO_COSINE)
{
std::istringstream iss(coord_spec);
double coord1=IOUtils::nodata, coord2=IOUtils::nodata;
int epsg=IOUtils::inodata;
iss >> std::skipws >> coord1;
iss >> std::skipws >> coord2;
iss >> std::skipws >> epsg;
if(coord1!=IOUtils::nodata && coord2!=IOUtils::nodata && epsg!=IOUtils::inodata) {
setEPSG(epsg);
setXY(coord1, coord2, IOUtils::nodata);
setProj(in_coordinatesystem, in_parameters);
} else if(coord1!=IOUtils::nodata && coord2!=IOUtils::nodata) {
setLatLon(coord1, coord2, IOUtils::nodata);
}
}
/**
* @brief Returns the East coordinate in the configured projection system
* @return easting
......
......@@ -57,6 +57,7 @@ class Coords {
//Constructors
Coords();
Coords(const std::string& in_coordinatesystem, const std::string& in_parameters="");
Coords(const std::string& in_coordinatesystem, const std::string& in_parameters, const std::string& coord_spec);
Coords(const double& in_lat_ref, const double& in_long_ref);
Coords(const Coords& c);
......
......@@ -24,12 +24,59 @@ namespace mio {
IOManager::IOManager(const Config& i_cfg) : cfg(i_cfg), rawio(cfg), bufferedio(rawio, cfg),
meteoprocessor(cfg), interpolator(cfg), dataGenerator(cfg),
v_params(), v_coords(), v_stations(),
proc_properties(), point_cache(), filtered_cache(),
fcache_start(Date(0.0, 0.)), fcache_end(Date(0.0, 0.)), //this should not matter, since 0 is still way back before any real data...
processing_level(IOManager::filtered | IOManager::resampled | IOManager::generated)
processing_level(IOManager::filtered | IOManager::resampled | IOManager::generated),
virtual_stations(false), skip_virtual_stations(false)
{
meteoprocessor.getWindowSize(proc_properties);
interpolator.setIOManager(*this); //because "*this" does not necessarily exist in the initialization list...
cfg.getValue("Virtual_stations", "Input", virtual_stations, IOUtils::nothrow);
if(virtual_stations) {
initVirtualStations();
}
}
void IOManager::initVirtualStations()
{
if(!cfg.keyExists("DEM", "Input"))
throw NoAvailableDataException("In order to use virtual stations, please provide a DEM!", AT);
DEMObject dem;
bufferedio.readDEM(dem);
//get virtual stations coordinates
std::string coordin, coordinparam, coordout, coordoutparam;
IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);
std::vector<std::string> vecStation;
cfg.getValues("Vstation", "INPUT", vecStation);
for(size_t ii=0; ii<vecStation.size(); ii++) {
Coords tmp(coordin, coordinparam, vecStation[ii]);
if(!tmp.isNodata())
v_coords.push_back( tmp );
}
//create stations' metadata
for(size_t ii=0; ii<v_coords.size(); ii++) {
if(!dem.gridify(v_coords[ii])) {
ostringstream ss;
ss << "Virtual station \"" << vecStation[ii] << "\" is not contained is provided DEM";
throw NoAvailableDataException(ss.str(), AT);
}
const size_t i = v_coords[ii].getGridI(), j = v_coords[ii].getGridJ();
v_coords[ii].setAltitude(dem(i,j), false);
ostringstream name;
name << "Virtual_Station_" << ii+1;
ostringstream id;
id << "VIR" << ii+1;
StationData sd(v_coords[ii], id.str(), name.str());
sd.setSlope(dem.slope(i,j), dem.azi(i,j));
v_stations.push_back( sd );
}
}
void IOManager::setProcessingLevel(const unsigned int& i_level)
......@@ -259,23 +306,59 @@ size_t IOManager::getTrueMeteoData(const Date& i_date, METEO_SET& vecMeteo)
return vecMeteo.size();
}
size_t IOManager::getVirtualMeteoData(const Date& i_date, METEO_SET& vecMeteo)
{
vecMeteo.clear();
METEO_SET vecTrueMeteo;
getTrueMeteoData(i_date, vecTrueMeteo);
if(vecTrueMeteo.empty()) return 0;
if(v_params.empty()) {
//get parameters to interpolate if not already done
//we need valid data in order to handle extra parameters
std::vector<std::string> vecStr;
cfg.getValue("Virtual_parameters", "Input", vecStr);
for(size_t ii=0; ii<vecStr.size(); ii++) {
v_params.push_back( vecTrueMeteo[0].getParameterIndex(vecStr[ii]) );
}
}
DEMObject dem;
bufferedio.readDEM(dem);
//create stations without measurements
for(size_t ii=0; ii<v_stations.size(); ii++) {
MeteoData md(i_date, v_stations[ii]);
vecMeteo.push_back( md );
}
//fill meteo parameters
for(size_t param=0; param<v_params.size(); param++) {
std::vector<double> result;
interpolate(i_date, dem, static_cast<MeteoData::Parameters>(v_params[param]), v_coords, result);
for(size_t ii=0; ii<v_coords.size(); ii++)
vecMeteo[ii](v_params[param]) = result[ii];
}
return vecMeteo.size();
}
//data can be raw or processed (filtered, resampled)
size_t IOManager::getMeteoData(const Date& i_date, METEO_SET& vecMeteo)
{
vecMeteo.clear();
getTrueMeteoData(i_date, vecMeteo);
if(!virtual_stations || skip_virtual_stations)
getTrueMeteoData(i_date, vecMeteo);
else
getVirtualMeteoData(i_date, vecMeteo);
//Store result in the local cache
add_to_cache(i_date, vecMeteo);
return vecMeteo.size();
}
#ifdef _POPC_ //HACK popc
void IOManager::writeMeteoData(/*const*/ std::vector< METEO_SET >& vecMeteo, /*const*/ std::string& name)
#else
void IOManager::writeMeteoData(const std::vector< METEO_SET >& vecMeteo, const std::string& name)
#endif
{
if (processing_level == IOManager::raw){
rawio.writeMeteoData(vecMeteo, name);
......@@ -284,65 +367,27 @@ void IOManager::writeMeteoData(const std::vector< METEO_SET >& vecMeteo, const s
}
}
#ifdef _POPC_ //HACK popc
bool IOManager::getMeteoData(/*const*/ Date& date, /*const*/ DEMObject& dem, /*const*/ MeteoData::Parameters& meteoparam,
Grid2DObject& result)
#else
bool IOManager::getMeteoData(const Date& date, const DEMObject& dem, const MeteoData::Parameters& meteoparam,
Grid2DObject& result)
#endif
{
string info_string;
interpolator.interpolate(date, dem, meteoparam, result, info_string);
const bool status = getMeteoData(date, dem, meteoparam, result, info_string);
cerr << "[i] Interpolating " << MeteoData::getParameterName(meteoparam);
cerr << " (" << info_string << ") " << endl;
return (!result.isEmpty());
return status;
}
#ifdef _POPC_ //HACK popc
bool IOManager::getMeteoData(/*const*/ Date& date, /*const*/ DEMObject& dem, /*const*/ MeteoData::Parameters& meteoparam,
Grid2DObject& result, std::string& info_string)
#else
bool IOManager::getMeteoData(const Date& date, const DEMObject& dem, const MeteoData::Parameters& meteoparam,
Grid2DObject& result, std::string& info_string)
#endif
{
skip_virtual_stations = true;
interpolator.interpolate(date, dem, meteoparam, result, info_string);
skip_virtual_stations = false;
return (!result.isEmpty());
}
#ifdef _POPC_ //HACK popc
void IOManager::interpolate(/*const*/ Date& date, /*const*/ DEMObject& dem, /*const*/ MeteoData::Parameters meteoparam,
Grid2DObject& result)
#else
void IOManager::interpolate(const Date& date, const DEMObject& dem, const MeteoData::Parameters& meteoparam,
Grid2DObject& result)
#endif
{
string info_string;
interpolate(date, dem, meteoparam, result, info_string);
cerr << "[i] Interpolating " << MeteoData::getParameterName(meteoparam);
cerr << " (" << info_string << ") " << endl;
}
#ifdef _POPC_ //HACK popc
void IOManager::interpolate(/*const*/ Date& date, /*const*/ DEMObject& dem, /*const*/ MeteoData::Parameters meteoparam,
Grid2DObject& result, std::string& info_string)
#else
void IOManager::interpolate(const Date& date, const DEMObject& dem, const MeteoData::Parameters& meteoparam,
Grid2DObject& result, std::string& info_string)
#endif
{
interpolator.interpolate(date, dem, meteoparam, result, info_string);
}
#ifdef _POPC_ //HACK popc
void IOManager::interpolate(/*const*/ Date& date, /*const*/ DEMObject& dem, /*const*/ MeteoData::Parameters& meteoparam,
/*const*/ std::vector<Coords>& in_coords, std::vector<double>& result)
#else
void IOManager::interpolate(const Date& date, const DEMObject& dem, const MeteoData::Parameters& meteoparam,
const std::vector<Coords>& in_coords, std::vector<double>& result)
#endif
{
string info_string;
interpolate(date, dem, meteoparam, in_coords, result, info_string);
......@@ -350,16 +395,11 @@ void IOManager::interpolate(const Date& date, const DEMObject& dem, const MeteoD
cerr << " (" << info_string << ") " << endl;
}
#ifdef _POPC_ //HACK popc
void IOManager::interpolate(/*const*/ Date& date, /*const*/ DEMObject& dem, /*const*/ MeteoData::Parameters& meteoparam,
/*const*/ std::vector<Coords>& in_coords, std::vector<double>& result,
std::string& info_string)
#else
void IOManager::interpolate(const Date& date, const DEMObject& dem, const MeteoData::Parameters& meteoparam,
const std::vector<Coords>& in_coords, std::vector<double>& result, std::string& info_string)
#endif
{
result.clear();
skip_virtual_stations = true;
vector<Coords> vec_coords = in_coords;
......@@ -389,6 +429,7 @@ void IOManager::interpolate(const Date& date, const DEMObject& dem, const MeteoD
result.push_back(result_grid.grid2D(0,0));
}
skip_virtual_stations = false;
}
void IOManager::read2DGrid(Grid2DObject& grid2D, const std::string& filename)
......
......@@ -23,6 +23,7 @@
#include <meteoio/BufferedIOHandler.h>
#include <meteoio/MeteoProcessor.h>
#include <meteoio/MeteoData.h>
#include <meteoio/Coords.h>
namespace mio {
......@@ -137,54 +138,18 @@ class IOManager {
* @param result grid returned filled with the requested data
* @return true if the grid got filled
*/
#ifdef _POPC_ //HACK popc
bool getMeteoData(/*const*/ Date& date, /*const*/ DEMObject& dem, /*const*/ MeteoData::Parameters& meteoparam,
Grid2DObject& result);
#else
bool getMeteoData(const Date& date, const DEMObject& dem, const MeteoData::Parameters& meteoparam,
Grid2DObject& result);
#endif
#ifdef _POPC_ //HACK popc
bool getMeteoData(/*const*/ Date& date, /*const*/ DEMObject& dem, /*const*/ MeteoData::Parameters& meteoparam,
Grid2DObject& result, std::string& info_string);
#else
bool getMeteoData(const Date& date, const DEMObject& dem, const MeteoData::Parameters& meteoparam,
Grid2DObject& result, std::string& info_string);
#endif
#ifdef _POPC_ //HACK popc
void interpolate(/*const*/ Date& date, /*const*/ DEMObject& dem, /*const*/ MeteoData::Parameters meteoparam,
Grid2DObject& result);
#else
void interpolate(const Date& date, const DEMObject& dem, const MeteoData::Parameters& meteoparam,
Grid2DObject& result);
#endif
#ifdef _POPC_ //HACK popc
void interpolate(/*const*/ Date& date, /*const*/ DEMObject& dem, /*const*/ MeteoData::Parameters meteoparam,
Grid2DObject& result, std::string& info_string);
#else
void interpolate(const Date& date, const DEMObject& dem, const MeteoData::Parameters& meteoparam,
bool getMeteoData(const Date& date, const DEMObject& dem, const MeteoData::Parameters& meteoparam,
Grid2DObject& result, std::string& info_string);
#endif
#ifdef _POPC_ //HACK popc
void interpolate(/*const*/ Date& date, /*const*/ DEMObject& dem, /*const*/ MeteoData::Parameters& meteoparam,
/*const*/ std::vector<Coords>& in_coords, std::vector<double>& result);
#else
void interpolate(const Date& date, const DEMObject& dem, const MeteoData::Parameters& meteoparam,
const std::vector<Coords>& in_coords, std::vector<double>& result);
#endif
#ifdef _POPC_ //HACK popc
void interpolate(/*const*/ Date& date, /*const*/ DEMObject& dem, /*const*/ MeteoData::Parameters& meteoparam,
/*const*/ std::vector<Coords>& in_coords, std::vector<double>& result,
std::string& info_string);
#else
void interpolate(const Date& date, const DEMObject& dem, const MeteoData::Parameters& meteoparam,
const std::vector<Coords>& in_coords, std::vector<double>& result, std::string& info_string);
#endif
/**
* @brief Set the desired ProcessingLevel of the IOManager instance
......@@ -221,11 +186,8 @@ class IOManager {
*/
double getAvgSamplingRate() const;
#ifdef _POPC_ //HACK popc
void writeMeteoData(/*const*/ std::vector< METEO_SET >& vecMeteo, /*const*/ std::string& name/*=""*/);
#else
void writeMeteoData(const std::vector< METEO_SET >& vecMeteo, const std::string& name="");
#endif
/**
* @brief Returns a copy of the internal Config object.
* This is convenient to clone an iomanager
......@@ -238,10 +200,12 @@ class IOManager {
void add_to_cache(const Date& i_date, const METEO_SET& vecMeteo);
private:
void initVirtualStations();
void fill_filtered_cache();
bool read_filtered_cache(const Date& start_date, const Date& end_date,
std::vector< METEO_SET >& vec_meteo);
size_t getTrueMeteoData(const Date& i_date, METEO_SET& vecMeteo);
size_t getVirtualMeteoData(const Date& i_date, METEO_SET& vecMeteo);
const Config cfg; ///< we keep this Config object as full copy, so the original one can get out of scope/be destroyed
IOHandler rawio;
......@@ -249,12 +213,18 @@ class IOManager {
MeteoProcessor meteoprocessor;
Meteo2DInterpolator interpolator;
DataGenerator dataGenerator;
ProcessingProperties proc_properties; ///< buffer constraints in order to be able to compute the requested values
std::vector<size_t> v_params; ///< Parameters for virtual stations
std::vector<Coords> v_coords; ///< Coordinates for virtual stations
std::vector<StationData> v_stations; ///< metadata for virtual stations
ProcessingProperties proc_properties; ///< buffer constraints in order to be able to compute the requested values
std::map<Date, METEO_SET > point_cache; ///< stores already resampled data points
std::vector< METEO_SET > filtered_cache; ///< stores already filtered data intervals
Date fcache_start, fcache_end; ///< store the beginning and the end date of the filtered_cache
unsigned int processing_level;
bool virtual_stations; ///< compute the meteo values at virtual stations
bool skip_virtual_stations; ///< skip virtual stations in subsequent calls to prevent recursive calls...
};
} //end namespace
#endif
......@@ -172,6 +172,10 @@ MeteoData::MeteoData(const Date& date_in)
: date(date_in), meta(), param_name(s_default_paramname), data(MeteoData::nrOfParameters, IOUtils::nodata), nrOfAllParameters(MeteoData::nrOfParameters), resampled(false)
{ }
MeteoData::MeteoData(const Date& date_in, const StationData& meta_in)
: date(date_in), meta(meta_in), param_name(s_default_paramname), data(MeteoData::nrOfParameters, IOUtils::nodata), nrOfAllParameters(MeteoData::nrOfParameters), resampled(false)
{ }
void MeteoData::setDate(const Date& in_date)
{
date = in_date;
......
......@@ -128,6 +128,13 @@ class MeteoData {
*/
MeteoData(const Date& in_date);
/**
* @brief A constructor that sets the measurment time and meta data
* @param in_date A Date object representing the time of the measurement
* @param meta_in A StationData object containing the meta data
*/
MeteoData(const Date& date_in, const StationData& meta_in);
/**
* @brief A setter function for the measurement date
* @param in_date A Date object representing the time of the measurement
......
......@@ -28,7 +28,7 @@ int main(int /*argc*/, char** argv) {
Grid2DObject param, ref;
int status = EXIT_SUCCESS;
io.interpolate(d1, dem, MeteoData::TA, param);
io.getMeteoData(d1, dem, MeteoData::TA, param);
if(gen_ref) io.write2DGrid(param, MeteoGrids::TA, d1);
else {
io.read2DGrid(ref, date_str+"_TA_ref.asc");
......@@ -37,7 +37,7 @@ int main(int /*argc*/, char** argv) {
}
}
io.interpolate(d1, dem, MeteoData::TSS, param);
io.getMeteoData(d1, dem, MeteoData::TSS, param);
if(gen_ref) io.write2DGrid(param, MeteoGrids::TSS, d1);
else {
io.read2DGrid(ref, date_str+"_TSS_ref.asc");
......@@ -46,7 +46,7 @@ int main(int /*argc*/, char** argv) {
}
}
io.interpolate(d1, dem, MeteoData::TSG, param);
io.getMeteoData(d1, dem, MeteoData::TSG, param);
if(gen_ref) io.write2DGrid(param, MeteoGrids::TSG, d1);
else {
io.read2DGrid(ref, date_str+"_TSG_ref.asc");
......@@ -55,7 +55,7 @@ int main(int /*argc*/, char** argv) {
}
}
io.interpolate(d1, dem, MeteoData::HNW, param);
io.getMeteoData(d1, dem, MeteoData::HNW, param);
if(gen_ref) io.write2DGrid(param, MeteoGrids::HNW, d1);
else {
io.read2DGrid(ref, date_str+"_HNW_ref.asc");
......@@ -64,7 +64,7 @@ int main(int /*argc*/, char** argv) {
}
}
io.interpolate(d1, dem, MeteoData::RH, param);
io.getMeteoData(d1, dem, MeteoData::RH, param);
if(gen_ref) io.write2DGrid(param, MeteoGrids::RH, d1);
else {
io.read2DGrid(ref, date_str+"_RH_ref.asc");
......@@ -73,7 +73,7 @@ int main(int /*argc*/, char** argv) {
}
}
io.interpolate(d1, dem, MeteoData::VW, param);
io.getMeteoData(d1, dem, MeteoData::VW, param);
if(gen_ref) io.write2DGrid(param, MeteoGrids::VW, d1);
else {
io.read2DGrid(ref, date_str+"_VW_ref.asc");
......@@ -82,7 +82,7 @@ int main(int /*argc*/, char** argv) {
}
}
io.interpolate(d1, dem, MeteoData::RSWR, param);
io.getMeteoData(d1, dem, MeteoData::RSWR, param);
if(gen_ref) io.write2DGrid(param, MeteoGrids::RSWR, d1);
else {
io.read2DGrid(ref, date_str+"_RSWR_ref.asc");
......@@ -91,7 +91,7 @@ int main(int /*argc*/, char** argv) {
}
}
io.interpolate(d1, dem, MeteoData::P, param);
io.getMeteoData(d1, dem, MeteoData::P, param);
if(gen_ref) io.write2DGrid(param, MeteoGrids::P, d1);
else {
io.read2DGrid(ref, date_str+"_P_ref.asc");
......
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