update interpolation output points

dbMeteoGrid/dbMeteoGrid.cpp

@@ -2374,7 +2374,7 @@ std::vector<float> Crit3DMeteoGridDbHandler::loadGridHourlyVar(QString *myError,
     QString tableH = _tableHourly.prefix + meteoPoint + _tableHourly.postFix;
     QDateTime dateTime, previousDateTime;
     dateTime.setTimeSpec(Qt::UTC);
-    previousDateTime.setTimeSpec((Qt::UTC));
+    previousDateTime.setTimeSpec(Qt::UTC);
 
     std::vector<float> hourlyVarList;
 

dbMeteoPoints/dbMeteoPointsHandler.cpp

@@ -224,6 +224,7 @@ QDateTime Crit3DMeteoPointsDbHandler::getFirstDate(frequencyType frequency)
         }
     }
 
+    QString tmp = firstDate.toString("yyyy-MM-dd");
     return firstDate;
 }
 

gis/gis.cpp

@@ -371,13 +371,15 @@ namespace gis
     }
 
 
+    // clean the grid (all NO DATA)
     void Crit3DRasterGrid::emptyGrid()
     {
         for (int row = 0; row < header->nrRows; row++)
             for (int col = 0; col < header->nrCols; col++)
                 value[row][col] = header->flag;
     }
 
+
     Crit3DRasterGrid::~Crit3DRasterGrid()
     {
         clear();

interpolation/interpolation.cpp

@@ -36,6 +36,7 @@
 #include "meteoPoint.h"
 #include "gis.h"
 #include "spatialControl.h"
+#include "interpolationPoint.h"
 #include "interpolation.h"
 #include <functional>
 
@@ -1048,25 +1049,18 @@ void localSelection(vector <Crit3DInterpolationDataPoint> &inputPoints, vector <
     mySettings.setLocalRadius(r1);
 }
 
-bool checkPrecipitationZero(std::vector <Crit3DInterpolationDataPoint> &myPoints, float precThreshold, int* nrPrecNotNull, bool* flatPrecipitation)
+
+bool checkPrecipitationZero(const std::vector<Crit3DInterpolationDataPoint> &myPoints, float precThreshold, int &nrNotNull)
 {
-    *flatPrecipitation = true;
-    *nrPrecNotNull = 0;
-    float myValue = NODATA;
+    nrNotNull = 0;
 
     for (unsigned int i = 0; i < myPoints.size(); i++)
         if (myPoints[i].isActive)
-            if (int(myPoints[i].value) != int(NODATA))
-                if (myPoints[i].value >= float(precThreshold))
-                {
-                    if (*nrPrecNotNull > 0 && myPoints[i].value != myValue)
-                        *flatPrecipitation = false;
+            if (! isEqual(myPoints[i].value, NODATA))
+                if (myPoints[i].value >= precThreshold)
+                    nrNotNull++;
 
-                    myValue = myPoints[i].value;
-                    (*nrPrecNotNull)++;
-                }
-
-    return (*nrPrecNotNull == 0);
+    return (nrNotNull == 0);
 }
 
 
@@ -1643,8 +1637,7 @@ bool preInterpolation(std::vector <Crit3DInterpolationDataPoint> &myPoints, Crit
     if (myVar == precipitation || myVar == dailyPrecipitation)
     {
         int nrPrecNotNull;
-        bool isFlatPrecipitation;
-        if (checkPrecipitationZero(myPoints, meteoSettings->getRainfallThreshold(), &nrPrecNotNull, &isFlatPrecipitation))
+        if (checkPrecipitationZero(myPoints, meteoSettings->getRainfallThreshold(), nrPrecNotNull))
         {
             mySettings->setPrecipitationAllZero(true);
             return true;
@@ -1727,15 +1720,16 @@ float interpolate(vector <Crit3DInterpolationDataPoint> &myPoints, Crit3DInterpo
         myResult = MAXVALUE(myResult, 0);
 
     return myResult;
-
 }
 
-bool getActiveProxyValues(Crit3DInterpolationSettings* mySettings, std::vector <double> & allProxyValues, std::vector <double> & activeProxyValues)
+
+bool getActiveProxyValues(Crit3DInterpolationSettings *mySettings, const std::vector<double> &allProxyValues, std::vector<double> &activeProxyValues)
 {
     Crit3DProxyCombination myCombination = mySettings->getCurrentCombination();
-    std::vector <double> myValues;
 
-    if (allProxyValues.size() != mySettings->getProxyNr()) return false;
+    if (allProxyValues.size() != mySettings->getProxyNr())
+        return false;
+
     activeProxyValues.clear();
 
     bool isComplete = true;
@@ -1744,12 +1738,14 @@ bool getActiveProxyValues(Crit3DInterpolationSettings* mySettings, std::vector <
         if (myCombination.getValue(i) && mySettings->getProxy(i)->getIsSignificant())
         {
             activeProxyValues.push_back(allProxyValues[i]);
-            if (allProxyValues[i] == NODATA) isComplete = false;
+            if (allProxyValues[i] == NODATA)
+                isComplete = false;
         }
 
     return isComplete;
 }
 
+
 void getProxyValuesXY(float x, float y, Crit3DInterpolationSettings* mySettings, std::vector<double> &myValues)
 {
     float myValue;

interpolation/interpolation.h

@@ -53,13 +53,15 @@
     bool krigLinearPrep(double *mySlope, double *myNugget, int nrPointData);
 
     void clearInterpolationPoints();
+    bool checkPrecipitationZero(const std::vector<Crit3DInterpolationDataPoint> &myPoints, float precThreshold, int &nrNotNull);
 
     bool neighbourhoodVariability(meteoVariable myVar, std::vector<Crit3DInterpolationDataPoint> &myInterpolationPoints, Crit3DInterpolationSettings *mySettings, float x, float y, float z, int nMax,
                                   float* devSt, float* avgDeltaZ, float* minDistance);
 
     float interpolate(std::vector<Crit3DInterpolationDataPoint> &myPoints, Crit3DInterpolationSettings *mySettings, Crit3DMeteoSettings *meteoSettings, meteoVariable myVar, float myX, float myY, float myZ, std::vector<double> myProxyValues, bool excludeSupplemental);
     void getProxyValuesXY(float x, float y, Crit3DInterpolationSettings* mySettings, std::vector<double> &myValues);
-    bool getActiveProxyValues(Crit3DInterpolationSettings* mySettings, std::vector<double> &allProxyValues, std::vector<double> &activeProxyValues);
+
+    bool getActiveProxyValues(Crit3DInterpolationSettings *mySettings, const std::vector<double> &allProxyValues, std::vector<double> &activeProxyValues);
 
     void detrending(std::vector <Crit3DInterpolationDataPoint> &myPoints,
                     Crit3DProxyCombination myCombination, Crit3DInterpolationSettings *mySettings, Crit3DClimateParameters *myClimate,

pragaProject/pragaMeteoMaps.cpp

@@ -38,6 +38,15 @@ void PragaHourlyMeteoMaps::clear()
 }
 
 
+void PragaHourlyMeteoMaps::initialize()
+{
+    mapHourlyWindVectorInt->emptyGrid();
+    mapHourlyWindVectorDir->emptyGrid();
+    mapHourlyWindVectorX->emptyGrid();
+    mapHourlyWindVectorY->emptyGrid();
+}
+
+
 gis::Crit3DRasterGrid* PragaHourlyMeteoMaps::getMapFromVar(meteoVariable myVar)
 {
     if (myVar == windVectorIntensity)

pragaProject/pragaMeteoMaps.h

@@ -22,6 +22,8 @@
         ~PragaHourlyMeteoMaps();
 
         void clear();
+        void initialize();
+
         gis::Crit3DRasterGrid* getMapFromVar(meteoVariable myVar);
         bool computeWindVector();
     };

pragaProject/pragaProject.cpp

@@ -1699,13 +1699,17 @@ QString getMapFileOutName(meteoVariable myVar, QDate myDate, int myHour)
     return name;
 }
 
+
 gis::Crit3DRasterGrid* PragaProject::getPragaMapFromVar(meteoVariable myVar)
 {
     gis::Crit3DRasterGrid* myGrid = nullptr;
 
+    // search in project maps (hourlyMeteoMaps and radiationMaps)
     myGrid = getHourlyMeteoRaster(myVar);
-    if (myGrid == nullptr) myGrid = pragaHourlyMaps->getMapFromVar(myVar);
+    // search in pragaDailyMaps
     if (myGrid == nullptr) myGrid = pragaDailyMaps->getMapFromVar(myVar);
+    // saerch in pragaHourlyMaps
+    if (myGrid == nullptr) myGrid = pragaHourlyMaps->getMapFromVar(myVar);
 
     return myGrid;
 }
@@ -1908,6 +1912,12 @@ bool PragaProject::interpolationOutputPointsPeriod(QDate firstDate, QDate lastDa
         return false;
     }
 
+    if (! DEM.isLoaded)
+    {
+        errorString = "Load a Digital Elevation Model before.";
+        return false;
+    }
+
     // check dates
     if (firstDate.isNull() || lastDate.isNull() || firstDate > lastDate)
     {
@@ -1944,8 +1954,34 @@ bool PragaProject::interpolationOutputPointsPeriod(QDate firstDate, QDate lastDa
             varToSave.push_back(windVectorIntensity);
     }
 
+    // initialize maps
+    if (pragaDailyMaps == nullptr) pragaDailyMaps = new Crit3DDailyMeteoMaps(DEM);
+    if (pragaHourlyMaps == nullptr) pragaHourlyMaps = new PragaHourlyMeteoMaps(DEM);
+
+    // check if topographic distance is needed
+    bool useTd = false;
+    if (interpolationSettings.getUseTD())
+    {
+        foreach (myVar, variables)
+        {
+            if (getUseTdVar(myVar))
+            {
+                useTd = true;
+                break;
+            }
+        }
+    }
+
+    // load topographic distance
+    if (useTd)
+    {
+        logInfoGUI("Loading topographic distance maps...");
+        if (! loadTopographicDistanceMaps(true, false))
+            return false;
+    }
+
     int nrDays = firstDate.daysTo(lastDate) + 1;
-    int nrDaysLoading = std::min(nrDays, 10);
+    int nrDaysLoading = std::min(nrDays, 30);
     int countDaysSaving = 0;
 
     QDate myDate = firstDate;
@@ -1966,6 +2002,57 @@ bool PragaProject::interpolationOutputPointsPeriod(QDate firstDate, QDate lastDa
             if (! loadMeteoPointsData(myDate.addDays(-1), lastLoadingDate, isHourly, isDaily, false))
                 return false;
         }
+
+        if (isHourly)
+        {
+            // initialize
+            hourlyMeteoMaps->initialize();
+            radiationMaps->initialize();
+            pragaHourlyMaps->initialize();
+
+            for (int hour = 1; hour <= 24; hour++)
+            {
+                logInfoGUI("Interpolating hourly variables for " + myDate.toString("yyyy-MM-dd") + " " + QString::number(hour) + ":00");
+                foreach (myVar, variables)
+                {
+                    if (getVarFrequency(myVar) == hourly)
+                    {
+                        // TODO RH and other special variables
+
+                        if (! interpolationDemMain(myVar, getCrit3DTime(myDate, hour), getPragaMapFromVar(myVar)))
+                            return false;
+                    }
+                }
+            }
+        }
+
+        if (isDaily)
+        {
+            logInfoGUI("Interpolating daily variables for " + myDate.toString("yyyy-MM-dd"));
+
+            // initialize
+            pragaDailyMaps->initialize();
+
+            foreach (myVar, variables)
+            {
+                if (getVarFrequency(myVar) == daily)
+                {
+                    // TODO special variables
+
+                    if (! interpolationDemMain(myVar, getCrit3DTime(myDate, 0), getPragaMapFromVar(myVar)))
+                        return false;
+
+                    // fix daily temperatures consistency
+                    if (myVar == dailyAirTemperatureMax || myVar == dailyAirTemperatureMin)
+                    {
+                        if (! pragaDailyMaps->fixDailyThermalConsistency())
+                            return false;
+                    }
+                }
+            }
+        }
+
+        myDate = myDate.addDays(1);
     }
 
     closeLogInfo();
@@ -2214,7 +2301,6 @@ bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QL
                     }
 
                     meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, daily, getCrit3DDate(myDate), 0, 0, &DEM, myGrid, interpolationSettings.getMeteoGridAggrMethod());
-
                 }
             }
         }
@@ -2246,9 +2332,9 @@ bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QL
         return false;
 
     return true;
-
 }
 
+
 bool PragaProject::interpolationMeteoGrid(meteoVariable myVar, frequencyType myFrequency, const Crit3DTime& myTime)
 {
     if (meteoGridDbHandler != nullptr)

project/meteoMaps.cpp

@@ -58,6 +58,20 @@ void Crit3DDailyMeteoMaps::clear()
 }
 
 
+void Crit3DDailyMeteoMaps::initialize()
+{
+    mapDailyTAvg->emptyGrid();
+    mapDailyTMax->emptyGrid();
+    mapDailyTMin->emptyGrid();
+    mapDailyPrec->emptyGrid();
+    mapDailyRHAvg->emptyGrid();
+    mapDailyRHMin->emptyGrid();
+    mapDailyRHMax->emptyGrid();
+    mapDailyLeafW->emptyGrid();
+    mapDailyET0HS->emptyGrid();
+}
+
+
 bool Crit3DDailyMeteoMaps::computeHSET0Map(gis::Crit3DGisSettings* gisSettings, Crit3DDate myDate)
 {
     float airTmin, airTmax;
@@ -84,30 +98,26 @@ bool Crit3DDailyMeteoMaps::computeHSET0Map(gis::Crit3DGisSettings* gisSettings,
     return gis::updateMinMaxRasterGrid(mapDailyET0HS);
 }
 
+
 bool Crit3DDailyMeteoMaps::fixDailyThermalConsistency()
 {
     if (! mapDailyTMax->isLoaded || ! mapDailyTMin->isLoaded) return true;
     if ( mapDailyTMax->getMapTime() != mapDailyTMin->getMapTime()) return true;
 
-    float TRange = NODATA;
-    unsigned row, col;
-
-    for (row = 0; row < unsigned(mapDailyTMax->header->nrRows); row++)
-        for (col = 0; col < unsigned(mapDailyTMax->header->nrCols); col++)
+    for (unsigned row = 0; row < unsigned(mapDailyTMax->header->nrRows); row++)
+    {
+        for (unsigned col = 0; col < unsigned(mapDailyTMax->header->nrCols); col++)
         {
             if (! isEqual(mapDailyTMax->value[row][col], mapDailyTMax->header->flag) &&
                 ! isEqual(mapDailyTMin->value[row][col], mapDailyTMin->header->flag))
             {
                 if (mapDailyTMin->value[row][col] > mapDailyTMax->value[row][col])
                 {
-                    //TRange = findNeighbourTRangeRaster(grdTmax, grdTmin, myRow, myCol)
-                    if (! isEqual(TRange, NODATA))
-                        mapDailyTMin->value[row][col] = mapDailyTMax->value[row][col] - TRange;
-                    else
-                        mapDailyTMin->value[row][col] = mapDailyTMax->value[row][col] - 0.1f;
+                    mapDailyTMin->value[row][col] = mapDailyTMax->value[row][col] - 0.1f;
                 }
             }
         }
+    }
 
     return true;
 }

project/meteoMaps.h

@@ -32,8 +32,10 @@
         ~Crit3DDailyMeteoMaps();
 
         void clear();
+        void initialize();
 
         gis::Crit3DRasterGrid* getMapFromVar(meteoVariable myVar);
+
         bool computeHSET0Map(gis::Crit3DGisSettings *gisSettings, Crit3DDate myDate);
         bool fixDailyThermalConsistency();
     };

project/project.cpp

@@ -63,6 +63,9 @@ void Project::initializeProject()
     outputPointsDbHandler = nullptr;
     meteoGridDbHandler = nullptr;
     aggregationDbHandler = nullptr;
+
+    meteoPointsDbFirstTime.setTimeSpec(Qt::UTC);
+    meteoPointsDbLastTime.setTimeSpec(Qt::UTC);
     meteoPointsDbFirstTime.setSecsSinceEpoch(0);
     meteoPointsDbLastTime.setSecsSinceEpoch(0);
 
@@ -2228,6 +2231,7 @@ bool Project::interpolationDem(meteoVariable myVar, const Crit3DTime& myTime, gi
     }
 
     myRaster->setMapTime(myTime);
+
     return true;
 }