WSL/SLF GitLab Repository

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

The Winstral spatial interpolation algorithm is getting into shape: all three...

The Winstral spatial interpolation algorithm is getting into shape: all three modes of operation have been implemented and tested and it is now properly documented. A new data generator has been created that computes a relative humidity from either a dew point temperature or a specific humidity.  Adding a "check_attribute" to libncpp.
parent b3bd1947
......@@ -36,6 +36,8 @@ GeneratorAlgorithm* GeneratorAlgorithmFactory::getAlgorithm(const std::string& i
return new SinGenerator(vecArgs, i_algoname);
} else if (algoname == "STD_PRESS"){
return new StandardPressureGenerator(vecArgs, i_algoname);
} else if (algoname == "RELHUM"){
return new RhGenerator(vecArgs, i_algoname);
} else if (algoname == "CLEARSKY_LW"){
return new ClearSkyLWGenerator(vecArgs, i_algoname);
} else if (algoname == "ALLSKY_LW"){
......@@ -177,6 +179,73 @@ bool StandardPressureGenerator::generate(const size_t& param, std::vector<MeteoD
}
bool RhGenerator::generate(const size_t& param, MeteoData& md)
{
double &value = md(param);
if(value == IOUtils::nodata) {
const double TA = md(MeteoData::TA);
if (TA==IOUtils::nodata) //nothing else we can do here
return false;
//first chance to compute RH
if (md.param_exists("TD")) {
const double TD = md("TD");
if (TD!=IOUtils::nodata)
value = Atmosphere::DewPointtoRh(TD, TA, false);
}
//second chance to try to compute RH
if (value==IOUtils::nodata && md.param_exists("SH")) {
const double SH = md("SH");
const double altitude = md.meta.position.getAltitude();
if (SH!=IOUtils::nodata && altitude!=IOUtils::nodata)
value = Atmosphere::specToRelHumidity(altitude, TA, SH);
}
if (value==IOUtils::nodata) return false;
}
return true; //all missing values could be filled
}
bool RhGenerator::generate(const size_t& param, std::vector<MeteoData>& vecMeteo)
{
if(vecMeteo.empty()) return true;
const double altitude = vecMeteo.front().meta.position.getAltitude(); //if the stations move, this has to be in the loop
bool all_filled = true;
for(size_t ii=0; ii<vecMeteo.size(); ii++) {
double &value = vecMeteo[ii](param);
if(value == IOUtils::nodata) {
const double TA = vecMeteo[ii](MeteoData::TA);
if (TA==IOUtils::nodata) { //nothing else we can do here
all_filled=false;
continue;
}
//first chance to compute RH
if (vecMeteo[ii].param_exists("TD")) {
const double TD = vecMeteo[ii]("TD");
if (TD!=IOUtils::nodata)
value = Atmosphere::DewPointtoRh(TD, TA, false);
}
//second chance to try to compute RH
if (value==IOUtils::nodata && vecMeteo[ii].param_exists("SH")) {
const double SH = vecMeteo[ii]("SH");
if (SH!=IOUtils::nodata && altitude!=IOUtils::nodata)
value = Atmosphere::specToRelHumidity(altitude, TA, SH);
}
if (value==IOUtils::nodata) all_filled=false;
}
}
return all_filled;
}
void ClearSkyLWGenerator::parse_args(const std::vector<std::string>& vecArgs)
{
//Get the optional arguments for the algorithm: constant value to use
......
......@@ -70,6 +70,7 @@ namespace mio {
* @section generators_keywords Available generators
* The keywords defining the algorithms are the following:
* - STD_PRESS: standard atmospheric pressure as a function of the elevation of each station (see StandardPressureGenerator)
* - RELHUM: relative humidity from other humidity measurements (see RhGenerator)
* - CST: constant value as provided in argument (see ConstGenerator)
* - SIN: sinusoidal variation (see SinGenerator)
* - CLEARSKYLW: use a clear sky model to generate ILWR from TA, RH (see ClearSkyLWGenerator)
......@@ -191,6 +192,24 @@ class StandardPressureGenerator : public GeneratorAlgorithm {
bool generate(const size_t& param, std::vector<MeteoData>& vecMeteo);
};
/**
* @class RhGenerator
* @brief Relative humidity generator.
* Generate the relative humidity from either dew point temperature or specific humidity and air temperature.
* The dew point temperature must be named "TD" and the specific humidity "SH"
* @code
* RH::generators = RELHUM
* @endcode
*/
class RhGenerator : public GeneratorAlgorithm {
public:
RhGenerator(const std::vector<std::string>& vecArgs, const std::string& i_algo)
: GeneratorAlgorithm(vecArgs, i_algo) { parse_args(vecArgs); }
bool generate(const size_t& param, MeteoData& md);
bool generate(const size_t& param, std::vector<MeteoData>& vecMeteo);
};
/**
* @class ClearSkyLWGenerator
* @brief ILWR clear sky parametrization
......
......@@ -627,59 +627,50 @@ void RyanAlgorithm::calculate(const DEMObject& dem, Grid2DObject& grid)
const double WinstralAlgorithm::dmax = 300.;
double WinstralAlgorithm::getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param)
{
//This algorithm is only valid for HNW (we could add HS later)
if (in_param!=MeteoData::HNW)
return 0.0;
date = i_date;
param = in_param;
nrOfMeasurments = getData(date, param, vecData, vecMeta);
if (nrOfMeasurments==0)
return 0.0;
return 0.99;
}
void WinstralAlgorithm::getParameters(std::string &base_algo, double &synoptic_bearing) const
WinstralAlgorithm::WinstralAlgorithm(Meteo2DInterpolator& i_mi, const std::vector<std::string>& i_vecArgs,
const std::string& i_algo, IOManager& iom)
: InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom), base_algo("IDW_LAPSE"), ref_station(), user_synoptic_bearing(IOUtils::nodata)
{
//option 0: compute DW from multiple stations
//option 1: fixed DW
//option 2: specify the station_id to take DW from
base_algo = "IDW_LAPSE";
const size_t nr_args = vecArgs.size();
if (nr_args==0) {
synoptic_bearing = getSynopticBearing(vecMeteo);
if (synoptic_bearing==IOUtils::nodata)
throw IOException("Synoptic wind could not be computed by the "+algo+" algorithm", AT);
return;
} else if (nr_args==1) {
if (nr_args==1) {
if (IOUtils::isNumeric(vecArgs[0]))
IOUtils::convertString(synoptic_bearing, vecArgs[0]);
IOUtils::convertString(user_synoptic_bearing, vecArgs[0]);
else
throw InvalidArgumentException("Please provide both the base interpolation method and the station_id to use for wind direction for the "+algo+" algorithm", AT);
return;
} else if (nr_args==2) {
if (IOUtils::isNumeric(vecArgs[0])) {
IOUtils::convertString(synoptic_bearing, vecArgs[0]);
IOUtils::convertString(user_synoptic_bearing, vecArgs[0]);
base_algo = IOUtils::strToUpper( vecArgs[1] );
} else if (IOUtils::isNumeric(vecArgs[1])) {
IOUtils::convertString(synoptic_bearing, vecArgs[1]);
IOUtils::convertString(user_synoptic_bearing, vecArgs[1]);
base_algo = IOUtils::strToUpper( vecArgs[0] );
} else {
base_algo = IOUtils::strToUpper( vecArgs[0] );
const string ref_station( vecArgs[1] );
synoptic_bearing = getSynopticBearing(vecMeteo, ref_station, algo);
ref_station = vecArgs[1];
}
return;
} else
} else if (nr_args>2)
throw InvalidArgumentException("Wrong number of arguments supplied for the "+algo+" algorithm", AT);
}
void WinstralAlgorithm::initGrid(const std::string& base_algo, const DEMObject& dem, Grid2DObject& grid)
double WinstralAlgorithm::getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param)
{
//This algorithm is only valid for HNW (we could add HS later)
if (in_param!=MeteoData::HNW)
return 0.0;
date = i_date;
param = in_param;
nrOfMeasurments = getData(date, param, vecData, vecMeta);
if (nrOfMeasurments==0)
return 0.0;
return 0.99;
}
void WinstralAlgorithm::initGrid(const DEMObject& dem, Grid2DObject& grid)
{
//initialize precipitation grid with user supplied algorithm (IDW_LAPSE by default)
vector<string> vecArgs2;
......@@ -700,12 +691,51 @@ double WinstralAlgorithm::getSynopticBearing(const std::vector<MeteoData>& vecMe
throw InvalidArgumentException("Station \""+ref_station+"\" should be used as reference station for the "+algo+" algorithm but has not been found", AT);
}
double WinstralAlgorithm::getSynopticBearing(const std::vector<MeteoData>& vecMeteo)
//this method scans a square centered on the station for summmits
//that would shelter the station for wind.
//the sheltering criteria is: if ( height_difference*shade_factor > distance ) sheltered=true
bool WinstralAlgorithm::isExposed(const DEMObject& dem, Coords location)
{
// 1) locate the stations in DEM and check if they are higher than their surroundings within a given radius
// 2) simply compute a mean or median direction
// (2) can be used on all the stations selected in (1)
if (!dem.gridify(location)) {
return false;
}
const int i_ref = location.getGridI();
const int j_ref = location.getGridJ();
const double alt_ref = location.getAltitude()+4.; //4 m mast added so flat terrain behaves properly
const size_t ii_ref = static_cast<size_t>(i_ref);
const size_t jj_ref = static_cast<size_t>(j_ref);
const double shade_factor = 5.;
const double cellsize = dem.cellsize;
const double min_dh = cellsize/shade_factor; //min_dist=cellsize -> min_dh
const double search_dist = (dem.grid2D.getMax() - alt_ref) * shade_factor;
if (search_dist<=cellsize) return true;
const int search_idx = Optim::ceil( search_dist/cellsize );
const size_t i_min = static_cast<size_t>(std::max( 0, i_ref-search_idx ));
const size_t i_max = static_cast<size_t>(std::min( static_cast<int>(dem.getNx()-1), i_ref+search_idx ));
const size_t j_min = static_cast<size_t>(std::max( 0, j_ref-search_idx ));
const size_t j_max = static_cast<size_t>(std::min( static_cast<int>(dem.getNy()-1), j_ref+search_idx ));
for (size_t jj=j_min; jj<=j_max; jj++) {
for (size_t ii=i_min; ii<=i_max; ii++) {
if (ii==ii_ref && jj==jj_ref) continue; //skip the cell containing the station!
const double dh = dem.grid2D(ii,jj) - alt_ref;
if (dh<=min_dh) continue; //for negative or too small dh
const double distance = Optim::fastSqrt_Q3( Optim::pow2(static_cast<int>(ii)-i_ref) + Optim::pow2(static_cast<int>(jj)-j_ref) ) * cellsize;
if (distance<dh*shade_factor) {
return false;
}
}
}
return true;
}
double WinstralAlgorithm::getSynopticBearing(const std::vector<MeteoData>& vecMeteo)
{
double ve=0.0, vn=0.0;
size_t count=0;
for (size_t ii=0; ii<vecMeteo.size(); ii++) {
......@@ -717,19 +747,41 @@ double WinstralAlgorithm::getSynopticBearing(const std::vector<MeteoData>& vecMe
count++;
}
}
//HACK: use median instead of mean for exposed stations?
if (count!=0) {
ve /= static_cast<double>(count);
vn /= static_cast<double>(count);
//const double meanspeed = sqrt(ve*ve + vn*vn);
const double meandirection = fmod( atan2(ve,vn) * Cst::to_deg + 360., 360.);
return meandirection;
return Optim::round(meandirection/10.)*10; //round to nearest 10 deg
}
return IOUtils::nodata;
}
double WinstralAlgorithm::getSynopticBearing(const DEMObject& dem, const std::vector<MeteoData>& vecMeteo)
{
// 1) locate the stations in DEM and check if they are higher than their surroundings within a given radius
// 2) simply compute a mean or median direction
// (2) can be used on all the stations selected in (1)
std::vector<MeteoData> stationsSubset;
for (size_t ii=0; ii<vecMeteo.size(); ii++) {
if (isExposed(dem, vecMeteo[ii].meta.position)) {
//std::cout << "Station " << vecMeteo[ii].meta.stationID << " - " << vecMeteo[ii].meta.stationName << " -> DW=" << vecMeteo[ii](MeteoData::DW) << "\n";
stationsSubset.push_back( vecMeteo[ii] );
}
}
if (!stationsSubset.empty()) {
return getSynopticBearing(stationsSubset);
} else {
std::cerr << "[W] Synoptic wind direction computed from wind-sheltered stations only, the quality will strongly depend on their spatial distribution\n";
return getSynopticBearing(vecMeteo);
}
}
void WinstralAlgorithm::calculate(const DEMObject& dem, Grid2DObject& grid)
{
info.clear(); info.str("");
......@@ -740,10 +792,17 @@ void WinstralAlgorithm::calculate(const DEMObject& dem, Grid2DObject& grid)
return;
}
string base_algo, ref_station;
double synoptic_bearing;
getParameters(base_algo, synoptic_bearing);
initGrid(base_algo, dem, grid);
//HACK: also get Synoptic_VW and use a threshold to trigger snow transport?
double synoptic_bearing = user_synoptic_bearing;
if (synoptic_bearing==IOUtils::nodata) {
if (!ref_station.empty())
synoptic_bearing = getSynopticBearing(vecMeteo, ref_station, algo);
else
synoptic_bearing = getSynopticBearing(dem, vecMeteo);
}
info << "DW=" << synoptic_bearing << " - ";
initGrid(dem, grid);
//get TA interpolation from call back to Meteo2DInterpolator
Grid2DObject ta;
......
......@@ -85,6 +85,7 @@ class Meteo2DInterpolator; // forward declaration, cyclic header include
* - ILWR: the incoming long wave radiation is converted to emissivity and then interpolated (see ILWRAlgorithm)
* - WIND_CURV: the wind field (VW and DW) is interpolated using IDW_LAPSE and then altered depending on the local curvature and slope (taken from the DEM, see SimpleWindInterpolationAlgorithm)
* - RYAN: the wind direction is interpolated using IDW and then altered depending on the local slope (see RyanAlgorithm)
* - WINSTRAL: the solid precipitation is redistributed by wind according to (Winstral, 2002) (see WinstralAlgorithm)
* - HNW_SNOW: precipitation interpolation according to (Magnusson, 2011) (see SnowHNWInterpolation)
* - ODKRIG: ordinary kriging (see OrdinaryKrigingAlgorithm)
* - ODKRIG_LAPSE: ordinary kriging with lapse rate (see LapseOrdinaryKrigingAlgorithm)
......@@ -409,11 +410,20 @@ class RyanAlgorithm : public InterpolationAlgorithm {
* The DEM is used to compute wind exposure factors that are used to alter the precipitation fields.
* It is usually a good idea to provide a DEM that also contain the accumulated snow height in order
* to get a progressive softening of the terrain features.
* It takes the following arguments: the base algorithm to generate the initial wind field, the
* synoptic wind direction or the station_id of the station to get the wind direction from. All these arguments
* are optional but in order to make it possible to distinguish the arguments, it is not possible to only
* provide one argument of type "string" (in such a case, the second one must be provided and the order
* defines which is which).
*
* This method must therefore first use another algorithm to generate an initial precipitation field,
* and then modifies this field accordingly. This base method is "idw_lapse" by default.
* Then it requires a synoptic wind direction that can be provided by different means:
* - without any extra argument, the stations are located in the DEM and their wind shading (or exposure)
* is computed. If at least one station is found that is not sheltered from the wind (in every direction), it
* provides the synoptic wind (in case of multiple stations, the vector average is used). Please note that
* the stations that are not included in the DEM are considered to be sheltered. If no such station
* is found, the vector average of all the available stations is used.
* - by providing a fixed synoptic wind bearing that is used for all time steps
* - by providing the station_id of the station to get the wind direction from. In this case, the base algorithm
* for generating the initial wind field must be specified in the first position.
*
* @remarks Only cells with an air temperature below freezing participate in the redistribution
* @code
* HNW::algorithms = WINSTRAL
* HNW::winstral = idw_lapse 180
......@@ -422,17 +432,19 @@ class RyanAlgorithm : public InterpolationAlgorithm {
class WinstralAlgorithm : public InterpolationAlgorithm {
public:
WinstralAlgorithm(Meteo2DInterpolator& i_mi,
const std::vector<std::string>& i_vecArgs,
const std::string& i_algo, IOManager& iom)
: InterpolationAlgorithm(i_mi, i_vecArgs, i_algo, iom) {}
const std::vector<std::string>& i_vecArgs,
const std::string& i_algo, IOManager& iom);
virtual double getQualityRating(const Date& i_date, const MeteoData::Parameters& in_param);
virtual void calculate(const DEMObject& dem, Grid2DObject& grid);
private:
void getParameters(std::string &base_algo, double &synoptic_bearing) const;
void initGrid(const std::string& base_algo, const DEMObject& dem, Grid2DObject& grid);
void initGrid(const DEMObject& dem, Grid2DObject& grid);
static bool isExposed(const DEMObject& dem, Coords location);
static double getSynopticBearing(const std::vector<MeteoData>& vecMeteo, const std::string& ref_station, const std::string& algo);
static double getSynopticBearing(const std::vector<MeteoData>& vecMeteo);
static double getSynopticBearing(const DEMObject& dem, const std::vector<MeteoData>& vecMeteo);
std::string base_algo, ref_station;
double user_synoptic_bearing;
static const double dmax;
};
......
......@@ -699,7 +699,10 @@ void Interpol2D::WinstralSX(const DEMObject& dem, const double& dmax, const doub
* A linear correlation between erosion coefficients and eroded mass is assumed, that is that the points with maximum erosion get all their
* precipitation removed. The eroded mass is then distributed on the cells with positive Sx (with a linear correlation between positive Sx and deposited
* mass) and enforcing mass conservation within the domain.
* @remarks Only cells with an air temperature below freezing participate in the redistribution
*
* @param dem digital elevation model
* @param TA air temperature grid (in order to discriminate between solid and liquid precipitation)
* @param dmax search radius
* @param in_bearing wind direction to consider
* @param grid 2D array of precipitation to fill
......
......@@ -83,6 +83,16 @@ void get_attribute(const int& ncid, const std::string& varname, const int& varid
delete[] value;
}
bool check_attribute(const int& ncid, const int& varid, const std::string& attr_name)
{
size_t attr_len;
const int status = nc_inq_attlen (ncid, varid, attr_name.c_str(), &attr_len);
if (status != NC_NOERR) return false;
return true;
}
bool check_variable(const int& ncid, const std::string& varname)
{
int varid;
......
......@@ -45,6 +45,7 @@ namespace ncpp {
void add_attribute(const int& ncid, const int& varid, const std::string& attr_name, const std::string& attr_value);
void add_attribute(const int& ncid, const int& varid, const std::string& attr_name, const double& attr_value);
void get_attribute(const int& ncid, const std::string& varname, const int& varid, const std::string& attr_name, std::string& attr_value);
bool check_attribute(const int& ncid, const int& varid, const std::string& attr_name);
//Adding dimensions
void add_dimension(const int& ncid, const std::string& dimname, const size_t& length, int& dimid);
......
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