WSL/SLF GitLab Repository

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

A new getMeteo() method has been implemented that transparently performs...

A new getMeteo() method has been implemented that transparently performs spatial interpolations and returns a grid for a given parameter. This should be prefered to interpolate(), since it allows the reading plugin to directly provide gridded data (such as coming out of GRIB).
Code cleanup in ResamplingAlgorithms, improved error message in A3DIO. Several keys have been updated/added to the plugins, in order to be more consistent between plugins. Therefore, meteo file extensions are given with METEOEXT, grid2d extensions with GRID2DEXT, grid2d prefix with GRID2DPREFIX (not all plugins support them, though), etc
parent eaaf7ff3
......@@ -22,7 +22,7 @@ using namespace std;
namespace mio {
IOManager::IOManager(const Config& i_cfg) : cfg(i_cfg), rawio(cfg), bufferedio(rawio, cfg),
IOManager::IOManager(const Config& i_cfg) : cfg(i_cfg), rawio(cfg), bufferedio(rawio, cfg),
meteoprocessor(cfg), interpolator(cfg, *this)
{
setProcessingLevel(IOManager::filtered | IOManager::resampled);
......@@ -255,6 +255,33 @@ void IOManager::writeMeteoData(const std::vector< METEO_TIMESERIE >& vecMeteo, c
}
}
#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);
cout << "[i] Interpolating " << MeteoData::getParameterName(meteoparam);
cout << " (" << info_string << ") " << endl;
return (!result.isEmpty());
}
#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
{
interpolator.interpolate(date, dem, meteoparam, result, info_string);
return (!result.isEmpty());
}
#ifdef _POPC_ //HACK popc
void IOManager::interpolate(/*const*/ Date& date, /*const*/ DEMObject& dem, /*const*/ MeteoData::Parameters meteoparam,
Grid2DObject& result)
......@@ -265,6 +292,8 @@ void IOManager::interpolate(const Date& date, const DEMObject& dem, const MeteoD
{
string info_string;
interpolate(date, dem, meteoparam, result, info_string);
cout << "[i] Interpolating " << MeteoData::getParameterName(meteoparam);
cout << " (" << info_string << ") " << endl;
}
#ifdef _POPC_ //HACK popc
......@@ -276,8 +305,6 @@ void IOManager::interpolate(const Date& date, const DEMObject& dem, const MeteoD
#endif
{
interpolator.interpolate(date, dem, meteoparam, result, info_string);
cout << "[i] Interpolating " << MeteoData::getParameterName(meteoparam);
cout << " (" << info_string << ") " << endl;
}
#ifdef _POPC_ //HACK popc
......@@ -290,6 +317,8 @@ void IOManager::interpolate(const Date& date, const DEMObject& dem, const MeteoD
{
string info_string;
interpolate(date, dem, meteoparam, in_coords, result, info_string);
cout << "[i] Interpolating " << MeteoData::getParameterName(meteoparam);
cout << " (" << info_string << ") " << endl;
}
#ifdef _POPC_ //HACK popc
......
......@@ -113,6 +113,44 @@ class IOManager {
void push_meteo_data(const ProcessingLevel& level, const Date& date_start, const Date& date_end,
const std::vector< METEO_TIMESERIE >& vecMeteo);
/**
* @brief Fill Grid2DObject with spatial data.
* Depending on which meteo plugin is in use, this might be spatially interpolated
* point measurements or grids as provided by the data source itself.
* Depending on the ProcessingLevel configured data will be either
* raw (read directly from the IOHandler)
*
* NOTE:
* - grid will be empty if there is no data found
*
* Example Usage:
* @code
* Grid2DObject grid; //empty grid
* IOManager iomanager(Config("io.ini"));
* iomanager.getMeteoData(Date(2008,06,21,11,00), MeteoData::TA, grid); //21.6.2008 11:00
* @endcode
* @param date A Date object representing the date/time for the sought MeteoData objects
* @param dem Digital Elevation Model data
* @param meteoparam which meteo parameter to return
* @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);
......
......@@ -17,6 +17,7 @@
*/
#include <meteoio/ResamplingAlgorithms.h>
#include <cmath>
#include <algorithm>
using namespace std;
......@@ -182,14 +183,14 @@ void ResamplingAlgorithms::LinearResampling(const size_t& pos, const size_t& par
if ((taskargs.size()==1) && (taskargs[0]=="extrapolate"))
extrapolate = true;
//Now find two points within the vecM (before and aft, that are not IOUtils::nodata)
//Now find two points within the vecM (before and after, that are not IOUtils::nodata)
//If that condition cannot be met, simply add nodata for the resampled value (exception: extrapolate)
size_t indexP1=IOUtils::npos, indexP2=IOUtils::npos;
bool foundP1=false, foundP2=false;
for (size_t ii=pos; (ii--) > 0; ){
if (vecM[ii](paramindex) != IOUtils::nodata){
indexP1 = (size_t)ii;
indexP1 = ii;
foundP1 = true;
break;
}
......@@ -197,7 +198,7 @@ void ResamplingAlgorithms::LinearResampling(const size_t& pos, const size_t& par
for (size_t ii=pos+1; ii<vecM.size(); ii++){
if (vecM[ii](paramindex) != IOUtils::nodata){
indexP2 = (size_t)ii;
indexP2 = ii;
foundP2 = true;
break;
}
......@@ -255,9 +256,7 @@ void ResamplingAlgorithms::LinearResampling(const size_t& pos, const size_t& par
void ResamplingAlgorithms::Accumulate(const size_t& pos, const size_t& paramindex,
const std::vector<std::string>& taskargs, std::vector<MeteoData>& vecM)
{
/*
* HACK TODO: Overall check IOUtils::nodata data path and test all scenarios with good test cases
*/
//HACK TODO: Overall check IOUtils::nodata data path and test all scenarios with good test cases
if (pos >= vecM.size())
throw IOException("The position of the resampled element is out of bounds", AT);
......@@ -279,10 +278,10 @@ void ResamplingAlgorithms::Accumulate(const size_t& pos, const size_t& paraminde
}
//find start of accumulation period
bool found_start=false;
size_t start_idx = pos+1;
Date dateStart(vecM[pos].date.getJulianDate() - accumulate_period/(24.*3600.), vecM[pos].date.getTimeZone());
bool found_start=false;
size_t start_idx;
for (start_idx=pos+1; (start_idx--) > 0; ){
if(vecM[start_idx].date <= dateStart) {
found_start=true;
......@@ -395,13 +394,13 @@ double ResamplingAlgorithms::linearInterpolation(const double& x1, const double&
const double& x2, const double& y2, const double& x3)
{
if (x1 == x2)
throw IOException("Attempted division by null", AT);
throw IOException("Attempted division by zero", AT);
//Solving y = kx +d
double k = (y1 - y2) / (x1 - x2);
double d = y2 - k*x2;
//Solving y = ax + b
const double a = (y2 - y1) / (x2 - x1);
const double b = y2 - a*x2;
return (k*x3 + d);
return (a*x3 + b);
}
} //namespace
......
......@@ -22,7 +22,6 @@
#include <meteoio/StationData.h>
#include <meteoio/meteostats/libinterpol1D.h>
#include <algorithm>
#include <iostream>
#include <string>
#include <vector>
......
......@@ -436,13 +436,16 @@ void A3DIO::read2DMeteo(std::vector< std::vector<MeteoData> >& vecMeteo)
constructMeteo2DFilenames(startDate, endDate, filenames);//get all files for all years
const size_t stations = getNrOfStations(filenames, hashStations);
constructMeteo2DFilenames(startDate, startDate, filenames);//get filenames for current year
std::cerr << "[I] Number of 2D meteo stations: " << stations << std::endl;
if (stations < 1) {
throw InvalidFormatException("[E] No StationData found in 2D Meteo Files", AT);
string tmp;
cfg.getValue("METEOPATH", "Input", tmp);
std::stringstream ss;
ss << "[E] No stations metadata found between between " << startDate.toString(Date::ISO) << " and " << endDate.toString(Date::ISO);
ss << " in header of 2D meteo files in " << tmp;
throw InvalidFormatException(ss.str(), AT);
}
constructMeteo2DFilenames(startDate, startDate, filenames);//get filenames for current year
std::vector<StationData> tmpvecS = std::vector<StationData>(stations); //stores unique stations
try {
......
......@@ -92,6 +92,7 @@ namespace mio {
* - COORDSYS: output coordinate system (see Coords) specified in the [Output] section
* - COORDPARAM: extra output coordinates parameters (see Coords) specified in the [Output] section
* - GRID2DPATH: meteo grids directory where to read/write the grids; [Input] and [Output] sections
* - GRID2DEXT: grid file extension, or <i>none</i> for no file extension (default: .asc)
* - A3D_VIEW: use Alpine3D's grid viewer naming scheme (default=false)? [Input] and [Output] sections.
* - DEMFILE: for reading the data as a DEMObject
* - LANDUSE: for interpreting the data as landuse codes
......@@ -138,6 +139,13 @@ void ARCIO::getGridPaths() {
cfg.getValue("GRID2D", "Output", tmp, Config::nothrow);
if (tmp == "ARC") //keep it synchronized with IOHandler.cc for plugin mapping!!
cfg.getValue("GRID2DPATH", "Output", grid2dpath_out);
grid2d_ext_in = ".asc";
cfg.getValue("GRID2DEXT", "Input", grid2d_ext_in, Config::nothrow);
if(grid2d_ext_in=="none") grid2d_ext_in="";
grid2d_ext_out = ".asc";
cfg.getValue("GRID2DEXT", "Output", grid2d_ext_out, Config::nothrow);
if(grid2d_ext_out=="none") grid2d_ext_out="";
}
ARCIO::~ARCIO() throw()
......@@ -259,7 +267,7 @@ void ARCIO::read2DGrid(Grid2DObject& grid_out, const MeteoGrids::Parameters& par
} else {
std::string date_str = date.toString(Date::ISO);
std::replace( date_str.begin(), date_str.end(), ':', '.');
read2DGrid_internal(grid_out, grid2dpath_in + "/" + date_str + "_" + MeteoGrids::getParameterName(parameter) + ".asc");
read2DGrid_internal(grid_out, grid2dpath_in + "/" + date_str + "_" + MeteoGrids::getParameterName(parameter) + grid2d_ext_in);
}
}
......@@ -372,11 +380,11 @@ void ARCIO::write2DGrid(const Grid2DObject& grid_in, const MeteoGrids::Parameter
write2DGrid(grid_in, date.toString(Date::NUM)+"."+ext );
} else {
if(parameter==MeteoGrids::DEM || parameter==MeteoGrids::AZI || parameter==MeteoGrids::SLOPE) {
write2DGrid(grid_in, MeteoGrids::getParameterName(parameter) + ".asc");
write2DGrid(grid_in, MeteoGrids::getParameterName(parameter) + grid2d_ext_out);
} else {
std::string date_str = date.toString(Date::ISO);
std::replace( date_str.begin(), date_str.end(), ':', '.');
write2DGrid(grid_in, date_str + "_" + MeteoGrids::getParameterName(parameter) + ".asc");
write2DGrid(grid_in, date_str + "_" + MeteoGrids::getParameterName(parameter) + grid2d_ext_out);
}
}
}
......
......@@ -77,6 +77,7 @@ class ARCIO : public IOInterface {
std::ofstream fout;//Output file streams
std::string coordin, coordinparam, coordout, coordoutparam; //projection parameters
std::string grid2dpath_in, grid2dpath_out;
std::string grid2d_ext_in, grid2d_ext_out; //file extension
};
} //end namespace mio
......
......@@ -27,20 +27,20 @@ namespace mio {
/**
* @page arps ARPSIO
* @section arps_format Format
* This is for reading grid data in the ARPS grid format (it transparently supports both true ARPS ascii grids and grids modified by the ARPSGRID utility). Currently, only DEM reading is implemented.
* This is for reading grid data in the ARPS grid format (it transparently supports both true ARPS ascii grids and grids modified by the ARPSGRID utility). DEM reading works well while reading meteo parameters might be a rough ride (since ARPS files do not always contain a consistent set of meteo fields).
*
* @section arps_units Units
*
* All units are assumed to be MKSA.
*
* @section arps_keywords Keywords
* This plugin uses the following keywords:
* - COORDSYS: coordinate system (see Coords); [Input] and [Output] section
* - COORDPARAM: extra coordinates parameters (see Coords); [Input] and [Output] section
* - DEMFILE: path and file containing the DEM; [Input] section
* - ARPS_X: x coordinate of the lower left corner of the grids; [Input] section
* - ARPS_Y: y coordinate of the lower left corner of the grids; [Input] section
* - GRID2DPATH: path to the input directory where to find the arps files to be read as grids; [Input] section //NOT USED YET
* - ARPS_EXT: arps file extension, or <i>none</i> for no file extension (default: .asc)
* - ARPS_XCOORD: x coordinate of the lower left corner of the grids; [Input] section
* - ARPS_YCOORD: y coordinate of the lower left corner of the grids; [Input] section
* - GRID2DPATH: path to the input directory where to find the arps files to be read as grids; [Input] section
* - GRID2DEXT: arps file extension, or <i>none</i> for no file extension (default: .asc)
*/
const double ARPSIO::plugin_nodata = -999.; //plugin specific nodata value
......@@ -86,11 +86,11 @@ void ARPSIO::setOptions()
cfg.getValue("GRID2DPATH", "Input", grid2dpath_in);
}
cfg.getValue("ARPS_X", "Input", xcoord, Config::dothrow);
cfg.getValue("ARPS_Y", "Input", ycoord, Config::dothrow);
cfg.getValue("ARPS_XCOORD", "Input", xcoord, Config::dothrow);
cfg.getValue("ARPS_YCOORD", "Input", ycoord, Config::dothrow);
ext = default_ext;
cfg.getValue("ARPS_EXT", "Input", ext, Config::nothrow);
cfg.getValue("GRID2DEXT", "Input", ext, Config::nothrow);
if(ext=="none") ext="";
}
......
......@@ -69,10 +69,12 @@ namespace mio {
* - COORDSYS: coordinate system (see Coords)
* - COORDPARAM: extra coordinates parameters (see Coords)
* - GRID2DPATH: path where to find the grids
* - GRID2DPREFIX: prefix to append when generating a file name for reading (ie: something like "laf" for Cosmo-Analysis-full domain), optional
* - GRID2DEXT: grib file extension, or <i>none</i> for no file extension (default: .grb)
* - GRIB_DEM_UPDATE: recompute slope/azimuth from the elevations when reading a DEM (default=false,
* that is we use the slope and azimuth included in the GRIB file)
* - GRIB_PREFIX: prefix to append when generating a file name for reading (ie: something like "laf" for Cosmo-Analysis-full domain), optional
* - GRIB_EXT: grib file extension, or <i>none</i> for no file extension (default: .grb)
* - METEOPATH: path where to find the grids for extracting time series at special points
* - METEOEXT: file extension, or <i>none</i> for no file extension (default: .grb)
* - STATION#: coordinates for virtual stations (if using GRIB as METEO plugin). Each station is given by its coordinates and the closest
* grid point will be chosen. Coordinates are given one one line as "lat lon" or "xcoord ycoord epsg_code". If a point leads to duplicate grid points,
* it will be removed from the list.
......@@ -138,11 +140,15 @@ void GRIBIO::setOptions()
cfg.getValue("GRID2DPATH", "Input", grid2dpath_in);
cfg.getValue("GRIB_DEM_UPDATE", "Input", update_dem, Config::nothrow);
}
cfg.getValue("GRIB_PREFIX", "Input", prefix, Config::nothrow);
cfg.getValue("GRID2DPREFIX", "Input", grid2d_prefix, Config::nothrow);
ext = default_ext;
cfg.getValue("GRIB_EXT", "Input", ext, Config::nothrow);
if(ext=="none") ext="";
meteo_ext = default_ext;
cfg.getValue("METEOEXT", "Input", meteo_ext, Config::nothrow);
if(meteo_ext=="none") meteo_ext="";
grid2d_ext = default_ext;
cfg.getValue("GRID2DEXT", "Input", grid2d_ext, Config::nothrow);
if(grid2d_ext=="none") grid2d_ext="";
}
void GRIBIO::readStations()
......@@ -507,7 +513,7 @@ void GRIBIO::read2DGrid(Grid2DObject& grid_out, const MeteoGrids::Parameters& pa
Date UTC_date = date;
UTC_date.setTimeZone(tz_in);
const std::string filename = grid2dpath_in+"/"+prefix+UTC_date.toString(Date::NUM).substr(0,10)+ext;
const std::string filename = grid2dpath_in+"/"+grid2d_prefix+UTC_date.toString(Date::NUM).substr(0,10)+grid2d_ext;
read2DGrid(filename, grid_out, parameter, UTC_date);
}
......@@ -724,7 +730,7 @@ void GRIBIO::readStationData(const Date&, std::vector<StationData>& /*vecStation
void GRIBIO::scanMeteoPath()
{
std::list<std::string> dirlist;
IOUtils::readDirectory(meteopath_in, dirlist, ext);
IOUtils::readDirectory(meteopath_in, dirlist, meteo_ext);
dirlist.sort();
//Check date in every filename and cache it
......
......@@ -90,7 +90,9 @@ class GRIBIO : public IOInterface {
std::string meteopath_in;
std::vector<Coords> vecPts; //points to use for virtual stations if METEO=GRIB
FILE *fp; //since passing fp always fail...
std::string ext; //file extension
std::string meteo_ext; //file extension
std::string grid2d_ext; //file extension
std::string grid2d_prefix; //filename prefix, like "laf"
bool indexed; //flag to know if the file has already been indexed
grib_index *idx; //because it needs to be kept between calls
std::string idx_filename; //matching file name for the index
......@@ -108,7 +110,6 @@ class GRIBIO : public IOInterface {
static const unsigned int MAX_VAL_LEN; //max value string lengthin GRIB
static const double plugin_nodata; //plugin specific nodata value, e.g. -999
static const double tz_in; //GRIB time zone
std::string prefix; //filename prefix, like "laf"
static const std::string default_ext;
std::string coordin, coordinparam; //projection parameters
bool update_dem;
......
......@@ -39,7 +39,7 @@ namespace mio {
* - COORDPARAM: extra input coordinates parameters (see Coords) specified in the [Input] section
* - COORDSYS: output coordinate system (see Coords) specified in the [Output] section
* - COORDPARAM: extra output coordinates parameters (see Coords) specified in the [Output] section
* - METAFILE: string containing the absolute filename of the geotop.inpts file in the [Input] section
* - METAFILE: string containing the absolute filename of the geotop.inpts file in the [Input] section
* - METEOPATH: string containing the path to the meteorological files
* - METEOPREFIX: file name prefix for meteorological files
* - METEOSEQ:specifiy in which order the columns should be printed out! ONLY relevant for writing out
......@@ -126,7 +126,7 @@ void GeotopIO::writeMeteoData(
map<string, size_t>::iterator it = mapParam.find(vecSequence[ii]);
if (it == mapParam.end())
throw InvalidFormatException("Key " + vecSequence[ii]
+ " invalid in io.ini:METEODESTSEQ", AT);
+ " invalid in io.ini:METEOSEQ", AT);
}
//write the meta data file _meteo.txt
......@@ -229,7 +229,7 @@ void GeotopIO::readMeteoData(const Date& dateStart, const Date& dateEnd,
if (vec_streampos.size() == 0) //the vec_streampos save file pointers for certain dates
vec_streampos = vector<map<Date, std::streampos> > (vecStation.size());
if (nr_of_stations == IOUtils::npos)
if (nr_of_stations == IOUtils::npos)
nr_of_stations = vecStation.size();
cout << "[i] GEOtopIO: Found " << nr_of_stations << " station(s)" << std::endl;
......@@ -303,7 +303,7 @@ void GeotopIO::readMeteoData(const Date& dateStart, const Date& dateEnd,
md(indices[jj-1]) = tmpdata[jj];
}
}
convertUnits(md);
vecMeteo[ii].push_back(md);
} else if (md.date > dateEnd) {
......@@ -464,7 +464,7 @@ void GeotopIO::readMetaData(const std::string& metafile) {
size_t pos = line.find("=");
string val = line.substr(pos+1);
IOUtils::trim(val);
if (!IOUtils::convertString(nr_of_stations, val, std::dec))
throw InvalidFormatException(metafile + ": " + line, AT);
}
......@@ -533,7 +533,7 @@ std::string GeotopIO::getValueForKey(const std::string& line)
size_t pos_start = line.find("\"");
size_t pos_end = line.find("\"", pos_start+1);
//cout << "t: pos_start = " << pos_start << " pos_end = " << pos_end << endl;
if ((pos_start != string::npos) && (pos_end != string::npos)) {
string param_name = line.substr(pos_start+1, pos_end - pos_start - 1);
IOUtils::trim(param_name);
......
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