Merge commit '68510b41084cfaac61cecb69e50bc5b0d125687b'

agrolib/dbMeteoGrid/dbMeteoGrid.cpp

@@ -34,7 +34,7 @@ bool Crit3DMeteoGridDbHandler::parseXMLFile(QString xmlFileName, QDomDocument* x
     if (!myFile.open(QIODevice::ReadOnly))
     {
         *error = "Open XML failed:\n" + xmlFileName + "\n" + myFile.errorString();
-        return (false);
+        return false;
     }
 
     QString errorStr;
@@ -46,7 +46,7 @@ bool Crit3DMeteoGridDbHandler::parseXMLFile(QString xmlFileName, QDomDocument* x
                 + " - Column: " + QString::number(myErrColumn)
                 + "\n" + errorStr;
         myFile.close();
-        return(false);
+        return false;
     }
 
     myFile.close();
@@ -1908,8 +1908,9 @@ bool Crit3DMeteoGridDbHandler::loadGridDailyDataEnsemble(QString &errorStr, QStr
     int numberOfDays = first.daysTo(last) + 1;
     _meteoGrid->meteoPointPointer(row,col)->initializeObsDataD(numberOfDays, getCrit3DDate(first));
 
-    QString statement = QString("SELECT * FROM `%1` WHERE `%2`>= '%3' AND `%2`<= '%4' AND MemberNr = '%5' ORDER BY `%2`").arg(tableD, _tableDaily.fieldTime).arg(first.toString("yyyy-MM-dd")).arg(last.toString("yyyy-MM-dd")).arg(memberNr);
-    if( !qry.exec(statement) )
+    QString statement = QString("SELECT * FROM `%1` WHERE `%2`>= '%3' AND `%2`<= '%4' AND MemberNr = '%5' ORDER BY `%2`")
+                            .arg(tableD, _tableDaily.fieldTime, first.toString("yyyy-MM-dd"), last.toString("yyyy-MM-dd")).arg(memberNr);
+    if(! qry.exec(statement) )
     {
         errorStr = qry.lastError().text();
         return false;
@@ -2088,10 +2089,8 @@ bool Crit3DMeteoGridDbHandler::loadGridHourlyDataEnsemble(QString &errorStr, QSt
     int varCode;
     float value;
 
-    unsigned row;
-    unsigned col;
-
-    if (!_meteoGrid->findMeteoPointFromId(&row, &col, meteoPoint.toStdString()) )
+    unsigned row, col;
+    if (! _meteoGrid->findMeteoPointFromId(&row, &col, meteoPoint.toStdString()) )
     {
         errorStr = "Missing MeteoPoint id";
         return false;
@@ -2101,9 +2100,9 @@ bool Crit3DMeteoGridDbHandler::loadGridHourlyDataEnsemble(QString &errorStr, QSt
     _meteoGrid->meteoPointPointer(row, col)->initializeObsDataH(1, numberOfDays, getCrit3DDate(first.date()));
 
     QString statement = QString("SELECT * FROM `%1` WHERE `%2` >= '%3' AND `%2` <= '%4' AND MemberNr = '%5' ORDER BY `%2`")
-                                .arg(tableH, _tableHourly.fieldTime).arg(first.toString("yyyy-MM-dd hh:mm")).arg(last.toString("yyyy-MM-dd hh:mm")).arg(memberNr);
+                                .arg(tableH, _tableHourly.fieldTime, first.toString("yyyy-MM-dd hh:mm"), last.toString("yyyy-MM-dd hh:mm")).arg(memberNr);
 
-    if( !qry.exec(statement) )
+    if(! qry.exec(statement) )
     {
         errorStr = qry.lastError().text();
     }
@@ -2139,7 +2138,7 @@ bool Crit3DMeteoGridDbHandler::loadGridHourlyDataEnsemble(QString &errorStr, QSt
 }
 
 
-bool Crit3DMeteoGridDbHandler::loadGridHourlyDataFixedFields(QString &errorStr, QString meteoPoint, QDateTime first, QDateTime last)
+bool Crit3DMeteoGridDbHandler::loadGridHourlyDataFixedFields(QString &errorStr, const QString &meteoPoint, const QDateTime &first, const QDateTime &last)
 {
     errorStr = "";
 
@@ -2149,10 +2148,8 @@ bool Crit3DMeteoGridDbHandler::loadGridHourlyDataFixedFields(QString &errorStr,
     int varCode;
     float value;
 
-    unsigned row;
-    unsigned col;
-
-    if (!_meteoGrid->findMeteoPointFromId(&row, &col, meteoPoint.toStdString()) )
+    unsigned row, col;
+    if (! _meteoGrid->findMeteoPointFromId(&row, &col, meteoPoint.toStdString()) )
     {
         errorStr = "Missing MeteoPoint id";
         return false;
@@ -2161,8 +2158,9 @@ bool Crit3DMeteoGridDbHandler::loadGridHourlyDataFixedFields(QString &errorStr,
     int numberOfDays = first.date().daysTo(last.date());
     _meteoGrid->meteoPointPointer(row, col)->initializeObsDataH(1, numberOfDays, getCrit3DDate(first.date()));
 
-    QString statement = QString("SELECT * FROM `%1` WHERE `%2` >= '%3' AND `%2`<= '%4' ORDER BY `%2`").arg(tableH, _tableHourly.fieldTime).arg(first.toString("yyyy-MM-dd hh:mm")).arg(last.toString("yyyy-MM-dd hh:mm"));
-    if( !qry.exec(statement) )
+    QString statement = QString("SELECT * FROM `%1` WHERE `%2` >= '%3' AND `%2`<= '%4' ORDER BY `%2`")
+                            .arg(tableH, _tableHourly.fieldTime, first.toString("yyyy-MM-dd hh:mm"), last.toString("yyyy-MM-dd hh:mm"));
+    if(! qry.exec(statement) )
     {
         errorStr = qry.lastError().text();
     }
@@ -2189,9 +2187,7 @@ bool Crit3DMeteoGridDbHandler::loadGridHourlyDataFixedFields(QString &errorStr,
                 if (! _meteoGrid->meteoPointPointer(row,col)->setMeteoPointValueH(getCrit3DDate(date.date()), date.time().hour(), date.time().minute(), variable, value))
                     return false;
             }
-
         }
-
     }
 
     return true;

agrolib/dbMeteoGrid/dbMeteoGrid.h

@@ -122,7 +122,7 @@
         bool loadGridDailyDataEnsemble(QString &errorStr, QString meteoPoint, int memberNr, QDate first, QDate last);
         bool loadGridDailyMeteoPrec(QString &errorStr, const QString &meteoPointId, const QDate &firstDate, const QDate &lastDate);
         bool loadGridHourlyData(QString &errorStr, QString meteoPoint, QDateTime firstDate, QDateTime lastDate);
-        bool loadGridHourlyDataFixedFields(QString &errorStr, QString meteoPoint, QDateTime first, QDateTime last);
+        bool loadGridHourlyDataFixedFields(QString &errorStr, const QString &meteoPoint, const QDateTime &first, const QDateTime &last);
         bool loadGridHourlyDataEnsemble(QString &errorStr, QString meteoPoint, int memberNr, QDateTime first, QDateTime last);
         bool loadGridMonthlyData(QString &errorStr, QString meteoPoint, QDate firstDate, QDate lastDate);
         bool loadGridAllMonthlyData(QString &errorStr, QDate firstDate, QDate lastDate);

agrolib/interpolation/interpolation.cpp

@@ -1081,7 +1081,7 @@ void localSelection(vector <Crit3DInterpolationDataPoint> &inputPoints, vector <
                     float x, float y, Crit3DInterpolationSettings& mySettings)
 {
     // search more stations to assure min points with all valid proxies
-    float ratioMinPoints = float(1.3);
+    float ratioMinPoints = float(1.2);
     unsigned minPoints = unsigned(mySettings.getMinPointsLocalDetrending() * ratioMinPoints);
     if (inputPoints.size() <= minPoints)
     {
@@ -1094,7 +1094,7 @@ void localSelection(vector <Crit3DInterpolationDataPoint> &inputPoints, vector <
         inputPoints[i].distance = gis::computeDistance(x, y, float((inputPoints[i]).point->utm.x), float((inputPoints[i]).point->utm.y));
 
     unsigned int nrValid = 0;
-    float stepRadius = 1000;           // [m]
+    float stepRadius = 7500;           // [m]
     float r0 = 0;                       // [m]
     float r1 = stepRadius;              // [m]
     unsigned int i;
@@ -1116,6 +1116,8 @@ void localSelection(vector <Crit3DInterpolationDataPoint> &inputPoints, vector <
                     nrPrimaries++;
             }
         }
+        if (nrValid > int(minPoints*0.8))
+            stepRadius = 1000;
         r0 = r1;
         r1 += stepRadius;
     }
@@ -2204,10 +2206,10 @@ float interpolate(vector <Crit3DInterpolationDataPoint> &myPoints, Crit3DInterpo
         myResult = modifiedShepardIdw(myPoints, mySettings, radius, myX, myY);
     }
 
-    if (int(myResult) != int(NODATA))
-        myResult += retrend(myVar, myProxyValues, mySettings);
-    else
+    if (int(myResult) == int(NODATA))
         return NODATA;
+    else if (!mySettings->getUseDoNotRetrend())
+        myResult += retrend(myVar, myProxyValues, mySettings);
 
     if (myVar == precipitation || myVar == dailyPrecipitation)
     {

agrolib/interpolation/interpolationSettings.cpp

@@ -510,6 +510,9 @@ bool Crit3DInterpolationSettings::getUseTD()
 bool Crit3DInterpolationSettings::getUseLocalDetrending()
 { return useLocalDetrending;}
 
+bool Crit3DInterpolationSettings::getUseDoNotRetrend()
+{ return useDoNotRetrend;}
+
 float Crit3DInterpolationSettings::getMaxHeightInversion()
 { return maxHeightInversion;}
 
@@ -528,6 +531,9 @@ void Crit3DInterpolationSettings::setUseTD(bool myValue)
 void Crit3DInterpolationSettings::setUseLocalDetrending(bool myValue)
 { useLocalDetrending = myValue;}
 
+void Crit3DInterpolationSettings::setUseDoNotRetrend(bool myValue)
+{ useDoNotRetrend = myValue;}
+
 void Crit3DInterpolationSettings::setUseDewPoint(bool myValue)
 { useDewPoint = myValue;}
 

agrolib/interpolation/interpolationSettings.h

@@ -150,6 +150,7 @@
         bool useMultipleDetrending;
         bool useDewPoint;
         bool useInterpolatedTForRH;
+        bool useDoNotRetrend;
         int minPointsLocalDetrending;
         bool meteoGridUpscaleFromDem;
         aggregationMethod meteoGridAggrMethod;
@@ -201,6 +202,9 @@
         void setUseLocalDetrending(bool myValue);
         bool getUseLocalDetrending();
 
+        void setUseDoNotRetrend(bool myValue);
+        bool getUseDoNotRetrend();
+
         void setUseDewPoint(bool myValue);
         bool getUseDewPoint();
 

agrolib/interpolation/spatialControl.cpp

@@ -125,7 +125,7 @@ bool computeResiduals(meteoVariable myVar, Crit3DMeteoPoint* meteoPoints, int nr
 
             if ((interpolatedValue != NODATA) && (myValue != NODATA))
             {
-                meteoPoints[i].residual = interpolatedValue - myValue;
+                meteoPoints[i].residual = myValue - interpolatedValue;
             }
         }
     }
@@ -191,7 +191,7 @@ bool computeResidualsLocalDetrending(meteoVariable myVar, Crit3DTime myTime, Cri
 
             if ((interpolatedValue != NODATA) && (myValue != NODATA))
             {
-                meteoPoints[i].residual = interpolatedValue - myValue;
+                meteoPoints[i].residual = myValue - interpolatedValue;
             }
         }
     }
@@ -213,7 +213,7 @@ float computeErrorCrossValidation(Crit3DMeteoPoint* myPoints, int nrMeteoPoints)
             if (value != NODATA && residual != NODATA)
             {
                 obsValues.push_back(value);
-                estValues.push_back(value + residual);
+                estValues.push_back(value - residual);
             }
         }
     }

agrolib/mathFunctions/statistics.cpp

@@ -242,7 +242,7 @@ namespace statistics
         return sigma;
     }
 
-    double compoundRelativeError(std::vector <float> measured, std::vector <float> simulated)
+    double NashSutcliffeEfficiency(std::vector <float> measured, std::vector <float> simulated)
     {
         if (measured.size() != simulated.size()) return NODATA;
 
@@ -264,7 +264,7 @@ namespace statistics
             if (!isEqual(measured[i], NODATA) && !isEqual(simulated[i], NODATA))
                 sumError += (simulated[i] - measured[i]) * (simulated[i] - measured[i]);
 
-        return sumError / sumDev;
+        return 1 - (sumError / sumDev);
     }
 
     float coefficientOfVariation(float *measured , float *simulated , int nrData)

agrolib/mathFunctions/statistics.h

@@ -27,7 +27,7 @@
         double rootMeanSquareError(std::vector <float> measured, std::vector <float> simulated);
         float meanError(std::vector<float> measured , std::vector<float> simulated );
         float meanAbsoluteError(std::vector <float> measured, std::vector <float> simulated);
-        double compoundRelativeError(std::vector <float> measured, std::vector <float> simulated);
+        double NashSutcliffeEfficiency(std::vector <float> measured, std::vector <float> simulated);
         float coefficientOfVariation(float *measured , float *simulated , int nrData);
         float weighedMean(float *data , float *weights, int nrData);
         float linearInterpolation(float x1, float y1, float x2, float y2, float xx);

agrolib/meteo/meteo.h

@@ -101,7 +101,7 @@
                     snowWaterEquivalent, snowFall, snowSurfaceTemperature, snowInternalEnergy, snowSurfaceEnergy,
                     snowAge, snowLiquidWaterContent, snowMelt, sensibleHeat, latentHeat,
                     dailyWaterTableDepth, leafAreaIndex,
-                    anomaly, elaboration, noMeteoTerrain,
+                    anomaly, elaboration, cvResidual, noMeteoTerrain,
                     noMeteoVar};
 
 

agrolib/pragaProject/pragaProject.cpp

@@ -1761,7 +1761,7 @@ bool PragaProject::averageSeriesOnZonesMeteoGrid(meteoVariable variable, meteoCo
     int infoStep = 0;
     if (showInfo)
     {
-        infoStep = setProgressBar("Aggregating data...", this->meteoGridDbHandler->gridStructure().header().nrRows);
+        infoStep = setProgressBar("Creating data array...", this->meteoGridDbHandler->gridStructure().header().nrRows);
     }
 
     Crit3DMeteoPoint* meteoPointTemp = new Crit3DMeteoPoint;
@@ -1803,8 +1803,16 @@ bool PragaProject::averageSeriesOnZonesMeteoGrid(meteoVariable variable, meteoCo
      int nrDays = int(startDate.daysTo(endDate) + 1);
      std::vector< std::vector<float> > dailyElabAggregation(nrDays, std::vector<float>(int(zoneGrid->maximum), NODATA));
 
+     if (showInfo)
+     {
+         infoStep = setProgressBar("Computing spatial elaborations...", nrDays);
+     }
+
      for (int day = 0; day < nrDays; day++)
      {
+         if (showInfo && (day % infoStep) == 0)
+             updateProgressBar(day);
+
          for (int zoneRow = 0; zoneRow < zoneGrid->header->nrRows; zoneRow++)
          {
              for (int zoneCol = 0; zoneCol < zoneGrid->header->nrCols; zoneCol++)
@@ -1886,8 +1894,10 @@ bool PragaProject::averageSeriesOnZonesMeteoGrid(meteoVariable variable, meteoCo
          }
      }
 
+     if (showInfo) closeProgressBar();
+
      // save dailyElabAggregation result into DB
-     if (showInfo) setProgressBar("Save data...", 0);
+     if (showInfo) setProgressBar("Saving data...", 0);
      if (! aggregationDbHandler->saveAggrData(int(zoneGrid->maximum), aggregationString, periodType,
                                              startDate, endDate, variable, dailyElabAggregation))
      {
@@ -2673,7 +2683,6 @@ bool PragaProject::interpolationCrossValidationPeriod(QDate dateIni, QDate dateF
 
     // order variables for derived computation
     std::string errString;
-    QString myError;
     int myHour;
     QDate myDate = dateIni;
     frequencyType myFreq = getVarFrequency(myVar);
@@ -2685,24 +2694,23 @@ bool PragaProject::interpolationCrossValidationPeriod(QDate dateIni, QDate dateF
             return false;
     }
 
-    // loading point data
-    logInfoGUI("Loading meteo points data from " + dateIni.addDays(-1).toString("yyyy-MM-dd") + " to " + dateFin.toString("yyyy-MM-dd"));
-    // load one day before (for transmissivity)
-    if (! loadMeteoPointsData(dateIni, dateFin, myFreq == hourly, myFreq == daily, false))
-        return false;
-
-    crossValidationStatistics stats;
     QFile file(filename);
     if (! file.open(QIODevice::WriteOnly | QIODevice::Text))
     {
-        myError = "Error writing to file " + filename;
+        errorString = "Error writing to file " + filename;
         return false;
     }
 
+    // loading point data
+    logInfoGUI("Loading meteo points data from " + dateIni.addDays(-1).toString("yyyy-MM-dd") + " to " + dateFin.toString("yyyy-MM-dd"));
+    // load one day before (for transmissivity)
+    if (! loadMeteoPointsData(dateIni, dateFin, myFreq == hourly, myFreq == daily, false))
+        return false;
+
     Crit3DTime myTime;
 
     QTextStream cvOutput(&file);
-    cvOutput << "Time,MAE,MBE,RMSE,CRE,R2" << '\n';
+    cvOutput << "Time,MAE,MBE,RMSE,NS,R2" << '\n';
 
     logInfoGUI("Cross validating " + QString::fromStdString(getMeteoVarName(myVar)) + " from " + dateIni.toString("yyyy-MM-dd") + " to " + dateFin.toString("yyyy-MM-dd"));
     while (myDate <= dateFin)
@@ -2715,29 +2723,29 @@ bool PragaProject::interpolationCrossValidationPeriod(QDate dateIni, QDate dateF
             {
                 myTime = getCrit3DTime(myDate, myHour);
 
-                if (interpolationCv(myVar, myTime, &stats))
+                if (interpolationCv(myVar, myTime))
                 {
                     cvOutput << getQDateTime(myTime).toString();
-                    cvOutput << "," << stats.getMeanAbsoluteError();
-                    cvOutput << "," << stats.getMeanBiasError();
-                    cvOutput << "," << stats.getRootMeanSquareError();
-                    cvOutput << "," << stats.getCompoundRelativeError();
-                    cvOutput << "," << stats.getR2() << '\n';
+                    cvOutput << "," << crossValidationStatistics.getMeanAbsoluteError();
+                    cvOutput << "," << crossValidationStatistics.getMeanBiasError();
+                    cvOutput << "," << crossValidationStatistics.getRootMeanSquareError();
+                    cvOutput << "," << crossValidationStatistics.getNashSutcliffeEfficiency();
+                    cvOutput << "," << crossValidationStatistics.getR2() << '\n';
                 }
             }
         }
         else
         {
             myTime = getCrit3DTime(myDate, 0);
 
-            if (interpolationCv(myVar, myTime, &stats))
+            if (interpolationCv(myVar, myTime))
             {
                 cvOutput << getQDateTime(myTime).date().toString();
-                cvOutput << "," << stats.getMeanAbsoluteError();
-                cvOutput << "," << stats.getMeanBiasError();
-                cvOutput << "," << stats.getRootMeanSquareError();
-                cvOutput << "," << stats.getCompoundRelativeError();
-                cvOutput << "," << stats.getR2() << '\n';
+                cvOutput << "," << crossValidationStatistics.getMeanAbsoluteError();
+                cvOutput << "," << crossValidationStatistics.getMeanBiasError();
+                cvOutput << "," << crossValidationStatistics.getRootMeanSquareError();
+                cvOutput << "," << crossValidationStatistics.getNashSutcliffeEfficiency();
+                cvOutput << "," << crossValidationStatistics.getR2() << '\n';
             }
         }