WSL/SLF GitLab Repository

Commit 3b884d89 authored by Mathias Bavay's avatar Mathias Bavay
Browse files

The CH1903+ coordinate system has been implemented. in GSNIO, setting the...

The CH1903+ coordinate system has been implemented. in GSNIO, setting the station name was not always properly done. The ADD and MULT filters now offer a better parsing of their arguments in order to avoid silently failing; their documentation has been improved. When compiling with DATA_QA, parameters that are filtered but have no value (nodata) are now reported.
parent bb009879
......@@ -40,7 +40,7 @@ namespace mio {
* There are two ways of supporting a given coordinate system: through the use of an adhoc implementation
* (that becomes part of MeteoIO) or through the use of an external library, Proj4 [ref: http://trac.osgeo.org/proj/].
* The current internal implementations are the following (given by their keyword):
* - CH1903 for coordinates in the Swiss Grid [ref: http://geomatics.ladetto.ch/ch1903_wgs84_de.pdf]
* - CH1903 or CH1903+ for coordinates in the Swiss Grid [ref: http://geomatics.ladetto.ch/ch1903_wgs84_de.pdf]
* - UTM for UTM coordinates (the zone must be specified in the parameters, for example 31T) [ref: http://www.oc.nps.edu/oc2902w/maps/utmups.pdf]
* - UPS for Universal Polar Stereographic coordinates (the zone, either N or S, must be specified in the parameters). [ref: J. Hager, J. Behensky, B. Drew, <i>THE UNIVERSAL GRIDS: Universal Transverse Mercator (UTM) and Universal Polar Stereographic (UPS)</i>, 1989, Defense Mapping Agency, DMATM 8358.2]
* - LOCAL for local coordinate system (using the horizontal and vertical distance from a reference point, see Coords::geo_distances for the available choice of distance algorithms)
......@@ -745,6 +745,7 @@ void Coords::copyProj(const Coords& source, const bool i_update) {
*/
short int Coords::getEPSG() const {
if(coordsystem=="CH1903") return 21781;
if(coordsystem=="CH1903+") return 2056;
if(coordsystem=="UTM") {
//UTM Zone information
short int zoneNumber;
......@@ -802,6 +803,11 @@ void Coords::setEPSG(const int epsg) {
coord_param.clear();
found=true;
}
if(!found && (epsg==2056)) {
coord_sys.assign("CH1903+");
coord_param.clear();
found=true;
}
if(!found && (epsg>=32601) && (epsg<=32660)) {
//northern hemisphere
coord_sys.assign("UTM");
......@@ -849,6 +855,7 @@ void Coords::convert_to_WGS84(double i_easting, double i_northing, double& o_lat
if(coordsystem=="UTM") UTM_to_WGS84(i_easting, i_northing, o_latitude, o_longitude);
else if(coordsystem=="UPS") UPS_to_WGS84(i_easting, i_northing, o_latitude, o_longitude);
else if(coordsystem=="CH1903") CH1903_to_WGS84(i_easting, i_northing, o_latitude, o_longitude);
else if(coordsystem=="CH1903+") CH1903_to_WGS84(i_easting-2e6, i_northing-1e6, o_latitude, o_longitude);
else if(coordsystem=="LOCAL") local_to_WGS84(i_easting, i_northing, o_latitude, o_longitude);
else if(coordsystem=="PROJ4") PROJ4_to_WGS84(i_easting, i_northing, o_latitude, o_longitude);
else if(coordsystem=="NULL") NULL_to_WGS84(i_easting, i_northing, o_latitude, o_longitude);
......@@ -872,6 +879,11 @@ void Coords::convert_from_WGS84(double i_latitude, double i_longitude, double& o
if(coordsystem=="UTM") WGS84_to_UTM(i_latitude, i_longitude, o_easting, o_northing);
else if(coordsystem=="UPS") WGS84_to_UPS(i_latitude, i_longitude, o_easting, o_northing);
else if(coordsystem=="CH1903") WGS84_to_CH1903(i_latitude, i_longitude, o_easting, o_northing);
else if(coordsystem=="CH1903+") {
WGS84_to_CH1903(i_latitude, i_longitude, o_easting, o_northing);
o_easting += 2e6;
o_northing += 1e6;
}
else if(coordsystem=="LOCAL") WGS84_to_local(i_latitude, i_longitude, o_easting, o_northing);
else if(coordsystem=="PROJ4") WGS84_to_PROJ4(i_latitude, i_longitude, o_easting, o_northing);
else if(coordsystem=="NULL") WGS84_to_NULL(i_latitude, i_longitude, o_easting, o_northing);
......
......@@ -74,7 +74,8 @@ void ProcAdd::parse_args(const std::vector<std::string>& vec_args)
const size_t nrArgs = vec_args.size();
if (nrArgs==1) {
type='c'; //constant
IOUtils::convertString(offset, vec_args[0]);
if (!IOUtils::convertString(offset, vec_args[0]))
throw InvalidArgumentException("Invalid offset \""+vec_args[0]+"\" specified for the "+getName()+" filter. If correcting for a period, please specify the period!", AT);
} else if (nrArgs==2) {
const string type_str=IOUtils::strToUpper( vec_args[0] );
if (type_str=="MONTHLY") type='m';
......
......@@ -33,8 +33,8 @@ namespace mio {
* This adds to all values a given offset. Either a fixed value is given as single argument or a period
* (hourly/daily/monthly) as well as a filename (and absolute or relative path) containing the offsets to apply.
* This file must contain in the first column the indices (months from 1 to 12 or days from 1 to 366 or hours from 0 to 23)
* and the matching offset in the second column. Comments following the same syntax as in the ini file are accepted, missing
* indices are treated as 0.
* and the matching offset in the second column (<a href="http://www.cplusplus.com/reference/cctype/isspace/">whitespace</a> delimited).
* Comments following the same syntax as in the ini file are accepted, missing indices are treated as 0.
* @code
* TA::filter1 = add
* TA::arg1 = 2.5
......@@ -42,6 +42,20 @@ namespace mio {
* TSG::filter1 = add
* TSG::arg1 = daily input/TSG_corr.dat
* @endcode
*
* Example of correction file (monthly correction, December will receive a correction of 0):
* @code
* 01 -0.375
* 02 -1.932
* 03 -4.304
* 04 -2.449
* 05 -1.629
* 06 -1.734
* 07 -2.414
* 09 -1.289
* 10 -1.086
* 11 -0.769
* @endcode
*/
class ProcAdd : public ProcessingBlock {
......
......@@ -75,7 +75,8 @@ void ProcMult::parse_args(const std::vector<std::string>& vec_args)
const size_t nrArgs = vec_args.size();
if (nrArgs==1) {
type='c'; //constant
IOUtils::convertString(factor, vec_args[0]);
if (!IOUtils::convertString(factor, vec_args[0]))
throw InvalidArgumentException("Invalid factor \""+vec_args[0]+"\" specified for the "+getName()+" filter. If correcting for a period, please specify the period!", AT);
} else if (nrArgs==2) {
const string type_str=IOUtils::strToUpper( vec_args[0] );
if (type_str=="MONTHLY") type='m';
......
......@@ -33,14 +33,29 @@ namespace mio {
* This multiplies all values by a given factor. Either a fixed value is given as single argument or a period
* (hourly/daily/monthly) as well as a filename (and absolute or relative path) containing the factors to apply.
* This file must contain in the first column the indices (months from 1 to 12 or days from 1 to 366 or hours from 0 to 23)
* and the matching factor in the second column. Comments following the same syntax as in the ini file are accepted, missing
* indices are treated as 1.
* and the matching factor in the second column (<a href="http://www.cplusplus.com/reference/cctype/isspace/">whitespace</a> delimited).
* Comments following the same syntax as in the ini file are accepted, missing indices are treated as 1.
* @code
* HNW::filter1 = mult
* HNW::arg1 = 1.3
*
* ISWR::filter1 = mult
* ISWR::arg1 = montlhy input/ISWR_corr.dat
* ISWR::arg1 = monthly input/ISWR_corr.dat
* @endcode
*
* Example of correction file (monthly correction, August will receive a correction of 1):
* @code
* 01 0.440593
* 02 0.815111
* 03 0.475562
* 04 0.674975
* 05 0.700086
* 06 0.886783
* 07 1.70733
* 09 1.26533
* 10 0.577152
* 11 0.394095
* 12 0.347335
* @endcode
*/
......
......@@ -81,7 +81,7 @@ void ProcessingStack::process(const std::vector< std::vector<MeteoData> >& ivec,
ovec.resize( nr_stations );
for (size_t ii=0; ii<nr_stations; ii++){ //for every station
if( ivec[ii].empty() ) continue; //no data, nothing to do!
if ( ivec[ii].empty() ) continue; //no data, nothing to do!
//pick one element and check whether the param_name parameter exists
const size_t param = ivec[ii].front().getParameterIndex(param_name);
......@@ -89,6 +89,15 @@ void ProcessingStack::process(const std::vector< std::vector<MeteoData> >& ivec,
//since all filters start with ovec=ivec, maybe we could just swap pointers instead for copying
std::vector<MeteoData> tmp( ivec[ii] );
#ifdef DATA_QA
for (size_t kk=0; kk<ivec[ii].size(); kk++) {
const double orig = ivec[ii][kk](param);
if (orig==IOUtils::nodata) {
cout << "[DATA_QA] Filtering " << param_name << " is nodata " << ivec[ii][kk].date.toString(Date::ISO_TZ) << "\n";
}
}
#endif
//Now call the filters one after another for the current station and parameter
bool appliedFilter = false;
for (size_t jj=0; jj<nr_of_filters; jj++){
......@@ -104,13 +113,12 @@ void ProcessingStack::process(const std::vector< std::vector<MeteoData> >& ivec,
if (tmp.size() == ovec[ii].size()){
#ifdef DATA_QA
for(size_t kk=0; kk<ovec[ii].size(); kk++) {
for (size_t kk=0; kk<ovec[ii].size(); kk++) {
const double orig = tmp[kk](param);
const double filtered = ovec[ii][kk](param);
if(orig!=filtered) {
const string parname = tmp[kk].getNameForParameter(param);
if (orig!=filtered) {
const string filtername = (*filter_stack[jj]).getName();
cout << "[DATA_QA] Filtering " << parname << "::" << filtername << " " << tmp[kk].date.toString(Date::ISO_TZ) << "\n";
cout << "[DATA_QA] Filtering " << param_name << "::" << filtername << " " << tmp[kk].date.toString(Date::ISO_TZ) << "\n";
}
}
#endif
......@@ -141,7 +149,7 @@ const std::string ProcessingStack::toString() const
//os << "<ProcessingStack>";
os << setw(10) << param_name << "::";
for(size_t ii=0; ii<filter_stack.size(); ii++) {
for (size_t ii=0; ii<filter_stack.size(); ii++) {
os << setw(10) << (*filter_stack[ii]).toString();
}
......
......@@ -229,6 +229,8 @@ bool GSNIO::parseMetadata(std::stringstream& ss, StationData &sd, std::string &f
}
if (!line.compare(0, vsname_str.size(), vsname_str)) { //sensor name
name = line.substr(vsname_str.size());
IOUtils::trim(name);
if (valid == 15) { // Last station was valid: store StationData
buildStation(id, name, lat, lon, alt, slope_angle, slope_azi, sd);
ss.seekg(streampos, std::ios_base::beg); //point to the start of new station
......
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