WSL/SLF GitLab Repository

Commit 8ecdc212 authored by Mathias Bavay's avatar Mathias Bavay
Browse files

Finally, another method for handling daily sums of solar radiation has been...

Finally, another method for handling daily sums of solar radiation has been implemented: when a point is requested, the potential solar radiation is computed for the whole day with a 20 minutes resolution, filling a vector and computing the daily sum. This enables to compute the loss factor. When another point within the same day is requested, the radiation is interpolated from the values stored in the vector and multiplied by the loss factor. Tested on WFJ data, this works very well (of course, days that experience large cloudiness variations are not reproduced well but at least the average flux is correct).

A typo has been fixed in the documentation of Date. 
parent 77c85422
......@@ -73,7 +73,7 @@ class Date {
///Keywords for selecting the date formats
typedef enum {
ISO, ///< ISO 8601 extended format combined date: YYYY-MM-DDTHH:mm:SS (fields might be dropped, in the least to the most significant order)
ISO_TZ, ///< ISO 8601 format (same as ISO) but with tie zone specification
ISO_TZ, ///< ISO 8601 format (same as ISO) but with time zone specification
FULL, ///< ISO 8601 followed by the julian date (in parenthesis)
NUM, ///< ISO 8601 basic format date: YYYYMMDDHHmmSS (fields might be dropped, in the least to the most significant order)
DIN ///<DIN5008 format: DD.MM.YYYY HH:MM
......
......@@ -19,6 +19,7 @@
#include <meteoio/MathOptim.h>
#include <meteoio/meteolaws/Atmosphere.h>
#include <meteoio/meteolaws/Sun.h>
#include <meteoio/meteostats/libinterpol1D.h>
#include <cmath>
#include <algorithm>
......@@ -147,7 +148,7 @@ std::string NoResampling::toString() const
}
void NoResampling::resample(const size_t& index, const ResamplingPosition& position, const size_t& paramindex,
const std::vector<MeteoData>& vecM, MeteoData& md) const
const std::vector<MeteoData>& vecM, MeteoData& md)
{
if (index >= vecM.size())
throw IOException("The index of the element to be resampled is out of bounds", AT);
......@@ -195,7 +196,7 @@ std::string NearestNeighbour::toString() const
}
void NearestNeighbour::resample(const size_t& index, const ResamplingPosition& position, const size_t& paramindex,
const std::vector<MeteoData>& vecM, MeteoData& md) const
const std::vector<MeteoData>& vecM, MeteoData& md)
{
if (index >= vecM.size())
throw IOException("The index of the element to be resampled is out of bounds", AT);
......@@ -276,7 +277,7 @@ std::string LinearResampling::toString() const
}
void LinearResampling::resample(const size_t& index, const ResamplingPosition& position, const size_t& paramindex,
const std::vector<MeteoData>& vecM, MeteoData& md) const
const std::vector<MeteoData>& vecM, MeteoData& md)
{
if (index >= vecM.size())
throw IOException("The index of the element to be resampled is out of bounds", AT);
......@@ -369,7 +370,7 @@ std::string Accumulate::toString() const
}
void Accumulate::resample(const size_t& index, const ResamplingPosition& position, const size_t& paramindex,
const std::vector<MeteoData>& vecM, MeteoData& md) const
const std::vector<MeteoData>& vecM, MeteoData& md)
{
if (index >= vecM.size())
throw IOException("The index of the element to be resampled is out of bounds", AT);
......@@ -435,8 +436,10 @@ void Accumulate::resample(const size_t& index, const ResamplingPosition& positio
const double Daily_solar::soil_albedo = .23; //grass
const double Daily_solar::snow_albedo = .85; //snow
const double Daily_solar::snow_thresh = .1; //if snow height greater than this threshold -> snow albedo
const size_t Daily_solar::samples_per_day = 24*3; //every 20 minutes
Daily_solar::Daily_solar(const std::string& i_algoname, const std::string& i_parname, const double& dflt_window_size, const std::vector<std::string>& vecArgs)
: ResamplingAlgorithms(i_algoname, i_parname, dflt_window_size, vecArgs)
: ResamplingAlgorithms(i_algoname, i_parname, dflt_window_size, vecArgs), radiation(), station_index(), dateStart(), dateEnd(), loss_factor()
{
const size_t nr_args = vecArgs.size();
if(nr_args>0) {
......@@ -451,18 +454,15 @@ std::string Daily_solar::toString() const
return ss.str();
}
size_t Daily_solar::getNearestValidPt(const size_t& pos, const size_t& paramindex, const std::vector<MeteoData>& vecM, const Date& resampling_date)
//look for the daily sum of solar radiation for the current day
size_t Daily_solar::getNearestValidPt(const std::vector<MeteoData>& vecM, const size_t& paramindex, const size_t& stat_idx, const size_t& pos) const
{
size_t indexP1=IOUtils::npos;
size_t indexP2=IOUtils::npos;
//look for daily sum before the current point
Date dateStart( resampling_date );
dateStart.rnd(24*3600, Date::DOWN);
if(dateStart==resampling_date) //if resampling_date=midnight GMT, the rounding lands on the exact same date
dateStart -= 1.;
for (size_t ii=pos; ii-- >0; ) {
if (vecM[ii].date < dateStart) break;
if (vecM[ii].date < dateStart[stat_idx]) break;
if (vecM[ii](paramindex) != IOUtils::nodata){
indexP1 = ii;
break;
......@@ -470,26 +470,97 @@ size_t Daily_solar::getNearestValidPt(const size_t& pos, const size_t& paraminde
}
//look for daily sum after the current point
Date dateEnd( resampling_date );
dateEnd.rnd(24*3600, Date::UP);
for (size_t ii=pos; ii<vecM.size(); ++ii) {
if (vecM[ii].date > dateEnd) break;
if (vecM[ii].date > dateEnd[stat_idx]) break;
if (vecM[ii](paramindex) != IOUtils::nodata) {
indexP2 = ii;
break;
}
}
if(indexP1!=IOUtils::npos && indexP2!=IOUtils::npos)
throw IOException("More than one daily sum of solar radiation found for "+resampling_date.toString(Date::ISO)+" !", AT);
if(indexP1!=IOUtils::npos && indexP2!=IOUtils::npos) {
const string msg = "More than one daily sum of solar radiation found between "+dateStart[stat_idx].toString(Date::ISO)+" and "+dateEnd[stat_idx].toString(Date::ISO);
throw IOException(msg, AT);
}
if(indexP1!=IOUtils::npos) return indexP1;
if(indexP2!=IOUtils::npos) return indexP2;
return IOUtils::npos;
}
void Daily_solar::resample(const size_t& index, const ResamplingPosition& position, const size_t& paramindex,
const std::vector<MeteoData>& vecM, MeteoData& md) const
//compute the daily sum of solar radiation as well as fill a vector containing the solar radiation at a regular sampling rate
double Daily_solar::compRadiation(const double& lat, const double& lon, const double& alt, const double& HS, const size_t& stat_idx)
{
const double P = Atmosphere::stdAirPressure(alt);
double albedo = 0.5;
if (HS!=IOUtils::nodata) //no big deal if we can not adapt the albedo
albedo = (HS>=snow_thresh)? snow_albedo : soil_albedo;
const double TA=274.98, RH=0.666; //the reduced precipitable water will get an average value
SunObject sun(lat, lon, alt);
double sum = 0.;
size_t index=0;
for(Date date(dateStart[stat_idx]); date<dateEnd[stat_idx]; date += 1./double(samples_per_day)) {
//compute potential solar radiation at this time step
sun.setDate(date.getJulian(), date.getTimeZone());
sun.calculateRadiation(TA, RH, P, albedo);
double toa, direct, diffuse;
sun.getHorizontalRadiation(toa, direct, diffuse);
const double global_h = direct+diffuse;
//compute the integral by a simple triangle method
sum += global_h * static_cast<double>(24*3600/samples_per_day);
//store the radiation for later reuse
radiation[stat_idx][index++] = global_h;
}
return sum;
}
//interpolate the solar radiation from the vector containing regularly sampled solar radiation
double Daily_solar::getSolarInterpol(const Date& resampling_date, const size_t& stat_idx) const
{
const double in_day = (resampling_date.getJulian() - dateStart[stat_idx].getJulian()) / (dateEnd[stat_idx].getJulian() - dateStart[stat_idx].getJulian());
const size_t vec_index = static_cast<size_t>( Optim::floor(in_day*samples_per_day) ); //so the sample will be between vec_index and vec_index+1
const double weight = in_day*static_cast<double>(samples_per_day) - static_cast<double>(vec_index);
return Interpol1D::weightedMean(radiation[stat_idx][vec_index], radiation[stat_idx][vec_index+1], weight);
}
//a new, previously unknown station has been found, allocate the memory
size_t Daily_solar::getStationIndex(const std::string& key)
{
const size_t nr_stations = station_index.size();
for(size_t ii=0; ii<nr_stations; ++ii) {
if(station_index[ii]==key)
return ii;
}
radiation.push_back( vector<double>(samples_per_day, 0.) );
loss_factor.push_back( 0. );
Date null_date(0., 0.);
dateStart.push_back( null_date );
dateEnd.push_back( null_date );
station_index.push_back( key );
return nr_stations; //the new index is the old nr_stations
}
void Daily_solar::setDayStartAndEnd(const Date& resampling_date, const size_t& stat_idx)
{
dateStart[stat_idx] = resampling_date;
dateStart[stat_idx].rnd(24*3600, Date::DOWN);
if(dateStart[stat_idx]==resampling_date) //if resampling_date=midnight GMT, the rounding lands on the exact same date
dateStart[stat_idx] -= 1.;
dateEnd[stat_idx] = resampling_date;
dateEnd[stat_idx].rnd(24*3600, Date::UP);
}
void Daily_solar::resample(const size_t& index, const ResamplingPosition& /*position*/, const size_t& paramindex,
const std::vector<MeteoData>& vecM, MeteoData& md)
{
if (index >= vecM.size())
throw IOException("The index of the element to be resampled is out of bounds", AT);
......@@ -497,44 +568,43 @@ void Daily_solar::resample(const size_t& index, const ResamplingPosition& positi
if(paramindex!=MeteoData::ISWR && paramindex!=MeteoData::RSWR)
throw IOException("This method only applies to short wave radiation! (either ISWR or RSWR)", AT);
md.setResampled(true);
const size_t indexP = getNearestValidPt(index, paramindex, vecM, md.date);
if(indexP==IOUtils::npos) //no daily sum found for the current day
return;
const double lat = md.meta.position.getLat();
const double lon = md.meta.position.getLon();
const double alt = md.meta.position.getAltitude();
if(lat==IOUtils::nodata || lon==IOUtils::nodata || alt==IOUtils::nodata) return;
const double HS = md(MeteoData::HS);
double P = md(MeteoData::P);
double TA = md(MeteoData::TA);
double RH = md(MeteoData::RH);
double HS = md(MeteoData::HS);
double albedo = 0.5;
//get station index
const size_t stat_idx = getStationIndex(md.meta.stationID);
if (P==IOUtils::nodata) P = Atmosphere::stdAirPressure(alt);
if (HS!=IOUtils::nodata) //no big deal if we can not adapt the albedo
albedo = (HS>=snow_thresh)? snow_albedo : soil_albedo;
if(TA==IOUtils::nodata || RH==IOUtils::nodata) {
//set TA & RH so the reduced precipitable water will get an average value
TA=274.98;
RH=0.666;
//has the radiation already been calculated for this day and station?
if(md.date<dateStart[stat_idx] || md.date>=dateEnd[stat_idx]) {
setDayStartAndEnd(md.date, stat_idx);
const size_t indexP = getNearestValidPt(vecM, paramindex, stat_idx, index);
if(indexP==IOUtils::npos) { //no daily sum found for the current day
loss_factor[stat_idx] = IOUtils::nodata;
return;
}
const double daily_sum = compRadiation(lat, lon, alt, HS, stat_idx);
loss_factor[stat_idx] = (daily_sum>0)? vecM[indexP](paramindex) / daily_sum : 0.; //in case of polar night...
}
SunObject sun(lat, lon, alt, md.date.getJulian(), md.date.getTimeZone());
sun.calculateRadiation(TA, RH, P, albedo);
double toa, direct, diffuse;
sun.getHorizontalRadiation(toa, direct, diffuse);
const double sum_toa_h = sun.approxTOADailySum();
//const double loss_factor = (sum_toa_h>0.)? vecM[indexP](paramindex) / sum_toa_h : 0.;
const double loss_factor = 1.;
if(loss_factor[stat_idx]==IOUtils::nodata) //the station could not be calculated for this day
return;
//interpolate radiation for this timestep and write it out
const double rad = getSolarInterpol(md.date, stat_idx);
if(paramindex==MeteoData::ISWR) {
md(paramindex) = loss_factor[stat_idx] * rad;
} else {
double albedo = 0.5;
if (HS!=IOUtils::nodata) //no big deal if we can not adapt the albedo
albedo = (HS>=snow_thresh)? snow_albedo : soil_albedo;
md(paramindex) = loss_factor[stat_idx] * rad * albedo;
}
if(paramindex==MeteoData::ISWR)
md(paramindex) = loss_factor * (toa);
else
md(paramindex) = loss_factor * (toa) * albedo;
md.setResampled(true);
}
} //namespace
......
......@@ -68,6 +68,7 @@ namespace mio {
* - nearest: nearest neighbor data resampling, see NearestNeighbour
* - linear: linear data resampling, see LinearResampling
* - accumulate: data re-accumulation as suitable for precipitations, see Accumulate
* - daily_solar: generate solar radiation (ISWR or RSWR) from daily sums, see Daily_solar
*/
/**
......@@ -98,7 +99,7 @@ class ResamplingAlgorithms {
std::string getAlgo() const {return algo;};
virtual void resample(const size_t& index, const ResamplingPosition& position, const size_t& paramindex,
const std::vector<MeteoData>& vecM, MeteoData& md) const = 0;
const std::vector<MeteoData>& vecM, MeteoData& md) = 0;
virtual std::string toString() const = 0;
......@@ -136,7 +137,7 @@ class NoResampling : public ResamplingAlgorithms {
NoResampling(const std::string& i_algoname, const std::string& i_parname, const double& dflt_window_size, const std::vector<std::string>& vecArgs);
void resample(const size_t& index, const ResamplingPosition& position, const size_t& paramindex,
const std::vector<MeteoData>& vecM, MeteoData& md) const;
const std::vector<MeteoData>& vecM, MeteoData& md);
std::string toString() const;
};
......@@ -158,7 +159,7 @@ class NearestNeighbour : public ResamplingAlgorithms {
NearestNeighbour(const std::string& i_algoname, const std::string& i_parname, const double& dflt_window_size, const std::vector<std::string>& vecArgs);
void resample(const size_t& index, const ResamplingPosition& position, const size_t& paramindex,
const std::vector<MeteoData>& vecM, MeteoData& md) const;
const std::vector<MeteoData>& vecM, MeteoData& md);
std::string toString() const;
private:
bool extrapolate;
......@@ -181,7 +182,7 @@ class LinearResampling : public ResamplingAlgorithms {
LinearResampling(const std::string& i_algoname, const std::string& i_parname, const double& dflt_window_size, const std::vector<std::string>& vecArgs);
void resample(const size_t& index, const ResamplingPosition& position, const size_t& paramindex,
const std::vector<MeteoData>& vecM, MeteoData& md) const;
const std::vector<MeteoData>& vecM, MeteoData& md);
std::string toString() const;
private:
bool extrapolate;
......@@ -204,7 +205,7 @@ class Accumulate : public ResamplingAlgorithms {
Accumulate(const std::string& i_algoname, const std::string& i_parname, const double& dflt_window_size, const std::vector<std::string>& vecArgs);
void resample(const size_t& index, const ResamplingPosition& position, const size_t& paramindex,
const std::vector<MeteoData>& vecM, MeteoData& md) const;
const std::vector<MeteoData>& vecM, MeteoData& md);
std::string toString() const;
private:
double accumulate_period; //internally, in julian days
......@@ -212,9 +213,11 @@ class Accumulate : public ResamplingAlgorithms {
};
/**
* @brief Generate solar radiation out of daily sum
* @brief Generate solar radiation out of daily sums.
* Daily sums of solar radiation (once, per day, any time during the day) are compared to the potential radiation, leading to an atmospheric loss factor.
* This loss factor is then applied to the potential solar radiation calculated at the requested time.
* When using this algorithm for RSWR, an albedo is required. A default value of 0.5 is used. If the snow height is available and greater than a 10cm threshold,
* a snow albedo is used. Below this threshold, a soil albedo is used.
* @code
* ISWR::resample = daily_solar
* @endcode
......@@ -224,11 +227,22 @@ class Daily_solar : public ResamplingAlgorithms {
Daily_solar(const std::string& i_algoname, const std::string& i_parname, const double& dflt_window_size, const std::vector<std::string>& vecArgs);
void resample(const size_t& index, const ResamplingPosition& position, const size_t& paramindex,
const std::vector<MeteoData>& vecM, MeteoData& md) const;
const std::vector<MeteoData>& vecM, MeteoData& md);
std::string toString() const;
private:
static size_t getNearestValidPt(const size_t& pos, const size_t& paramindex, const std::vector<MeteoData>& vecM, const Date& resampling_date);
size_t getNearestValidPt(const std::vector<MeteoData>& vecM, const size_t& paramindex, const size_t& stat_idx, const size_t& pos) const;
double getSolarInterpol(const Date& resampling_date, const size_t& stat_idx) const;
double compRadiation(const double& lat, const double& lon, const double& alt, const double& HS, const size_t& stat_idx);
size_t getStationIndex(const std::string& key);
void setDayStartAndEnd(const Date& resampling_date, const size_t& stat_idx);
std::vector< std::vector<double> > radiation;
std::vector<std::string> station_index;
std::vector<Date> dateStart, dateEnd;
std::vector<double> loss_factor;
static const double soil_albedo, snow_albedo, snow_thresh;
static const size_t samples_per_day;
};
} //end namespace
......
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