WSL/SLF GitLab Repository

Commit 105aad4b authored by Mathias Bavay's avatar Mathias Bavay
Browse files

The MOD operator has been added to Date (it takes a number of seconds in...

The MOD operator has been added to Date (it takes a number of seconds in argument). This is used by the new VSTATIONS_REFRESH_RATE key to avoid calling the spatial interpolations too often in the virtual stations. More checks have been implemented in the virtual stations and dowscaling for added safety. 

A major bug has been found when reading 3D grids in ARPS: the indexes when not in the right order.
parent c94c08e9
......@@ -26,7 +26,7 @@ Meteo2DInterpolator::Meteo2DInterpolator(const Config& i_cfg, TimeSeriesManager&
: cfg(i_cfg), tsmanager(i_tsmanager), gridsmanager(i_gridsmanager),
grid_buffer(0), mapAlgorithms(),
v_params(), v_coords(), v_stations(), virtual_point_cache(),
algorithms_ready(false), use_full_dem(false)
vstations_refresh_rate(IOUtils::unodata), algorithms_ready(false), use_full_dem(false)
{
size_t max_grids = 10; //default number of grids to keep in buffer
cfg.getValue("BUFF_GRIDS", "Interpolations2D", max_grids, IOUtils::nothrow);
......@@ -35,8 +35,15 @@ Meteo2DInterpolator::Meteo2DInterpolator(const Config& i_cfg, TimeSeriesManager&
setAlgorithms();
bool virtual_stations = false; ///< compute the meteo values at virtual stations
cfg.getValue("Virtual_stations", "Input", virtual_stations, IOUtils::nothrow);
if (virtual_stations) {
cfg.getValue("VSTATIONS_REFRESH_RATE", "Input", vstations_refresh_rate, IOUtils::nothrow);
}
bool downscaling = false; ///< Are we downscaling meteo grids instead of interpolating stations' data?
cfg.getValue("Downscaling", "Input", downscaling, IOUtils::nothrow);
if (virtual_stations && downscaling)
throw InvalidArgumentException("It is not possible to use both Virtual_stations and Downscaling!", AT);
if (virtual_stations || downscaling) {
initVirtualStations(downscaling); //adjust the coordinates if downscaling
cfg.getValue("Interpol_Use_Full_DEM", "Input", use_full_dem, IOUtils::nothrow);
......@@ -123,9 +130,8 @@ void Meteo2DInterpolator::interpolate(const Date& date, const DEMObject& dem, co
//Show algorithms to be used for this parameter
const map<string, vector<InterpolationAlgorithm*> >::iterator it = mapAlgorithms.find(param_name);
if (it==mapAlgorithms.end()) {
if (it==mapAlgorithms.end())
throw IOException("No interpolation algorithms configured for parameter "+param_name, AT);
}
//look for algorithm with the highest quality rating
const vector<InterpolationAlgorithm*>& vecAlgs = it->second;
......@@ -275,13 +281,14 @@ void Meteo2DInterpolator::check_projections(const DEMObject& dem, const std::vec
size_t Meteo2DInterpolator::getVirtualMeteoData(const vstations_policy& strategy, const Date& i_date, METEO_SET& vecMeteo)
{
if (strategy==VSTATIONS) {
//this reads station data, interpolates the stations and extract points from the interpolated grids
return getVirtualStationsData(i_date, vecMeteo);
} else if (strategy==SMART_DOWNSCALING) {
//This reads already gridded data and extract points from the grids
//Questions/tricks:
// throw exception if grids don't match between parameters
// should we recompute which points to take between each time steps? (ie more flexibility)
// should the generated virtual stations data be filtered? (useful for applying corrections)
// if so: have an own iomanager and use iomanager.push_meteo_data()
// move this to the IOHandler so the data could be filtered, corrected, interpolated
return getVirtualStationsFromGrid(i_date, vecMeteo);
} else if (strategy==DOWNSCALING) {
//extract all grid points
......@@ -299,9 +306,9 @@ void Meteo2DInterpolator::initVirtualStations(const bool& adjust_coordinates)
if (!cfg.keyExists("DEM", "Input"))
throw NoDataException("In order to use virtual stations, please provide a DEM!", AT);
DEMObject dem;
dem.setUpdatePpt( DEMObject::SLOPE ); //we only need the elevation
gridsmanager.readDEM(dem);
//get virtual stations coordinates
std::string coordin, coordinparam, coordout, coordoutparam;
IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);
......@@ -387,6 +394,7 @@ size_t Meteo2DInterpolator::getVirtualStationsData(const Date& i_date, METEO_SET
MeteoData md(i_date, v_stations[ii]);
vecMeteo.push_back( md );
}
if ( vstations_refresh_rate!=IOUtils::unodata && Date::mod(i_date, vstations_refresh_rate) != 0) return vecMeteo.size();
//fill meteo parameters
DEMObject dem;
......@@ -430,11 +438,18 @@ size_t Meteo2DInterpolator::getVirtualStationsFromGrid(const Date& i_date, METEO
vecMeteo.push_back( md );
}
DEMObject dem;
dem.setUpdatePpt( DEMObject::NO_UPDATE ); //we only need the elevation
gridsmanager.readDEM(dem); //this is not a big deal since it will be in the buffer
for (size_t param=0; param<v_params.size(); param++) { //loop over required parameters
const MeteoGrids::Parameters grid_param = static_cast<MeteoGrids::Parameters>(v_params[param]);
Grid2DObject grid;
gridsmanager.read2DGrid(grid, grid_param, i_date);
if (!grid.isSameGeolocalization(dem))
throw InvalidArgumentException("In SMART_DOWNSCALING, the DEM and the source grid don't match for '"+MeteoGrids::getParameterName(grid_param)+"' on "+i_date.toString(Date::ISO));
for (size_t ii=0; ii<v_stations.size(); ii++) { //loop over all virtual stations
const size_t grid_i = v_stations[ii].position.getGridI(); //this should work since invalid stations have been removed in init
const size_t grid_j = v_stations[ii].position.getGridJ();
......
......@@ -185,6 +185,8 @@ class Meteo2DInterpolator {
std::vector<Coords> v_coords; ///< Coordinates for virtual stations
std::vector<StationData> v_stations; ///< metadata for virtual stations
std::map<Date, METEO_SET > virtual_point_cache; ///< stores already resampled virtual data points
unsigned int vstations_refresh_rate; ///< how often to refresh the spatial interpolations for virtual stations? (in seconds)
bool algorithms_ready; ///< Have the algorithms objects been constructed?
bool use_full_dem; ///< use full dem for point-wise spatial interpolations
......
......@@ -761,6 +761,31 @@ bool Date::isLeapYear() const {
return (isLeapYear( getYear() ));
}
/**
* @brief Modulus of a julian date by a given number of seconds.
* This returns the modulus (in seconds) of a given date by the given number of seconds.
* @param[in] julian input date
* @param[in] seconds period (in seconds)
* @return modulus in seconds
*/
unsigned int Date::mod(const double& julian, const unsigned int& seconds)
{
const long int julian_rnd =static_cast<long int>( rnd(julian*3600.*24., 1) );
return static_cast<unsigned int>( julian_rnd % static_cast<long int>(seconds) );
}
/**
* @brief Modulus of a date by a given number of seconds.
* This returns the modulus (in seconds) of a given date by the given number of seconds.
* @param[in] julian input date
* @param[in] seconds period (in seconds)
* @return modulus in seconds
*/
unsigned int Date::mod(const Date& indate, const unsigned int& seconds)
{
return mod(indate.getJulian(), seconds);
}
/**
* @brief Round a julian date to a given precision.
* If you want to round a local date, do NOT provide it as gmt julian but as local julian,
......
......@@ -132,6 +132,8 @@ class Date {
int getJulianDayNumber(const bool& gmt=false) const;
bool isLeapYear() const;
static unsigned int mod(const double& julian, const unsigned int& seconds);
static unsigned int mod(const Date& indate, const unsigned int& seconds);
static double rnd(const double& julian, const unsigned int& precision, const RND& type=CLOSEST);
void rnd(const unsigned int& precision, const RND& type=CLOSEST);
static const Date rnd(const Date& indate, const unsigned int& precision, const RND& type=CLOSEST);
......
......@@ -301,9 +301,9 @@ void ARPSIO::read3DGrid(Grid3DObject& grid_out, const std::string& i_name)
moveToMarker(fin, filename, parameter);
//read the data we are interested in
for (size_t ix = 0; ix < dimx; ix++) {
for (size_t iz = 0; iz < dimz; iz++) {
for (size_t iy = 0; iy < dimy; iy++) {
for (size_t iz = 0; iz < dimz; iz++) {
for (size_t ix = 0; ix < dimx; ix++) {
double tmp;
if (fscanf(fin," %16lf%*[\n]",&tmp)==1) {
grid_out.grid3D(ix,iy,iz) = tmp;
......
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