diff --git a/DATA/PROJECT/testFrostForecast/testFrostForecastSettings.ini b/DATA/PROJECT/testFrostForecast/testFrostForecastSettings.ini index 18b2a68c..9eca55eb 100644 --- a/DATA/PROJECT/testFrostForecast/testFrostForecastSettings.ini +++ b/DATA/PROJECT/testFrostForecast/testFrostForecastSettings.ini @@ -7,19 +7,25 @@ utm_zone=32 [project] csvPath=./outputData/ -id="4,10,2303,6257,6265" +id=4,10,2303,6257,6265 meteoGrid=../../METEOGRID/DBGridXML_ERG5_forecast.xml meteo_points=../../METEOPOINT/frost.db name=testFrostForecast -var="TAVG,RHAVG,W_SCAL_INT,W_VEC_DIR,PREC,RAD" +var=TAVG,RHAVG,W_SCAL_INT,W_VEC_DIR,PREC,RAD [reuter_param] -SE_intercept="0.182930, 0.249670, 0.169418, 0.332520, 0.274008" -SE_parRHss="0.001691, 0.002694, 0.001781, 0.003540, 0.002943" -SE_parTss="0.012427, 0.014025, 0.009732, 0.018650, 0.014905" -intercept="1.969822, 1.775462, 3.016488, 2.271830, 2.367134" -parRHss="-0.014889, -0.013332, -0.014434, -0.020210, -0.011935" -parTss="0.093746, 0.129308, 0.079078, 0.133310, 0.106645" +SE_intercept=0.182930, 0.249670, 0.169418, 0.332520, 0.274008 +SE_parRHss=0.001691, 0.002694, 0.001781, 0.003540, 0.002943 +SE_parTss=0.012427, 0.014025, 0.009732, 0.018650, 0.014905 +intercept=1.969822, 1.775462, 3.016488, 2.271830, 2.367134 +parRHss=-0.014889, -0.013332, -0.014434, -0.020210, -0.011935 +parTss=0.093746, 0.129308, 0.079078, 0.133310, 0.106645 + +[calibration] +monthIni=3 +monthFin=4 +thresholdTmin=5 +thresholdTrange=6 [software] software=frostForecast diff --git a/bin/frostForecast/frost.cpp b/bin/frostForecast/frost.cpp index 7302777b..22bbc34a 100644 --- a/bin/frostForecast/frost.cpp +++ b/bin/frostForecast/frost.cpp @@ -30,6 +30,11 @@ void Frost::initialize() SE_intercept.clear(); SE_parTss.clear(); SE_parRHss.clear(); + + monthIni = 1; + monthFin = 12; + thresholdTmin = 0; + thresholdTrange = 10; } void Frost::setSettingsFileName(const QString &value) @@ -143,92 +148,101 @@ int Frost::readSettings() csvFilePath = path + QDir::cleanPath(csvFilePath); } - QString idTemp = projectSettings->value("id","").toString(); - if (idTemp.isEmpty()) + idList = projectSettings->value("id").toStringList(); + if (idList.isEmpty()) { logger.writeError ("missing id station"); return ERROR_MISSINGPARAMETERS; } - else - { - idList = idTemp.split(","); - } - QString varTemp = projectSettings->value("var","").toString(); - if (varTemp.isEmpty()) + varList = projectSettings->value("var").toStringList(); + if (varList.isEmpty()) { logger.writeError ("missing var"); return ERROR_MISSINGPARAMETERS; } - else - { - varList = varTemp.split(","); - } projectSettings->endGroup(); + QStringList myList; + // reuter_param projectSettings->beginGroup("reuter_param"); - QString interceptTemp = projectSettings->value("intercept","").toString(); - if (interceptTemp.isEmpty()) + + myList = projectSettings->value("intercept").toStringList(); + if (myList.isEmpty()) { logger.writeError ("missing reuter param intercept"); return ERROR_MISSINGPARAMETERS; } else { - intercept = interceptTemp.split(","); + intercept = StringListToFloat(myList); } - QString parTssTemp = projectSettings->value("parTss","").toString(); - if (parTssTemp.isEmpty()) + + myList = projectSettings->value("parTss").toStringList(); + if (myList.isEmpty()) { logger.writeError ("missing reuter param parTss"); return ERROR_MISSINGPARAMETERS; } else { - parTss = parTssTemp.split(","); + parTss = StringListToFloat(myList); } - QString parRHssTemp = projectSettings->value("parRHss","").toString(); - if (parRHssTemp.isEmpty()) + + myList = projectSettings->value("parRHss").toStringList(); + if (myList.isEmpty()) { logger.writeError ("missing reuter param parRHss"); return ERROR_MISSINGPARAMETERS; } else { - parRHss = parRHssTemp.split(","); + parRHss = StringListToFloat(myList); } - QString SEinterceptTemp = projectSettings->value("SE_intercept","").toString(); - if (SEinterceptTemp.isEmpty()) + + myList = projectSettings->value("SE_intercept").toStringList(); + if (myList.isEmpty()) { logger.writeError ("missing reuter param SE_intercept"); return ERROR_MISSINGPARAMETERS; } else { - SE_intercept = SEinterceptTemp.split(","); + SE_intercept = StringListToFloat(myList); } - QString SEparTssTemp = projectSettings->value("SE_parTss","").toString(); - if (SEparTssTemp.isEmpty()) + + myList = projectSettings->value("SE_parTss").toStringList(); + if (myList.isEmpty()) { logger.writeError ("missing reuter param SE_parTss"); return ERROR_MISSINGPARAMETERS; } else { - SE_parTss = SEparTssTemp.split(","); + SE_parTss = StringListToFloat(myList); } - QString SEparRHssTemp = projectSettings->value("SE_parRHss","").toString(); - if (SEparRHssTemp.isEmpty()) + + myList = projectSettings->value("SE_parRHss").toStringList(); + if (myList.isEmpty()) { logger.writeError ("missing reuter param SE_parRHss"); return ERROR_MISSINGPARAMETERS; } else { - SE_parRHss = SEparRHssTemp.split(","); + SE_parRHss = StringListToFloat(myList); } + projectSettings->endGroup(); + + // calibration + projectSettings->beginGroup("calibration"); + monthIni = projectSettings->value("monthIni").toInt(); + monthFin = projectSettings->value("monthFin").toInt(); + thresholdTmin = projectSettings->value("thresholdTmin").toFloat(); + thresholdTrange = projectSettings->value("thresholdTrange").toFloat(); + projectSettings->endGroup(); return FROSTFORECAST_OK; } @@ -385,12 +399,12 @@ int Frost::getForecastData(QString id, int posIdList) int myTmpHour; QDateTime dateTimeTmp; - float myIntercept = intercept[posIdList].toFloat(); - float myParTss = parTss[posIdList].toFloat(); - float myParRHss = parRHss[posIdList].toFloat(); - float mySEintercept = SE_intercept[posIdList].toFloat(); - float mySEparTss = SE_parTss[posIdList].toFloat(); - float mySEparRHss = SE_parRHss[posIdList].toFloat(); + float myIntercept = intercept[posIdList]; + float myParTss = parTss[posIdList]; + float myParRHss = parRHss[posIdList]; + float mySEintercept = SE_intercept[posIdList]; + float mySEparTss = SE_parTss[posIdList]; + float mySEparRHss = SE_parRHss[posIdList]; for (int i = 0; i <= 1; i++) { @@ -550,8 +564,51 @@ QList Frost::getIdList() const return idList; } -bool Frost::getRadiativeCoolingHistory(QString id, float thresholdTmin, float thresholdTRange, int monthIni, int monthFin, int timeZone, - std::vector>& outData, std::vector >& sunsetData) +void fitCoolingCoefficient(std::vector >& hourly_series, int nrIterations, float tolerance, float minValue, float maxValue, std::vector & coeff) +{ + float myMin, myMax; + float mySum1, mySum2; + float myFirstError = 0; + float mySecondError = 0;; + int iteration; + float Tsunset; + int time; + + for(std::vector hourlyT : hourly_series) + { + + myMin = minValue; + myMax = maxValue; + iteration = 0; + time = 0; + do { + + mySum1 = 0; + mySum2 = 0; + Tsunset = hourlyT[0]; + time = 0; + for (float T : hourlyT) + { + mySum1 += pow((Tsunset - myMin * sqrt(time) - T), 2); + mySum2 += pow((Tsunset - myMax * sqrt(time) - T), 2); + time++; + } + + myFirstError = sqrt(mySum1); + mySecondError = sqrt(mySum2); + + if (myFirstError < mySecondError) + myMax = (myMin + myMax) / 2; + else + myMin = (myMin + myMax) / 2; + + } while (fabs(myFirstError - mySecondError) > tolerance && iteration <= nrIterations); + + coeff.push_back((myMin + myMax) / 2); + } +} + +bool Frost::getRadiativeCoolingHistory(QString id, std::vector>& outData, std::vector >& sunsetData) { unsigned pos = 0; bool found = false; @@ -616,7 +673,7 @@ bool Frost::getRadiativeCoolingHistory(QString id, float thresholdTmin, float th { hourSunSetLocal = round(sunPosition.set/3600); hourSunRiseLocal = round(sunPosition.rise/3600); - hourSunRiseUtc = hourSunRiseLocal - timeZone; + hourSunRiseUtc = hourSunRiseLocal - gisSettings.timeZone; dateIndexSunRiseUtc = i; if (hourSunRiseUtc < 1) @@ -630,7 +687,7 @@ bool Frost::getRadiativeCoolingHistory(QString id, float thresholdTmin, float th dateIndexSunRiseUtc--; } - h = hourSunSetLocal - timeZone; + h = hourSunSetLocal - gisSettings.timeZone; d = i-1; isSunRise = false; @@ -691,7 +748,7 @@ bool Frost::getRadiativeCoolingHistory(QString id, float thresholdTmin, float th if (h == hourSunRiseUtc && d == dateIndexSunRiseUtc) { isSunRise = true; - if (dataPresent && minT < thresholdTmin && (maxT - minT >= thresholdTRange)) + if (dataPresent && minT < thresholdTmin && (maxT - minT >= thresholdTrange)) { outData.push_back(hourlyT); sunsetData.push_back(ssData); @@ -709,47 +766,26 @@ bool Frost::getRadiativeCoolingHistory(QString id, float thresholdTmin, float th return (outData.size() > 0); } - -void Frost::fitCoolingCoefficient(std::vector >& hourly_series, int nrIterations, float tolerance, float minValue, float maxValue, std::vector & coeff) +bool Frost::calibrateModel(QString id) { - float myMin, myMax; - float mySum1, mySum2; - float myFirstError = 0; - float mySecondError = 0;; - int iteration; - float Tsunset; - int time; + std::vector > outData; + std::vector > sunsetData; + std::vector outCoeff; + std::vector weights; + float regrConst, R2, stdError; + std::vector regrCoeff; - for(std::vector hourlyT : hourly_series) + if (getRadiativeCoolingHistory(id, outData, sunsetData)) { + fitCoolingCoefficient(outData, 10000, float(0.001), float(0), float(10), outCoeff); - myMin = minValue; - myMax = maxValue; - iteration = 0; - time = 0; - do { + for (int j=0; j < sunsetData.size(); j++) + weights.push_back(1); - mySum1 = 0; - mySum2 = 0; - Tsunset = hourlyT[0]; - time = 0; - for (float T : hourlyT) - { - mySum1 += pow((Tsunset - myMin * sqrt(time) - T), 2); - mySum2 += pow((Tsunset - myMax * sqrt(time) - T), 2); - time++; - } + statistics::weightedMultiRegressionLinearWithStats(sunsetData, outCoeff, weights, ®rConst, regrCoeff, false, true, &R2, &stdError); - myFirstError = sqrt(mySum1); - mySecondError = sqrt(mySum2); - - if (myFirstError < mySecondError) - myMax = (myMin + myMax) / 2; - else - myMin = (myMin + myMax) / 2; - - } while (fabs(myFirstError - mySecondError) > tolerance && iteration <= nrIterations); - - coeff.push_back((myMin + myMax) / 2); + return true; } + + return false; } diff --git a/bin/frostForecast/frost.h b/bin/frostForecast/frost.h index 3c8d65c4..cf57623f 100644 --- a/bin/frostForecast/frost.h +++ b/bin/frostForecast/frost.h @@ -39,9 +39,8 @@ class Frost int createCsvFile(QString id); QList getIdList() const; - // calibration functions - bool getRadiativeCoolingHistory(QString id, float thresholdTmin, float thresholdTRange, int monthIni, int monthFin, int timeZone, std::vector > &outData, std::vector > &sunsetData); - void fitCoolingCoefficient(std::vector > &hourly_series, int nrIterations, float tolerance, float minValue, float maxValue, std::vector &coeff); + bool getRadiativeCoolingHistory(QString id, std::vector>& outData, std::vector >& sunsetData); + bool calibrateModel(QString id); private: gis::Crit3DGisSettings gisSettings; @@ -58,16 +57,24 @@ class Frost QList meteoPointsList; QList idList; QList varList; - QList intercept; - QList parTss; - QList parRHss; - QList SE_intercept; - QList SE_parTss; - QList SE_parRHss; + + std::vector intercept; + std::vector parTss; + std::vector parRHss; + std::vector SE_intercept; + std::vector SE_parTss; + std::vector SE_parRHss; + int indexSunSet; int indexSunRise; QString errorString; + //calibration parameters + int monthIni; + int monthFin; + float thresholdTmin; + float thresholdTrange; + QList myForecast; QList myForecastMin; QList myForecastMax; diff --git a/bin/frostForecast/main.cpp b/bin/frostForecast/main.cpp index fa7fa155..b6b59f52 100644 --- a/bin/frostForecast/main.cpp +++ b/bin/frostForecast/main.cpp @@ -66,37 +66,24 @@ int main(int argc, char *argv[]) return result; } - /*result = frost.downloadMeteoPointsData(); + /* + result = frost.downloadMeteoPointsData(); if (result!=FROSTFORECAST_OK) { return result; }*/ + bool calibrate = true; + QList idList = frost.getIdList(); for (int i = 0; i< idList.size(); i++) { - std::vector > outData; - std::vector > sunsetData; - std::vector outCoeff; - std::vector weights; - float regrConst, R2, stdError; - std::vector regrCoeff; - - if (frost.getRadiativeCoolingHistory(idList[i], 5, 6, 3, 4, 1, outData, sunsetData)) - { - frost.fitCoolingCoefficient(outData, 10000, float(0.001), float(0), float(10), outCoeff); + if (calibrate) frost.calibrateModel(idList[i]); - for (int j=0; j < sunsetData.size(); j++) - weights.push_back(1); - - statistics::weightedMultiRegressionLinearWithStats(sunsetData, outCoeff, weights, ®rConst, regrCoeff, false, true, &R2, &stdError); - - /* - result = frost.getForecastData(idList[i], i); + /*result = frost.getForecastData(idList[i], i); if (result == FROSTFORECAST_OK) { frost.createCsvFile(idList[i]); }*/ - } } }