WSL/SLF GitLab Repository

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

The K_TO_C and C_TO_K macros have been replaced by static methods in the...

The K_TO_C and C_TO_K macros have been replaced by static methods in the IOUtils namespace. They now properly check for nodata in order to avoid weird values and outputs...
parent 8b715791
......@@ -31,13 +31,6 @@
#include <meteoio/IOExceptions.h>
#include <meteoio/meteoLaws/Meteoconst.h>
#ifndef C_TO_K
#define C_TO_K( T ) ( T + Cst::t_water_freezing_pt ) // degree Celsius to kelvin
#endif
#ifndef K_TO_C
#define K_TO_C( T ) ( T - Cst::t_water_freezing_pt ) // kelvin to degree Celsius
#endif
namespace mio {
#ifdef _MSC_VER
......@@ -75,6 +68,12 @@ namespace IOUtils {
const double lon_epsilon = grid_epsilon / Cst::earth_R0 * Cst::to_deg; ///<in degrees. Small angle for longitudes, so sin(x)=x
const double lat_epsilon = lon_epsilon; ///<in degrees. Small angle for latitudes.
/*inline double C_TO_K(const double& T) {return ((T==nodata)? T : T + Cst::t_water_freezing_pt);}
inline double K_TO_C(const double& T) {return ((T==nodata)? T : T - Cst::t_water_freezing_pt);}*/
inline double C_TO_K(const double& T) {return (T + Cst::t_water_freezing_pt);}
inline double K_TO_C(const double& T) {return (T - Cst::t_water_freezing_pt);}
/**
* @brief Check whether two values are equal regarding a certain epsilon environment (within certain radius of each other)
* @param val1
......
......@@ -28,7 +28,7 @@ Meteo1DInterpolator::Meteo1DInterpolator(const Config& in_cfg)
//default window_size is 2 julian days
cfg.getValue("WINDOW_SIZE", "Interpolations1D", window_size, IOUtils::nothrow);
if (window_size <= 1.)
throw IOException("WINDOW_SIZE not valid", AT);
throw IOException("WINDOW_SIZE not valid, it should be a duration in seconds at least greater than 1", AT);
window_size /= 86400.; //user uses seconds, internally julian day is used
//read the Config object to create the resampling algorithms for each
......
......@@ -46,38 +46,38 @@ class MeteoGrids {
public:
/// \anchor meteogrids this enum provides names for possible meteogrids (from an ARPS file, etc)
enum Parameters {firstparam=0,
TA=firstparam, ///< Air temperature
RH, ///< Relative humidity
QI, ///< Specific humidity
TD, ///< Dew Point temperature
VW, ///< Wind velocity
DW, ///< Wind direction
VW_MAX, ///< Maximum wind velocity
ISWR, ///< Incoming short wave radiation
RSWR, ///< Reflected short wave radiation
ISWR_DIFF, ///< Incoming short wave, diffuse
ISWR_DIR, ///< Incoming short wave, direct
ILWR, ///< Incoming long wave radiation
TAU_CLD, ///< Cloud transmissivity or ISWR/ISWR_clear_sky
HS, ///< Height of snow
HNW, ///< Water equivalent of precipitations, either solid or liquid
HNW_S, ///< Water equivalent of precipitations, solid
HNW_L, ///< Water equivalent of precipitations, liquid
TSG, ///< Temperature ground surface
TSS, ///< Temperature snow surface
P, ///< Air pressure
P_SEA, ///< Sea level air pressure
U, ///< East component of wind
V, ///< North component of wind
W, ///< Vertical component of wind
SWE, ///< Snow Water Equivalent
ROT, ///< Total generated runoff
ALB, ///< Albedo
DEM, ///< Digital Elevation Model
SHADE, ///< Hillshade
SLOPE, ///< DEM slope angle
AZI, ///< DEM slope azimuth
lastparam=AZI};
TA=firstparam, ///< Air temperature
RH, ///< Relative humidity
QI, ///< Specific humidity
TD, ///< Dew Point temperature
VW, ///< Wind velocity
DW, ///< Wind direction
VW_MAX, ///< Maximum wind velocity
ISWR, ///< Incoming short wave radiation
RSWR, ///< Reflected short wave radiation
ISWR_DIFF, ///< Incoming short wave, diffuse
ISWR_DIR, ///< Incoming short wave, direct
ILWR, ///< Incoming long wave radiation
TAU_CLD, ///< Cloud transmissivity or ISWR/ISWR_clear_sky
HS, ///< Height of snow
HNW, ///< Water equivalent of precipitations, either solid or liquid
HNW_S, ///< Water equivalent of precipitations, solid
HNW_L, ///< Water equivalent of precipitations, liquid
TSG, ///< Temperature ground surface
TSS, ///< Temperature snow surface
P, ///< Air pressure
P_SEA, ///< Sea level air pressure
U, ///< East component of wind
V, ///< North component of wind
W, ///< Vertical component of wind
SWE, ///< Snow Water Equivalent
ROT, ///< Total generated runoff
ALB, ///< Albedo
DEM, ///< Digital Elevation Model
SHADE, ///< Hillshade
SLOPE, ///< DEM slope angle
AZI, ///< DEM slope azimuth
lastparam=AZI};
static const size_t nrOfParameters; ///<holds the number of meteo parameters stored in MeteoData
static const std::string& getParameterName(const size_t& parindex);
......
......@@ -74,7 +74,7 @@ void ProcUndercatch_Forland::process(const unsigned int& param, const std::vecto
//TA in celsius
double ProcUndercatch_Forland::solidPrecipitation(double TA, double VW)
{
TA = K_TO_C(TA); //convert to celsius
TA = IOUtils::K_TO_C(TA); //convert to celsius
if(type!=wfj)
VW = Atmosphere::windLogProfile(VW, 10., 2.); //impact seems minimal
else
......
......@@ -44,7 +44,7 @@ void ProcUndercatch_Hamon::process(const unsigned int& param, const std::vector<
VW = Atmosphere::windLogProfile(VW, 10., 2.); //impact seems minimal
double t = ovec[ii](MeteoData::TA);
if(t==IOUtils::nodata) continue; //we MUST have air temperature in order to filter
t = K_TO_C(t); //t in celsius
t = IOUtils::K_TO_C(t); //t in celsius
double k=0.;
if (tmp == IOUtils::nodata || tmp==0.) {
......
......@@ -46,7 +46,7 @@ void ProcUndercatch_WMO::process(const unsigned int& param, const std::vector<Me
if(VW!=IOUtils::nodata) VW = Atmosphere::windLogProfile(VW, 10., 2.); //impact seems minimal
double t = ovec[ii](MeteoData::TA);
if(t==IOUtils::nodata) continue; //we MUST have air temperature in order to filter
t=K_TO_C(t); //t in celsius
t=IOUtils::K_TO_C(t); //t in celsius
precip_type precip = (t<=Tsnow)? snow : (t>=Train)? rain : mixed;
//We don't use Tmax, Tmin, Tmean but only the current temperature instead
......@@ -88,7 +88,7 @@ void ProcUndercatch_WMO::process(const unsigned int& param, const std::vector<Me
const double alt = ovec[ii].meta.position.getAltitude();
double k=100.;
if(rh!=IOUtils::nodata && alt!=IOUtils::nodata) {
const double t_wb = K_TO_C(Atmosphere::wetBulbTemperature(ovec[ii](MeteoData::TA), rh, alt));
const double t_wb = IOUtils::K_TO_C(Atmosphere::wetBulbTemperature(ovec[ii](MeteoData::TA), rh, alt));
double ts_rate;
if(t_wb<1.1) ts_rate = 1. - .5*exp(-2.2*pow(1.1-t_wb, 1.3));
else ts_rate = .5*exp(-2.2*pow(t_wb-1.1, 1.3));
......@@ -154,9 +154,9 @@ void ProcUndercatch_WMO::parse_args(std::vector<std::string> filter_args)
IOUtils::convertString(factor_mixed, filter_args[2]);
if(filter_args.size()==5) {
IOUtils::convertString(Tsnow, filter_args[3]);
Tsnow = K_TO_C(Tsnow);
Tsnow = IOUtils::K_TO_C(Tsnow);
IOUtils::convertString(Train, filter_args[4]);
Train = K_TO_C(Train);
Train = IOUtils::K_TO_C(Train);
}
} else if(filter_args[0]=="nipher") {
type=nipher;
......
......@@ -186,12 +186,12 @@ double Atmosphere::windChill(const double& TA, const double& VW)
{
if(TA>(273.15+10.) || VW<5.) return TA; //not applicable in this case
const double t = K_TO_C(TA); //in Celsius
const double t = IOUtils::K_TO_C(TA); //in Celsius
const double v = VW*3.6; //in km/h
const double v016 = pow(v, 0.16);
const double WCT = 13.12 + 0.6215*t - 11.37*v016 + 0.3965*t*v016;
return C_TO_K(WCT);
return IOUtils::C_TO_K(WCT);
}
/**
......@@ -207,7 +207,7 @@ double Atmosphere::heatIndex(const double& TA, const double& RH)
{
if(TA<(273.15+26.7) || RH<0.4) return TA; //not applicable in this case
const double t = K_TO_C(TA); //in Celsius
const double t = IOUtils::K_TO_C(TA); //in Celsius
const double t2 = t*t;
const double rh = RH*100.; //in percent
const double rh2 = rh*rh;
......@@ -217,7 +217,7 @@ double Atmosphere::heatIndex(const double& TA, const double& RH)
const double HI = c1 + c2*t + c3*rh + c4*t*rh + c5*t2 + c6*rh2 + c7*t2*rh + c8*t*rh2 + c9*t2*rh2;
return C_TO_K(HI);
return IOUtils::C_TO_K(HI);
}
/**
......@@ -757,7 +757,7 @@ double Atmosphere::ILWR_parametrized(const double& lat, const double& lon, const
*/
double Atmosphere::RhtoDewPoint(double RH, double TA, const bool& force_water)
{
TA = K_TO_C(TA);
TA = IOUtils::K_TO_C(TA);
double Es, E, Tdw, Tdi; //saturation and current water vapor pressure
const double Aw = 611.21, Bw = 17.502, Cw = 240.97; //parameters for water
const double Ai = 611.15, Bi = 22.452, Ci = 272.55; //parameters for ice
......@@ -773,13 +773,13 @@ double Atmosphere::RhtoDewPoint(double RH, double TA, const bool& force_water)
Es = Aw * exp( (Bw * TA) / (Cw + TA) );
E = RH * Es;
Tdw = ( Cw * log(E / Aw) ) / ( Bw - log(E / Aw) );
return C_TO_K(Tdw);
return IOUtils::C_TO_K(Tdw);
}
if (TA < Tnucl) { //below nucleation, ice
Es = Ai * exp( (Bi * TA) / (Ci + TA) );
E = RH * Es;
Tdi = ( Ci * log(E / Ai) ) / ( Bi - log(E / Ai) );
return C_TO_K(Tdi);
return IOUtils::C_TO_K(Tdi);
}
//no clear state, we do a smooth interpolation between water and ice
......@@ -791,7 +791,7 @@ double Atmosphere::RhtoDewPoint(double RH, double TA, const bool& force_water)
E = RH * Es;
Tdw = ( Cw * log(E / Aw) ) / ( Bw - log(E / Aw) );
return C_TO_K( (di / (di + dw) * Tdi + dw / (di + dw) * Tdw) );
return IOUtils::C_TO_K( (di / (di + dw) * Tdi + dw / (di + dw) * Tdw) );
}
/**
......@@ -805,8 +805,8 @@ double Atmosphere::DewPointtoRh(double TD, double TA, const bool& force_water)
{
//Convert a dew point temperature into a Relative Humidity
//TA, TD are in Kelvins, RH is returned between 0 and 1
TA = K_TO_C(TA);
TD = K_TO_C(TD);
TA = IOUtils::K_TO_C(TA);
TD = IOUtils::K_TO_C(TD);
double Es, E, Rhi, Rhw, Rh; //saturation and current water vapro pressure
const double Aw = 611.21, Bw = 17.502, Cw = 240.97; //parameters for water
const double Ai = 611.15, Bi = 22.452, Ci = 272.55; //parameters for ice
......
......@@ -180,12 +180,10 @@ void A3DIO::convertUnits(MeteoData& meteo)
//converts C to Kelvin, converts RH to [0,1]
double& ta = meteo(MeteoData::TA);
if (ta != IOUtils::nodata)
ta = C_TO_K(ta);
ta = IOUtils::C_TO_K(ta);
double& tsg = meteo(MeteoData::TSG);
if (tsg != IOUtils::nodata)
tsg = C_TO_K(tsg);
tsg = IOUtils::C_TO_K(tsg);
double& rh = meteo(MeteoData::RH);
if (rh != IOUtils::nodata)
......@@ -755,7 +753,7 @@ int A3DIO::create1DFile(const std::vector< std::vector<MeteoData> >& data)
if(data[ii][j](MeteoData::TA) == IOUtils::nodata)
file << setw(6) << setprecision(0) << IOUtils::nodata << " ";
else
file << setw(6) << setprecision(2) << K_TO_C(data[ii][j](MeteoData::TA)) << " ";
file << setw(6) << setprecision(2) << IOUtils::K_TO_C(data[ii][j](MeteoData::TA)) << " ";
if(data[ii][j](MeteoData::ISWR) == IOUtils::nodata)
file << setw(6) << setprecision(0) << IOUtils::nodata << " ";
else
......@@ -866,7 +864,7 @@ int A3DIO::write2DmeteoFile(const std::vector< std::vector<MeteoData> >& data,
if(value==IOUtils::nodata) {
file << " " << setw(7) << setprecision(0) << IOUtils::nodata;
} else {
if(parindex==mio::MeteoData::TA) value = K_TO_C(value);
if(parindex==mio::MeteoData::TA) value = IOUtils::K_TO_C(value);
if(parindex==mio::MeteoData::RH) value = value*100.;
file << " " << setw(7) << setprecision(2) << value;
}
......
......@@ -431,16 +431,13 @@ void BormaIO::convertUnits(MeteoData& meteo)
//converts C to Kelvin, converts RH to [0,1]
double& ta = meteo(MeteoData::TA);
if (ta != IOUtils::nodata)
ta = C_TO_K(ta);
ta = IOUtils::C_TO_K(ta);
double& tsg = meteo(MeteoData::TSG);
if (tsg != IOUtils::nodata)
tsg = C_TO_K(tsg);
tsg = IOUtils::C_TO_K(tsg);
double& tss = meteo(MeteoData::TSS);
if (tss != IOUtils::nodata)
tss = C_TO_K(tss);
tss = IOUtils::C_TO_K(tss);
double& rh = meteo(MeteoData::RH);
if (rh != IOUtils::nodata)
......
......@@ -454,16 +454,13 @@ void GSNIO::convertUnits(MeteoData& meteo) const
{
//converts C to Kelvin, converts RH to [0,1]
double& ta = meteo(MeteoData::TA);
if (ta != IOUtils::nodata)
ta = C_TO_K(ta);
ta = IOUtils::C_TO_K(ta);
double& tsg = meteo(MeteoData::TSG);
if (tsg != IOUtils::nodata)
tsg = C_TO_K(tsg);
tsg = IOUtils::C_TO_K(tsg);
double& tss = meteo(MeteoData::TSS);
if (tss != IOUtils::nodata)
tss = C_TO_K(tss);
tss = IOUtils::C_TO_K(tss);
double& rh = meteo(MeteoData::RH);
if (rh != IOUtils::nodata)
......
......@@ -568,16 +568,13 @@ void GeotopIO::write2DGrid(const Grid2DObject&, const MeteoGrids::Parameters&, c
void GeotopIO::convertUnitsBack(MeteoData& meteo) {
//converts Kelvin to C, converts RH to [0,100]
double& ta = meteo(MeteoData::TA);
if (ta != IOUtils::nodata)
ta = K_TO_C(ta);
ta = IOUtils::K_TO_C(ta);
double& tsg = meteo(MeteoData::TSG);
if (tsg != IOUtils::nodata)
tsg = K_TO_C(tsg);
tsg = IOUtils::K_TO_C(tsg);
double& tss = meteo(MeteoData::TSS);
if (tss != IOUtils::nodata)
tss = K_TO_C(tss);
tss = IOUtils::K_TO_C(tss);
double& rh = meteo(MeteoData::RH);
if (rh != IOUtils::nodata)
......@@ -593,16 +590,13 @@ void GeotopIO::convertUnits(MeteoData& meteo) {
//converts C to Kelvin, converts RH to [0,1]
double& ta = meteo(MeteoData::TA);
if (ta != IOUtils::nodata)
ta = C_TO_K(ta);
ta = IOUtils::C_TO_K(ta);
double& tsg = meteo(MeteoData::TSG);
if (tsg != IOUtils::nodata)
tsg = C_TO_K(tsg);
tsg = IOUtils::C_TO_K(tsg);
double& tss = meteo(MeteoData::TSS);
if (tss != IOUtils::nodata)
tss = C_TO_K(tss);
tss = IOUtils::C_TO_K(tss);
double& rh = meteo(MeteoData::RH);
if (rh != IOUtils::nodata)
......
......@@ -1056,16 +1056,13 @@ void ImisIO::convertUnits(MeteoData& meteo)
//converts C to Kelvin, converts RH to [0,1]
double& ta = meteo(MeteoData::TA);
if (ta != IOUtils::nodata)
ta = C_TO_K(ta);
ta = IOUtils::C_TO_K(ta);
double& tsg = meteo(MeteoData::TSG);
if (tsg != IOUtils::nodata)
tsg = C_TO_K(tsg);
tsg = IOUtils::C_TO_K(tsg);
double& tss = meteo(MeteoData::TSS);
if (tss != IOUtils::nodata)
tss = C_TO_K(tss);
tss = IOUtils::C_TO_K(tss);
double& rh = meteo(MeteoData::RH);
if (rh != IOUtils::nodata)
......
......@@ -803,16 +803,13 @@ void PSQLIO::convertUnitsBack(MeteoData& meteo)
{
//converts Kelvin to °C, converts RH to [0,100]
double& ta = meteo(MeteoData::TA);
if (ta != IOUtils::nodata)
ta = K_TO_C(ta);
ta = IOUtils::K_TO_C(ta);
double& tsg = meteo(MeteoData::TSG);
if (tsg != IOUtils::nodata)
tsg = K_TO_C(tsg);
tsg = IOUtils::K_TO_C(tsg);
double& tss = meteo(MeteoData::TSS);
if (tss != IOUtils::nodata)
tss = K_TO_C(tss);
tss = IOUtils::K_TO_C(tss);
double& rh = meteo(MeteoData::RH);
if (rh != IOUtils::nodata)
......@@ -831,16 +828,13 @@ void PSQLIO::convertUnits(MeteoData& meteo)
{
//converts °C to Kelvin, converts RH to [0,1]
double& ta = meteo(MeteoData::TA);
if (ta != IOUtils::nodata)
ta = C_TO_K(ta);
ta = IOUtils::C_TO_K(ta);
double& tsg = meteo(MeteoData::TSG);
if (tsg != IOUtils::nodata)
tsg = C_TO_K(tsg);
tsg = IOUtils::C_TO_K(tsg);
double& tss = meteo(MeteoData::TSS);
if (tss != IOUtils::nodata)
tss = C_TO_K(tss);
tss = IOUtils::C_TO_K(tss);
double& rh = meteo(MeteoData::RH);
if (rh != IOUtils::nodata)
......
......@@ -281,12 +281,10 @@ void SASEIO::convertUnits(MeteoData& meteo)
//converts C to Kelvin, converts RH to [0,1]
double& ta = meteo(MeteoData::TA);
if (ta != IOUtils::nodata)
ta = C_TO_K(ta);
ta = IOUtils::C_TO_K(ta);
double& tss = meteo(MeteoData::TSS);
if (tss != IOUtils::nodata)
tss = C_TO_K(tss);
tss = IOUtils::C_TO_K(tss);
double& rh = meteo(MeteoData::RH);
if (rh != IOUtils::nodata)
......
......@@ -517,9 +517,9 @@ bool SNIO::parseMeteoLine(const std::vector<std::string>& vecLine, const std::st
if ((ea <= 1) && (ea != plugin_nodata)){
if ((md(MeteoData::TA) != plugin_nodata) && (md(MeteoData::RH) != plugin_nodata)) {
if(ea==0.)
ea = Atmosphere::Brutsaert_ilwr(md(MeteoData::RH)/100., C_TO_K(md(MeteoData::TA)));
ea = Atmosphere::Brutsaert_ilwr(md(MeteoData::RH)/100., IOUtils::C_TO_K(md(MeteoData::TA)));
else
ea = Atmosphere::Omstedt_ilwr(md(MeteoData::RH)/100., C_TO_K(md(MeteoData::TA)), ea); //calculate ILWR from cloudiness
ea = Atmosphere::Omstedt_ilwr(md(MeteoData::RH)/100., IOUtils::C_TO_K(md(MeteoData::TA)), ea); //calculate ILWR from cloudiness
} else {
ea = plugin_nodata;
}
......@@ -657,7 +657,7 @@ void SNIO::writeStationMeteo(const std::vector<MeteoData>& vecmd, const std::str
failure_count++;
fout << setw(6) << setprecision(0) << ta << " ";
} else
fout << setw(6) << setprecision(2) << K_TO_C(ta) << " ";
fout << setw(6) << setprecision(2) << IOUtils::K_TO_C(ta) << " ";
if(rh==IOUtils::nodata) {
failure_count++;
fout << setw(5) << setprecision(0) << rh << " ";
......@@ -701,13 +701,13 @@ void SNIO::writeStationMeteo(const std::vector<MeteoData>& vecmd, const std::str
Dirichlet_failure_count++;
fout << setw(7) << setprecision(1) << "0.0" << " ";
} else {
fout << setw(7) << setprecision(2) << K_TO_C(tss) << " ";
fout << setw(7) << setprecision(2) << IOUtils::K_TO_C(tss) << " ";
}
if(tsg==IOUtils::nodata) {
Dirichlet_failure_count++;
fout << setw(6) << setprecision(1) << "0.0" << " ";
} else {
fout << setw(6) << setprecision(2) << K_TO_C(tsg) << " ";
fout << setw(6) << setprecision(2) << IOUtils::K_TO_C(tsg) << " ";
}
//HNW, HS
......@@ -737,7 +737,7 @@ void SNIO::writeStationMeteo(const std::vector<MeteoData>& vecmd, const std::str
optional_failure_count++;
fout << setw(7) << setprecision(0) << ts << " ";
} else {
fout << setw(7) << setprecision(2) << K_TO_C(ts) << " ";
fout << setw(7) << setprecision(2) << IOUtils::K_TO_C(ts) << " ";
}
} else {
break;
......@@ -817,19 +817,19 @@ void SNIO::convertUnits(MeteoData& meteo)
double& ta = meteo(MeteoData::TA);
if (ta != IOUtils::nodata) {
if (ta < 100)
ta = C_TO_K(ta);
ta = IOUtils::C_TO_K(ta);
}
double& tsg = meteo(MeteoData::TSG);
if ( tsg != IOUtils::nodata) {
if (tsg < 100)
tsg = C_TO_K(tsg);
tsg = IOUtils::C_TO_K(tsg);
}
double& tss = meteo(MeteoData::TSS);
if (tss != IOUtils::nodata) {
if (tss < 100)
tss = C_TO_K(tss);
tss = IOUtils::C_TO_K(tss);
}
double& rh = meteo(MeteoData::RH);
......@@ -844,7 +844,7 @@ void SNIO::convertUnits(MeteoData& meteo)
ss << "TS" << ii;
if (meteo.param_exists(ss.str())){
double& value = meteo(ss.str());
value = C_TO_K(value);
value = IOUtils::C_TO_K(value);
} else {
break;
}
......
Markdown is supported
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