WSL/SLF GitLab Repository

Commit fcd1c23d authored by Mathias Bavay's avatar Mathias Bavay
Browse files

Small improvements when computing a variogram (in order to be more tolerant to...

Small improvements when computing a variogram (in order to be more tolerant to nodata values), new key for virtual stations (to force using the whole DEM instead of a one point dem), a results cache has been created for the virtual stations (mixing virtual and real stations is a really, really bad idea), more flexibility for the Ryan and Liston wind interpolations.
parent 5574f95c
......@@ -25,10 +25,10 @@ namespace mio {
IOManager::IOManager(const std::string& filename_in) : cfg(filename_in), rawio(cfg), bufferedio(rawio, cfg),
meteoprocessor(cfg), interpolator(cfg), dataGenerator(cfg),
v_params(), v_coords(), v_stations(),
proc_properties(), point_cache(), filtered_cache(),
proc_properties(), virtual_point_cache(), 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),
virtual_stations(false), skip_virtual_stations(false)
virtual_stations(false), skip_virtual_stations(false), interpol_use_full_dem(false)
{
initIOManager();
}
......@@ -36,10 +36,10 @@ IOManager::IOManager(const std::string& filename_in) : cfg(filename_in), rawio(c
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(),
proc_properties(), virtual_point_cache(), 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),
virtual_stations(false), skip_virtual_stations(false)
virtual_stations(false), skip_virtual_stations(false), interpol_use_full_dem(false)
{
initIOManager();
}
......@@ -93,6 +93,8 @@ void IOManager::initVirtualStations()
v_stations.push_back( sd );
}
cfg.getValue("Interpol_Use_Full_DEM", "Input", interpol_use_full_dem, IOUtils::nothrow);
}
void IOManager::setProcessingLevel(const unsigned int& i_level)
......@@ -256,7 +258,7 @@ bool IOManager::read_filtered_cache(const Date& start_date, const Date& end_date
void IOManager::add_to_cache(const Date& i_date, const METEO_SET& vecMeteo)
{
//Check cache size, delete oldest elements if necessary
if (point_cache.size() > 2000){
if (point_cache.size() > 2000) {
point_cache.clear();
}
......@@ -332,35 +334,43 @@ size_t IOManager::getTrueMeteoData(const Date& i_date, METEO_SET& vecMeteo)
size_t IOManager::getVirtualMeteoData(const Date& i_date, METEO_SET& vecMeteo)
{
vecMeteo.clear();
// Check if data is available in cache
const map<Date, vector<MeteoData> >::const_iterator it = virtual_point_cache.find(i_date);
if (it != virtual_point_cache.end()){
vecMeteo = it->second;
return vecMeteo.size();
}
//get data from real input stations
METEO_SET vecTrueMeteo;
getTrueMeteoData(i_date, vecTrueMeteo);
if(vecTrueMeteo.empty()) return 0;
if (vecTrueMeteo.empty()) return 0;
if(v_params.empty()) {
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++) {
for (size_t ii=0; ii<vecStr.size(); ii++) {
v_params.push_back( vecTrueMeteo[0].getParameterIndex(vecStr[ii]) );
}
}
DEMObject dem;
bufferedio.readDEM(dem); //this is not a big deal since it will be in the buffer
//create stations without measurements
for(size_t ii=0; ii<v_stations.size(); ii++) {
for (size_t ii=0; ii<v_stations.size(); ii++) {
MeteoData md(i_date, v_stations[ii]);
vecMeteo.push_back( md );
}
//fill meteo parameters
DEMObject dem;
bufferedio.readDEM(dem); //this is not a big deal since it will be in the buffer
string info_string;
for(size_t param=0; param<v_params.size(); param++) {
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, info_string);
for(size_t ii=0; ii<v_coords.size(); ii++)
for (size_t ii=0; ii<v_coords.size(); ii++)
vecMeteo[ii](v_params[param]) = result[ii];
}
......@@ -372,13 +382,18 @@ size_t IOManager::getMeteoData(const Date& i_date, METEO_SET& vecMeteo)
{
vecMeteo.clear();
if(!virtual_stations || skip_virtual_stations)
if(!virtual_stations || skip_virtual_stations) {
getTrueMeteoData(i_date, vecMeteo);
else
add_to_cache(i_date, vecMeteo); //Store result in the local cache
} else {
getVirtualMeteoData(i_date, vecMeteo);
//Store result in the local cache
if (virtual_point_cache.size() > 2000) {
virtual_point_cache.clear();
}
virtual_point_cache[i_date] = vecMeteo;
}
//Store result in the local cache
add_to_cache(i_date, vecMeteo);
return vecMeteo.size();
}
......@@ -427,25 +442,44 @@ void IOManager::interpolate(const Date& date, const DEMObject& dem, const MeteoD
vector<Coords> vec_coords(in_coords);
for (size_t ii=0; ii<vec_coords.size(); ii++) {
const bool gridify_success = dem.gridify(vec_coords[ii]);
if (interpol_use_full_dem) {
Grid2DObject result_grid;
interpolator.interpolate(date, dem, meteoparam, result_grid, info_string);
const bool gridify_success = dem.gridify(vec_coords);
if (!gridify_success)
throw InvalidArgumentException("Coordinate given to interpolate is outside of dem", AT);
//Make new DEM with just one point, namely the one specified by vec_coord[ii]
//Copy all other properties of the big DEM into the new one
DEMObject one_point_dem(dem, (unsigned)vec_coords[ii].getGridI(), (unsigned)vec_coords[ii].getGridJ(), 1, 1, false);
one_point_dem.min_altitude = dem.min_altitude;
one_point_dem.max_altitude = dem.max_altitude;
one_point_dem.min_slope = dem.min_slope;
one_point_dem.max_slope = dem.max_slope;
one_point_dem.min_curvature = dem.min_curvature;
one_point_dem.max_curvature = dem.max_curvature;
Grid2DObject result_grid;
interpolator.interpolate(date, one_point_dem, meteoparam, result_grid, info_string);
result.push_back(result_grid(0,0));
for (size_t ii=0; ii<vec_coords.size(); ii++) {
//we know the i,j are positive because of gridify_success
const size_t pt_i = static_cast<size_t>( vec_coords[ii].getGridI() );
const size_t pt_j = static_cast<size_t>( vec_coords[ii].getGridJ() );
result.push_back( result_grid(pt_i,pt_j) );
}
} else {
for (size_t ii=0; ii<vec_coords.size(); ii++) {
const bool gridify_success = dem.gridify(vec_coords[ii]);
if (!gridify_success)
throw InvalidArgumentException("Coordinate given to interpolate is outside of dem", AT);
//we know the i,j are positive because of gridify_success
const size_t pt_i = static_cast<size_t>( vec_coords[ii].getGridI() );
const size_t pt_j = static_cast<size_t>( vec_coords[ii].getGridJ() );
//Make new DEM with just one point, namely the one specified by vec_coord[ii]
//Copy all other properties of the big DEM into the new one
DEMObject one_point_dem(dem, pt_i, pt_j, 1, 1, false);
one_point_dem.min_altitude = dem.min_altitude;
one_point_dem.max_altitude = dem.max_altitude;
one_point_dem.min_slope = dem.min_slope;
one_point_dem.max_slope = dem.max_slope;
one_point_dem.min_curvature = dem.min_curvature;
one_point_dem.max_curvature = dem.max_curvature;
Grid2DObject result_grid;
interpolator.interpolate(date, one_point_dem, meteoparam, result_grid, info_string);
result.push_back(result_grid(0,0));
}
}
skip_virtual_stations = false;
}
......
......@@ -234,12 +234,14 @@ class IOManager {
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 > virtual_point_cache; ///< stores already resampled virtual data points
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...
bool interpol_use_full_dem; ///< use full dem for point-wise spatial interpolations
};
} //end namespace
#endif
......@@ -240,19 +240,24 @@ void InterpolationAlgorithm::simpleWindInterpolate(const DEMObject& dem, const s
if (vecAltitudes.empty())
throw IOException("Not enough data for spatially interpolating wind", AT);
Fit1D trend;
getTrend(vecAltitudes, Ve, trend);
info << trend.getInfo();
detrend(trend, vecAltitudes, Ve);
Interpol2D::IDW(Ve, vecMeta, dem, VW);
retrend(dem, trend, VW);
getTrend(vecAltitudes, Vn, trend);
info << trend.getInfo();
detrend(trend, vecAltitudes, Vn);
Interpol2D::IDW(Vn, vecMeta, dem, DW);
retrend(dem, trend, DW);
if (vecDataVW.size()>=4) { //at least for points to perform detrending
Fit1D trend;
getTrend(vecAltitudes, Ve, trend);
info << trend.getInfo();
detrend(trend, vecAltitudes, Ve);
Interpol2D::IDW(Ve, vecMeta, dem, VW);
retrend(dem, trend, VW);
getTrend(vecAltitudes, Vn, trend);
info << trend.getInfo();
detrend(trend, vecAltitudes, Vn);
Interpol2D::IDW(Vn, vecMeta, dem, DW);
retrend(dem, trend, DW);
} else {
Interpol2D::IDW(Ve, vecMeta, dem, VW);
Interpol2D::IDW(Vn, vecMeta, dem, DW);
}
//recompute VW, DW in each cell
for (size_t ii=0; ii<VW.getNx()*VW.getNy(); ii++) {
......@@ -619,7 +624,6 @@ void ListonWindAlgorithm::calculate(const DEMObject& dem, Grid2DObject& grid)
Interpol2D::constant(0., dem, grid);
return;
}
if (param==MeteoData::VW) {
Grid2DObject DW;
simpleWindInterpolate(dem, vecDataVW, vecDataDW, grid, DW);
......@@ -999,30 +1003,45 @@ void OrdinaryKrigingAlgorithm::getDataForEmpiricalVariogram(std::vector<double>
}
//this gets the full data over a preceeding period
//we can not rely on the data in vecData/vecMeta since they filter nodata points
//and we need to do our own nodata handling here (to keep stations' indices constant)
void OrdinaryKrigingAlgorithm::getDataForVariogram(std::vector<double> &distData, std::vector<double> &variData, const bool& detrend_data)
{
distData.clear();
variData.clear();
vector<double> vecAltitudes;
getStationAltitudes(vecMeta, vecAltitudes);
//get all the data between "date" and "date-daysBefore"
const double daysBefore = 1.5;
const double Tstep = 1./24.;
Date d1 = date - daysBefore;
std::vector< std::vector<double> > vecVecData;
vecVecData.insert(vecVecData.begin(), nrOfMeasurments, std::vector<double>()); //allocation for the vectors
std::vector<MeteoData> Meteo;
iomanager.getMeteoData(d1, Meteo);
const size_t nrMeteo = Meteo.size();
vecVecData.insert(vecVecData.begin(), nrMeteo, std::vector<double>()); //allocation for the vectors
for(Date d1 = date - daysBefore; d1<=date; d1+=Tstep) {
std::vector<MeteoData> Meteo;
//get the stations altitudes
vector<double> vecAltitudes;
for (size_t ii=0; ii<nrMeteo; ii++){
const double& alt = Meteo[ii].meta.position.getAltitude();
if (alt != IOUtils::nodata) {
vecAltitudes.push_back(alt);
}
}
for(; d1<=date; d1+=Tstep) {
iomanager.getMeteoData(d1, Meteo);
if (Meteo.size()!=nrOfMeasurments)
throw InvalidArgumentException("Only a fixed number of stations is currently suppported for "+algo, AT);
if (Meteo.size()!=nrMeteo) {
std::ostringstream ss;
ss << "Number of stations varying between " << nrMeteo << " and " << Meteo.size();
ss << ". This is currently not supported for " << algo << "!";
throw InvalidArgumentException(ss.str(), AT);
}
if (detrend_data) {
//find trend
std::vector<double> vecDat;
for(size_t ii=0; ii<Meteo.size(); ii++)
for(size_t ii=0; ii<nrMeteo; ii++)
vecDat.push_back( Meteo[ii](param) );
Fit1D trend;
......@@ -1032,7 +1051,7 @@ void OrdinaryKrigingAlgorithm::getDataForVariogram(std::vector<double> &distData
throw InvalidArgumentException("Could not fit variogram model to the data", AT);
//detrend the data
for(size_t ii=0; ii<Meteo.size(); ii++) {
for(size_t ii=0; ii<nrMeteo; ii++) {
double val = Meteo[ii](param);
if(val!=IOUtils::nodata)
val -= trend( vecAltitudes[ii] );
......@@ -1040,7 +1059,7 @@ void OrdinaryKrigingAlgorithm::getDataForVariogram(std::vector<double> &distData
vecVecData.at(ii).push_back( val );
}
} else {
for(size_t ii=0; ii<Meteo.size(); ii++)
for(size_t ii=0; ii<nrMeteo; ii++)
vecVecData.at(ii).push_back( Meteo[ii](param) );
}
}
......@@ -1048,7 +1067,7 @@ void OrdinaryKrigingAlgorithm::getDataForVariogram(std::vector<double> &distData
//for each station, compute distance to other stations and
// variance of ( Y(current) - Y(other station) )
for(size_t j=0; j<nrOfMeasurments; j++) {
for(size_t j=0; j<nrMeteo; j++) {
const Coords& st1 = vecMeta[j].position;
const double x1 = st1.getEasting();
const double y1 = st1.getNorthing();
......
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