diff --git a/agrolib/commonChartElements/dialogChangeAxis.cpp b/agrolib/commonChartElements/dialogChangeAxis.cpp index 6fd91e0b..c14f64fe 100644 --- a/agrolib/commonChartElements/dialogChangeAxis.cpp +++ b/agrolib/commonChartElements/dialogChangeAxis.cpp @@ -1,18 +1,24 @@ #include "dialogChangeAxis.h" -DialogChangeAxis::DialogChangeAxis(bool isLeftAxis_) +DialogChangeAxis::DialogChangeAxis(int nrAxis, bool isDate_) { - isLeftAxis = isLeftAxis_; + isDateAxis = isDate_; + QString title; - if (isLeftAxis) + if (nrAxis == 0) + { + title = "Change X Axis"; + } + else if (nrAxis == 1) { title = "Change Left Axis"; } - else + else if (nrAxis == 2) { title = "Change Right Axis"; } this->setWindowTitle(title); + QVBoxLayout* mainLayout = new QVBoxLayout; this->resize(200, 100); @@ -20,17 +26,32 @@ DialogChangeAxis::DialogChangeAxis(bool isLeftAxis_) QHBoxLayout *layoutEdit = new QHBoxLayout; QLabel minValueLabel("Minimum value:"); - minValueLabel.setBuddy(&minVal); - minVal.setValidator(new QDoubleValidator(-999.0, 999.0, 3)); + layoutEdit->addWidget(&minValueLabel); + if (isDateAxis) + { + minValueLabel.setBuddy(&minDate); + layoutEdit->addWidget(&minDate); + } + else + { + minValueLabel.setBuddy(&minVal); + minVal.setValidator(new QDoubleValidator(-9999.0, 9999.0, 3)); + layoutEdit->addWidget(&minVal); + } QLabel maxValueLabel("Maximum value:"); - maxValueLabel.setBuddy(&maxVal); - maxVal.setValidator(new QDoubleValidator(-999.0, 999.0, 3)); - - layoutEdit->addWidget(&minValueLabel); - layoutEdit->addWidget(&minVal); layoutEdit->addWidget(&maxValueLabel); - layoutEdit->addWidget(&maxVal); + if (isDateAxis) + { + minValueLabel.setBuddy(&maxDate); + layoutEdit->addWidget(&maxDate); + } + else + { + maxValueLabel.setBuddy(&maxVal); + maxVal.setValidator(new QDoubleValidator(-9999.0, 9999.0, 3)); + layoutEdit->addWidget(&maxVal); + } QDialogButtonBox buttonBox(QDialogButtonBox::Ok | QDialogButtonBox::Cancel); @@ -55,16 +76,28 @@ void DialogChangeAxis::done(bool res) { if (res) // ok { - if (minVal.text().isEmpty()) + if (isDateAxis) { - QMessageBox::information(nullptr, "Missing min value", "Insert min val"); - return; + if (minDate.date() >= maxDate.date()) + { + QMessageBox::warning(nullptr, "Wrong dates!", "Insert correct dates."); + return; + } } - if (maxVal.text().isEmpty()) + else { - QMessageBox::information(nullptr, "Missing max value", "Insert max val"); - return; + if (minVal.text().isEmpty()) + { + QMessageBox::warning(nullptr, "Missing min value", "Insert minimum value."); + return; + } + if (maxVal.text().isEmpty()) + { + QMessageBox::warning(nullptr, "Missing max value", "Insert maximum value."); + return; + } } + QDialog::done(QDialog::Accepted); return; } @@ -84,3 +117,13 @@ float DialogChangeAxis::getMaxVal() const { return maxVal.text().toFloat(); } + +QDate DialogChangeAxis::getMinDate() const +{ + return minDate.date(); +} + +QDate DialogChangeAxis::getMaxDate() const +{ + return maxDate.date(); +} diff --git a/agrolib/commonChartElements/dialogChangeAxis.h b/agrolib/commonChartElements/dialogChangeAxis.h index a6339371..bf230494 100644 --- a/agrolib/commonChartElements/dialogChangeAxis.h +++ b/agrolib/commonChartElements/dialogChangeAxis.h @@ -8,16 +8,24 @@ Q_OBJECT private: - bool isLeftAxis; + bool isDateAxis; + + QDateEdit minDate; QLineEdit minVal; + + QDateEdit maxDate; QLineEdit maxVal; public: - DialogChangeAxis(bool isLeftAxis); + DialogChangeAxis(int nrAxis, bool isDate_); ~DialogChangeAxis() override; void done(bool res); + float getMinVal() const; float getMaxVal() const; + + QDate getMinDate() const; + QDate getMaxDate() const; }; #endif // DIALOGCHANGEAXIS_H diff --git a/agrolib/crop/crop.cpp b/agrolib/crop/crop.cpp index d89405dc..5bdd9fb4 100644 --- a/agrolib/crop/crop.cpp +++ b/agrolib/crop/crop.cpp @@ -49,7 +49,7 @@ Crit3DCrop::Crit3DCrop() void Crit3DCrop::clear() { idCrop = ""; - type = HERBACEOUS_ANNUAL; + type = BARESOIL; roots.clear(); diff --git a/agrolib/crop/landUnit.cpp b/agrolib/crop/landUnit.cpp index be279418..d39f754d 100644 --- a/agrolib/crop/landUnit.cpp +++ b/agrolib/crop/landUnit.cpp @@ -59,13 +59,3 @@ bool loadLandUnitList(const QSqlDatabase &dbCrop, std::vector &l return true; } - - -int getLandUnitIndex(const std::vector &landUnitList, int idLandUnit) -{ - for (int index = 0; index < int(landUnitList.size()); index++) - if (landUnitList[index].id == idLandUnit) - return index; - - return NODATA; -} diff --git a/agrolib/crop/landUnit.h b/agrolib/crop/landUnit.h index be563b3d..fb11f46c 100644 --- a/agrolib/crop/landUnit.h +++ b/agrolib/crop/landUnit.h @@ -23,7 +23,5 @@ bool loadLandUnitList(const QSqlDatabase &dbCrop, std::vector &landUnitList, QString &errorStr); - int getLandUnitIndex(const std::vector &landUnitList, int idLandUnit); - #endif // LANDUNIT_H diff --git a/agrolib/dbMeteoGrid/dbMeteoGrid.cpp b/agrolib/dbMeteoGrid/dbMeteoGrid.cpp index a39a9b2d..22f804c7 100644 --- a/agrolib/dbMeteoGrid/dbMeteoGrid.cpp +++ b/agrolib/dbMeteoGrid/dbMeteoGrid.cpp @@ -2501,7 +2501,7 @@ std::vector Crit3DMeteoGridDbHandler::loadGridDailyVar(QString *errorStr, // read last date qry.last(); - if (!getValue(qry.value(_tableDaily.fieldTime), &lastDateDB)) + if (! getValue(qry.value(_tableDaily.fieldTime), &lastDateDB)) { *errorStr = "Missing last date"; return dailyVarList; diff --git a/agrolib/dbMeteoPoints/dbAggregationsHandler.cpp b/agrolib/dbMeteoPoints/dbAggregationsHandler.cpp index 1b343b9f..731670b4 100644 --- a/agrolib/dbMeteoPoints/dbAggregationsHandler.cpp +++ b/agrolib/dbMeteoPoints/dbAggregationsHandler.cpp @@ -243,7 +243,6 @@ bool Crit3DAggregationsDbHandler::writePointProperties(int numZones, QString agg QSqlQuery qry(_db); for (int i = 1; i <= numZones; i++) - { QString id = QString::number(i) + "_" + aggrType; QString name = id; diff --git a/agrolib/gis/gis.cpp b/agrolib/gis/gis.cpp index eb099139..7bd2265d 100644 --- a/agrolib/gis/gis.cpp +++ b/agrolib/gis/gis.cpp @@ -287,13 +287,13 @@ namespace gis bool Crit3DRasterGrid::initializeParameters(const Crit3DRasterHeader& initHeader) { - parametersCell.clear(); - parametersCell.resize(initHeader.nrRows*initHeader.nrCols); - for (int i = 0; i < int(parametersCell.size()); i++) + singleCell.clear(); + singleCell.resize(initHeader.nrRows*initHeader.nrCols); + for (int i = 0; i < int(singleCell.size()); i++) { - parametersCell[i].row = i / initHeader.nrCols; - parametersCell[i].col = i % initHeader.nrCols; - parametersCell[i].fittingParameters.clear(); + singleCell[i].row = i / initHeader.nrCols; + singleCell[i].col = i % initHeader.nrCols; + singleCell[i].fittingParameters.clear(); } return true; @@ -301,12 +301,12 @@ namespace gis bool Crit3DRasterGrid::initializeParametersLatLonHeader(const Crit3DLatLonHeader& latLonHeader) { - parametersCell.resize(latLonHeader.nrRows*latLonHeader.nrCols); - for (int i = 0; i < int(parametersCell.size()); i++) + singleCell.resize(latLonHeader.nrRows*latLonHeader.nrCols); + for (int i = 0; i < int(singleCell.size()); i++) { - parametersCell[i].row = i / latLonHeader.nrCols; - parametersCell[i].col = i % latLonHeader.nrCols; - parametersCell[i].fittingParameters.clear(); + singleCell[i].row = i / latLonHeader.nrCols; + singleCell[i].col = i % latLonHeader.nrCols; + singleCell[i].fittingParameters.clear(); } return true; @@ -497,8 +497,8 @@ namespace gis int index = row * header->nrCols + col; - if (index < int(parametersCell.size())) - parameters = parametersCell[index].fittingParameters; + if (index < int(singleCell.size())) + parameters = singleCell[index].fittingParameters; return parameters; @@ -511,7 +511,7 @@ namespace gis return false; int index = row * header->nrCols + col; - parametersCell[index].fittingParameters = parameters; + singleCell[index].fittingParameters = parameters; return true; } @@ -542,7 +542,7 @@ namespace gis for (m = col-1; m < col+2; m++) { index = l * header->nrCols + m; - if (index >= 0 && index < int(parametersCell.size()) && (parametersCell[index].fittingParameters.size() > i && !parametersCell[index].fittingParameters[i].empty()) && (l != row || m !=col)) { + if (index >= 0 && index < int(singleCell.size()) && (singleCell[index].fittingParameters.size() > i && !singleCell[index].fittingParameters[i].empty()) && (l != row || m !=col)) { findFirst = 1; } if (findFirst==1) break; @@ -555,7 +555,7 @@ namespace gis //you're on a specific proxy rn. cycle through the cells, calculate the avg avg.clear(); - avg.resize(parametersCell[index].fittingParameters[i].size()); + avg.resize(singleCell[index].fittingParameters[i].size()); counter = 0; for (k = l; k < row+2; k++) @@ -563,10 +563,10 @@ namespace gis for (p = m; p < col+2; p++) { index = k * header->nrCols + p; - if (index >= 0 && index < int(parametersCell.size()) && parametersCell[index].fittingParameters.size() > i && !parametersCell[index].fittingParameters[i].empty()) { + if (index >= 0 && index < int(singleCell.size()) && singleCell[index].fittingParameters.size() > i && !singleCell[index].fittingParameters[i].empty()) { for (unsigned int o = 0; o < avg.size(); o++) { - avg[o] += parametersCell[index].fittingParameters[i][o]; + avg[o] += singleCell[index].fittingParameters[i][o]; } counter++; @@ -583,63 +583,6 @@ namespace gis return tempProxyVector; } - std::vector> Crit3DRasterGrid::prepareParametersOld(int row, int col, unsigned int activeProxyNr) - { - std::vector> tempProxyVector; - std::vector tempParVector; - tempProxyVector.clear(); - tempParVector.clear(); - double avg; - int counter; - int l, i, j, k; - int index; - bool findFirst = 0; - - if (isOutOfGrid(row, col)) - return tempProxyVector; - - //look for the first cell that has data in it. if there isn't any, return empty vector - for (i = row-1; i < row+2; i++) - { - for (j = col-1; j < col+2; j++) - { - index = i * header->nrCols + j; - if (index >= 0 && index < int(parametersCell.size()) && (parametersCell[index].fittingParameters.size() == activeProxyNr) && (i != row || j !=col)) - findFirst = 1; - if (findFirst==1) break; - } - if (findFirst==1) break; - } - - for (k = 0; k < int(parametersCell[index].fittingParameters.size()); k++) - { - for (l = 0; l < int(parametersCell[index].fittingParameters[k].size()); l++) - { - avg = 0; - counter = 0; - for (int h = i; h < row+2; h++) - { - for (int m = j; m < col+2; m++) - { - index = h * header->nrCols + m; - if (index >= 0 && index < int(parametersCell.size()) && (parametersCell[index].fittingParameters.size() == activeProxyNr) && (i != row || j !=col)) - { - avg += parametersCell[index].fittingParameters[k][l]; - counter++; - } - } - } - tempParVector.push_back(avg/counter); - index = i*header->nrCols+j; - } - tempProxyVector.push_back(tempParVector); - tempParVector.clear(); - } - - return tempProxyVector; - - } - void convertFlagToNodata(Crit3DRasterGrid& myGrid) { if (myGrid.header->flag == NODATA) diff --git a/agrolib/gis/gis.h b/agrolib/gis/gis.h index 9130802e..9527dd51 100644 --- a/agrolib/gis/gis.h +++ b/agrolib/gis/gis.h @@ -139,7 +139,7 @@ Crit3DRasterCell(); }; - struct RasterGridParameters { + struct RasterGridCell { public: int row; int col; @@ -155,7 +155,7 @@ float minimum, maximum; bool isLoaded; Crit3DTime mapTime; - std::vector parametersCell; + std::vector singleCell; Crit3DUtmPoint* utmPoint(int myRow, int myCol); void getXY(int myRow, int myCol, double &x, double &y) const; @@ -192,7 +192,6 @@ std::vector> getParametersFromRowCol(int row, int col); bool setParametersForRowCol(int row, int col, std::vector> parameters); std::vector> prepareParameters(int row, int col, std::vector activeList); - std::vector> prepareParametersOld(int row, int col, unsigned int activeProxyNr); Crit3DTime getMapTime() const; void setMapTime(const Crit3DTime &value); @@ -255,8 +254,8 @@ bool openRaster(std::string fileName, Crit3DRasterGrid *rasterGrid, int currentUtmZone, std::string &errorStr); - bool readEsriGrid(std::string fileName, Crit3DRasterGrid* rasterGrid, std::string &errorStr); - bool writeEsriGrid(std::string fileName, Crit3DRasterGrid *rasterGrid, std::string &errorStr); + bool readEsriGrid(const std::string &fileName, Crit3DRasterGrid* rasterGrid, std::string &errorStr); + bool writeEsriGrid(const std::string &fileName, Crit3DRasterGrid *rasterGrid, std::string &errorStr); bool readEnviGrid(std::string fileName, Crit3DRasterGrid* rasterGrid, int currentUtmZone, std::string &errorStr); bool writeEnviGrid(std::string fileName, int utmZone, Crit3DRasterGrid *rasterGrid, std::string &errorStr); diff --git a/agrolib/gis/gisIO.cpp b/agrolib/gis/gisIO.cpp index c4e82984..62949938 100644 --- a/agrolib/gis/gisIO.cpp +++ b/agrolib/gis/gisIO.cpp @@ -112,17 +112,17 @@ namespace gis * \param error string * \return true on success, false otherwise */ - bool readEsriGridHeader(string fileName, gis::Crit3DRasterHeader *header, string &errorStr) + bool readEsriGridHeader(const string &fileName, gis::Crit3DRasterHeader *header, string &errorStr) { string myLine, myKey, upKey, valueStr; int nrKeys = 0; - fileName += ".hdr"; - ifstream myFile (fileName.c_str()); + string myFileName = fileName + ".hdr"; + ifstream myFile (myFileName.c_str()); - if (!myFile.is_open()) + if (! myFile.is_open()) { - errorStr = "Missing file: " + fileName; + errorStr = "Missing file: " + myFileName; return false; } @@ -315,7 +315,7 @@ namespace gis * \param error string * \return true on success, false otherwise */ - bool readRasterFloatData(string fileName, gis::Crit3DRasterGrid *rasterGrid, string &error) + bool readRasterFloatData(const string &fileName, gis::Crit3DRasterGrid *rasterGrid, string &error) { FILE* filePointer; @@ -382,12 +382,12 @@ namespace gis * \param error string pointer * \return true on success, false otherwise */ - bool writeEsriGridHeader(string myFileName, gis::Crit3DRasterHeader *myHeader, string &error) + bool writeEsriGridHeader(const string &fileName, gis::Crit3DRasterHeader *myHeader, string &error) { - myFileName += ".hdr"; + string myFileName = fileName + ".hdr"; ofstream myFile (myFileName.c_str()); - if (!myFile.is_open()) + if (! myFile.is_open()) { error = "File .hdr error."; return false; @@ -429,9 +429,9 @@ namespace gis * \param error string pointer * \return true on success, false otherwise */ - bool writeEsriGridFlt(string myFileName, gis::Crit3DRasterGrid* myGrid, string &error) + bool writeEsriGridFlt(const string &fileName, gis::Crit3DRasterGrid* myGrid, string &error) { - myFileName += ".flt"; + string myFileName = fileName + ".flt"; FILE* filePointer; @@ -454,13 +454,15 @@ namespace gis * \brief Write a ESRI float raster (.hdr and .flt) * \return true on success, false otherwise */ - bool writeEsriGrid(string fileName, Crit3DRasterGrid* rasterGrid, string &error) + bool writeEsriGrid(const string &fileName, Crit3DRasterGrid* rasterGrid, string &errorStr) { - if (gis::writeEsriGridHeader(fileName, rasterGrid->header, error)) - if (gis::writeEsriGridFlt(fileName, rasterGrid, error)) - return true; + if (! gis::writeEsriGridHeader(fileName, rasterGrid->header, errorStr)) + return false; - return false; + if (! gis::writeEsriGridFlt(fileName, rasterGrid, errorStr)) + return false; + + return true; } @@ -468,15 +470,15 @@ namespace gis * \brief Read a ESRI float raster (.hdr and .flt) * \return true on success, false otherwise */ - bool readEsriGrid(string fileName, Crit3DRasterGrid* rasterGrid, string &errorStr) + bool readEsriGrid(const string &fileName, Crit3DRasterGrid* rasterGrid, string &errorStr) { if (rasterGrid == nullptr) return false; rasterGrid->clear(); if(gis::readEsriGridHeader(fileName, rasterGrid->header, errorStr)) { - fileName += ".flt"; - if (gis::readRasterFloatData(fileName, rasterGrid, errorStr)) + string fltFileName = fileName + ".flt"; + if (gis::readRasterFloatData(fltFileName, rasterGrid, errorStr)) { gis::updateMinMaxRasterGrid(rasterGrid); rasterGrid->isLoaded = true; diff --git a/agrolib/grapevine/grapevine.cpp b/agrolib/grapevine/grapevine.cpp index 22125020..83e04c5d 100644 --- a/agrolib/grapevine/grapevine.cpp +++ b/agrolib/grapevine/grapevine.cpp @@ -21,9 +21,7 @@ const float LAIMIN = 0.01f; Vine3D_Grapevine::Vine3D_Grapevine() -{ -} - +{} bool Vine3D_Grapevine::compute(bool computeDaily, int secondsPerStep, Crit3DModelCase* modelCase, double chlorophyll) { @@ -106,9 +104,9 @@ bool Vine3D_Grapevine::compute(bool computeDaily, int secondsPerStep, Crit3DMode grassTranspiration(modelCase); return(true); - } + void Vine3D_Grapevine::resetLayers() { for (int i=0 ; i < nrMaxLayers ; i++) @@ -121,6 +119,7 @@ void Vine3D_Grapevine::resetLayers() } } + void Vine3D_Grapevine::initializeLayers(int myMaxLayers) { nrMaxLayers = myMaxLayers; @@ -156,22 +155,25 @@ bool Vine3D_Grapevine::setWeather(double meanDailyTemp, double temp, double irra myMeanDailyTemperature = meanDailyTemp; double deltaRelHum = MAXVALUE(100.0 - myRelativeHumidity, 0.01); myVaporPressureDeficit = 0.01 * deltaRelHum * 613.75 * exp(17.502 * myInstantTemp / (240.97 + myInstantTemp)); - //globalRadiation = globRad; - if ((int(prec) != NODATA) && (int(temp) != NODATA) && (int(windSpeed) != NODATA) && (int(irradiance) != NODATA) && (int(relativeHumidity) != NODATA) && (int(atmosphericPressure) != NODATA)) isReadingOK = true ; - return isReadingOK ; + + if ((int(prec) != NODATA) && (int(temp) != NODATA) && (int(windSpeed) != NODATA) + && (int(irradiance) != NODATA) && (int(relativeHumidity) != NODATA) && (int(atmosphericPressure) != NODATA)) + isReadingOK = true; + + return isReadingOK; } + bool Vine3D_Grapevine::setDerivedVariables(double diffuseIrradiance, double directIrradiance, double cloudIndex, double sunElevation) { bool isReadingOK = false; myDiffuseIrradiance = diffuseIrradiance ; myDirectIrradiance = directIrradiance ; - //myLongWaveIrradiance = longWaveIrradiance ; myCloudiness = MINVALUE(1,MAXVALUE(0,cloudIndex)) ; - //myAirVapourPressure = airVapourPressure ; mySunElevation = sunElevation; - if (int(sunElevation) != NODATA && int(diffuseIrradiance) != NODATA && int(directIrradiance) != NODATA - && int(cloudIndex) != NODATA) isReadingOK = true ; + if (int(sunElevation) != NODATA && int(diffuseIrradiance) != NODATA && int(directIrradiance) != NODATA && int(cloudIndex) != NODATA) + isReadingOK = true ; + return isReadingOK ; } @@ -185,32 +187,45 @@ void Vine3D_Grapevine::initializeWaterStress(Crit3DModelCase* modelCase) } -bool Vine3D_Grapevine::setSoilProfile(Crit3DModelCase* modelCase, double* myWiltingPoint, double* myFieldCapacity, - double* myPsiSoilProfile, double* mySoilWaterContentProfile, - double* mySoilWaterContentFC, double* mySoilWaterContentWP) +bool Vine3D_Grapevine::setSoilProfile(Crit3DModelCase* modelCase, std::vector& myWiltingPoint, std::vector& myFieldCapacity, + std::vector& myPsiSoilProfile, std::vector& mySoilWaterContentProfile, + std::vector& mySoilWaterContentFC, std::vector& mySoilWaterContentWP) { - double psiSoilProfile; - double logPsiSoilAverage = 0.; - double logPsiFCAverage = 0.; - double soilFieldCapacity; + // check last Soil layer + int lastLayer = modelCase->soilLayersNr - 1; + for (int i = 1; i < modelCase->soilLayersNr; i++) + { + if (isEqual(myPsiSoilProfile[i], NODATA) || isEqual(myFieldCapacity[i], NODATA) || isEqual(mySoilWaterContentProfile[i], NODATA) + || isEqual(mySoilWaterContentProfile[i], NODATA) || isEqual(mySoilWaterContentWP[i], NODATA)) + { + lastLayer = i-1; + break; + } + } + if (lastLayer == 0) + return false; psiSoilAverage = 0.; psiFieldCapacityAverage = 0.; - if (int(myWiltingPoint[int(modelCase->soilLayersNr / 2)]) == NODATA) - return false; + // initialize vector + for (int i = 1; i < modelCase->soilLayersNr; i++) + { + transpirationLayer[i] = 0.; + transpirationCumulatedGrass[i] = 0.; + fractionTranspirableSoilWaterProfile[i] = 0.; + } - wiltingPoint = myWiltingPoint[int(modelCase->soilLayersNr / 2)] / 101.97; // conversion from mH2O to MPa + int midLayer = (lastLayer + 1) / 2; + wiltingPoint = myWiltingPoint[midLayer] / 101.97; // conversion from mH2O to MPa //layer 0: surface, no soil - for (int i = 1; i < modelCase->soilLayersNr; i++) + double logPsiSoilAverage = 0.; + double logPsiFCAverage = 0.; + for (int i = 1; i <= lastLayer; i++) { - if (isEqual(myPsiSoilProfile[i], NODATA) || isEqual(myFieldCapacity[i], NODATA) || isEqual(mySoilWaterContentProfile[i], NODATA) - || isEqual(mySoilWaterContentProfile[i], NODATA) || isEqual(mySoilWaterContentWP[i], NODATA)) - return false; - - soilFieldCapacity = myFieldCapacity[i]/101.97; // conversion from mH2O to MPa - psiSoilProfile = MINVALUE(myPsiSoilProfile[i],-1.)/101.97 ; // conversion from mH2O to MPa + double soilFieldCapacity = myFieldCapacity[i] / 101.97; // conversion from mH2O to MPa + double psiSoilProfile = MINVALUE(myPsiSoilProfile[i], -1.) / 101.97 ; // conversion from mH2O to MPa logPsiSoilAverage += log(-psiSoilProfile) * modelCase->rootDensity[i]; logPsiFCAverage += log(-soilFieldCapacity) * modelCase->rootDensity[i]; } @@ -219,28 +234,26 @@ bool Vine3D_Grapevine::setSoilProfile(Crit3DModelCase* modelCase, double* myWilt psiFieldCapacityAverage = -exp(logPsiFCAverage); fractionTranspirableSoilWaterAverage = 0; - double waterContent, waterContentFC, waterContentWP; - - for (int i = 0; i < modelCase->soilLayersNr; i++) + for (int i = 0; i <= lastLayer; i++) { - waterContent = mySoilWaterContentProfile[i]; - waterContentFC = mySoilWaterContentFC[i]; - waterContentWP = mySoilWaterContentWP[i]; + double waterContent = mySoilWaterContentProfile[i]; + double waterContentFC = mySoilWaterContentFC[i]; + double waterContentWP = mySoilWaterContentWP[i]; fractionTranspirableSoilWaterProfile[i] = MAXVALUE(0, MINVALUE(1, (waterContent - waterContentWP) / (waterContentFC - waterContentWP))); fractionTranspirableSoilWaterAverage += fractionTranspirableSoilWaterProfile[i] * modelCase->rootDensity[i]; - transpirationLayer[i] = 0.; - transpirationCumulatedGrass[i] = 0. ; } + return true ; } -bool Vine3D_Grapevine::setStatePlant(TstatePlant myStatePlant, bool isVineyard) + +bool Vine3D_Grapevine::setStatePlant(TstatePlant myStatePlant, bool isVineyard_) { statePlant = myStatePlant; this->statePlant.outputPlant.transpirationNoStress = 0.; - if (! isVineyard) + if (! isVineyard_) { statePlant.outputPlant.brixBerry = NODATA; statePlant.statePheno.stage = NODATA; @@ -248,17 +261,8 @@ bool Vine3D_Grapevine::setStatePlant(TstatePlant myStatePlant, bool isVineyard) statePlant.stateGrowth.leafAreaIndex = NODATA; statePlant.stateGrowth.fruitBiomass = NODATA; } - return true; -} - -TstatePlant Vine3D_Grapevine::getStatePlant() -{ - return(this->statePlant); -} -ToutputPlant Vine3D_Grapevine::getOutputPlant() -{ - return (this->statePlant.outputPlant); + return true; } @@ -268,13 +272,13 @@ void Vine3D_Grapevine::getFixSimulationParameters() parameterBindiMigliettaFix.a = -0.28; parameterBindiMigliettaFix.b = 0.04; parameterBindiMigliettaFix.c = -0.015; - //parameterBindiMigliettaFix.baseTemperature = 10; - //parameterBindiMigliettaFix.tempMaxThreshold = 35; parameterBindiMigliettaFix.extinctionCoefficient = 0.5; parameterBindiMigliettaFix.shadedSurface = 0.8; + // Wang Leuning parameters parameterWangLeuningFix.stomatalConductanceMin = 0.008; parameterWangLeuningFix.optimalTemperatureForPhotosynthesis = 298.15; + // fenovitis parameters parameterPhenoVitisFix.a= 0.005; parameterPhenoVitisFix.optimalChillingTemp = 2.8; @@ -525,6 +529,8 @@ void Vine3D_Grapevine::radiationAbsorption() { sunlit.leafAreaIndex = 0.0 ; sunlit.absorbedPAR = 0.0 ; + + // TODO: non servono? sunlitAbsorbedNIR = 0.0 ; sunlitAbsorbedLW = 0.0 ; sunlit.isothermalNetRadiation = 0.0 ; @@ -537,12 +543,13 @@ void Vine3D_Grapevine::radiationAbsorption() shadedAbsorbedLW= dum[16] * (UPSCALINGFUNC(diffuseLightK,statePlant.stateGrowth.leafAreaIndex) - UPSCALINGFUNC(directLightK+diffuseLightK,statePlant.stateGrowth.leafAreaIndex)) ; shaded.isothermalNetRadiation = shaded.absorbedPAR + shadedAbsorbedNIR + shadedAbsorbedLW ; } + // Convert absorbed PAR into units of mol m-2 s-1 sunlit.absorbedPAR *= 4.57E-6 ; shaded.absorbedPAR *= 4.57E-6 ; - } + void Vine3D_Grapevine::leafTemperature() { // taken from Hydrall Model, Magnani UNIBO @@ -695,7 +702,8 @@ void Vine3D_Grapevine::aerodynamicalCoupling() { sunlitDeltaTemp = 0.0; } - sunlitDeltaTemp = 0.0; + sunlitDeltaTemp = 0.0; // TODO: check + sunlit.leafTemperature = myInstantTemp + sunlitDeltaTemp + ZEROCELSIUS ; //sunlit big-leaf stomatalConductanceWater = 10.0/shaded.leafAreaIndex ; //dummy stom res for shaded big-leaf //if (shaded.isothermalNetRadiation > 100) stomatalConductanceWater *= pow(100/shaded.isothermalNetRadiation,0.5); @@ -802,10 +810,9 @@ void Vine3D_Grapevine::upscale(TVineCultivar *cultivar) sunlit.darkRespiration = 0.0; sunlit.maximalCarboxylationRate = 0.0; } - - } + void Vine3D_Grapevine::photosynthesisKernel(TVineCultivar* cultivar, double COMP,double GAC,double GHR,double GSCD,double J,double KC,double KO ,double RD,double RNI,double STOMWL,double VCmax,double *ASS,double *GSC,double *TR) { diff --git a/agrolib/grapevine/grapevine.h b/agrolib/grapevine/grapevine.h index a0d18faa..1ff7442c 100644 --- a/agrolib/grapevine/grapevine.h +++ b/agrolib/grapevine/grapevine.h @@ -41,9 +41,9 @@ enum phenoStage {endoDormancy, ecoDormancy , budBurst , flowering , fruitSet, veraison, physiologicalMaturity, vineSenescence}; enum TfieldOperation {irrigationOperation, grassSowing, grassRemoving, trimming, leafRemoval, clusterThinning, harvesting, tartaricAnalysis}; - enum Crit3DLanduse {landuse_nodata, landuse_bare, landuse_vineyard}; + enum GrapevineLanduse {landuse_nodata, landuse_bare, landuse_vineyard}; - const std::map landuseNames = { + const std::map landuseNames = { { "UNDEFINED", landuse_nodata }, { "BARESOIL", landuse_bare }, { "VINEYARD", landuse_vineyard} @@ -191,7 +191,7 @@ struct Crit3DModelCase { int id; - Crit3DLanduse landuse; + GrapevineLanduse landuse; int soilIndex; float shootsPerPlant; @@ -300,8 +300,7 @@ double psiSoilAverage; double psiFieldCapacityAverage; - //double* layerRootDensity; - double totalStomatalConductance, totalStomatalConductanceNoStress ; + double totalStomatalConductance, totalStomatalConductanceNoStress; double transpirationInstant; double* currentProfile; double* transpirationInstantLayer; //molH2O m^-2 s^-1 @@ -327,8 +326,6 @@ bool isAmphystomatic ; double specificLeafArea ; double alphaLeuning ; - //double leafNitrogen ; - //double entropicFactorCarboxyliation,entropicFactorElectronTransporRate ; Vine3D_SunShade shaded ; Vine3D_SunShade sunlit ; @@ -374,7 +371,7 @@ double getWaterStressByPsiSoil(double myPsiSoil,double psiSoilStressParameter,double exponentialFactorForPsiRatio); double getWaterStressSawFunction(int index, TVineCultivar *cultivar); - //bool getExtractedWaterFromGrassTranspirationandEvaporation(double* myWaterExtractionProfile); + double getWaterStressSawFunctionAverage(TVineCultivar* cultivar); double getGrassTranspiration(double stress, double laiGrassMax, double sensitivityToVPD, double fieldCoverByPlant); double getFallowTranspiration(double stress, double laiGrassMax, double sensitivityToVPD); @@ -392,6 +389,11 @@ public: Vine3D_Grapevine(); + TstatePlant getStatePlant(){ return statePlant; } + ToutputPlant getOutputPlant() { return statePlant.outputPlant; } + + bool setStatePlant(TstatePlant myStatePlant, bool isVineyard_); + //void initializeGrapevineModel(TVineCultivar* cultivar, double secondsPerStep); void initializeLayers(int myMaxLayers); bool initializeStatePlant(int doy, Crit3DModelCase *vineField); @@ -409,15 +411,12 @@ double prec , double relativeHumidity , double windSpeed, double atmosphericPressure); bool setDerivedVariables (double diffuseIrradiance, double directIrradiance, double cloudIndex, double sunElevation); - bool setSoilProfile(Crit3DModelCase *modelCase, double* myWiltingPoint, double *myFieldCapacity, - double *myPsiSoilProfile , double *mySoilWaterContentProfile, - double* mySoilWaterContentFC, double* mySoilWaterContentWP); - bool setStatePlant(TstatePlant myStatePlant, bool isVineyard); - TstatePlant getStatePlant(); - ToutputPlant getOutputPlant(); + bool setSoilProfile(Crit3DModelCase* modelCase, std::vector& myWiltingPoint, std::vector& myFieldCapacity, + std::vector& myPsiSoilProfile, std::vector& mySoilWaterContentProfile, + std::vector& mySoilWaterContentFC, std::vector& mySoilWaterContentWP); + double* getExtractedWater(Crit3DModelCase* modelCase); - //bool getOutputPlant(int hour, ToutputPlant *outputPlant); double getStressCoefficient(); double getRealTranspirationGrapevine(Crit3DModelCase *modelCase); double getRealTranspirationGrass(Crit3DModelCase *modelCase); diff --git a/agrolib/graphics/mapGraphicsRasterUtm.cpp b/agrolib/graphics/mapGraphicsRasterUtm.cpp index a45a602f..a8937d86 100644 --- a/agrolib/graphics/mapGraphicsRasterUtm.cpp +++ b/agrolib/graphics/mapGraphicsRasterUtm.cpp @@ -370,7 +370,7 @@ bool RasterUtmObject::drawRaster(QPainter* painter) // check outliers (transparent) if (_rasterPointer->colorScale->isHideOutliers()) { - if (value <= _rasterPointer->colorScale->minimum() || value >= _rasterPointer->colorScale->maximum()) + if (value <= _rasterPointer->colorScale->minimum() || value > _rasterPointer->colorScale->maximum()) continue; } diff --git a/agrolib/homogeneityWidget/homogeneityWidget.cpp b/agrolib/homogeneityWidget/homogeneityWidget.cpp index 7fdc71b6..690b20b6 100644 --- a/agrolib/homogeneityWidget/homogeneityWidget.cpp +++ b/agrolib/homogeneityWidget/homogeneityWidget.cpp @@ -565,7 +565,7 @@ void Crit3DHomogeneityWidget::updateYears() void Crit3DHomogeneityWidget::on_actionChangeLeftAxis() { - DialogChangeAxis changeAxisDialog(true); + DialogChangeAxis changeAxisDialog(1, false); if (changeAxisDialog.result() == QDialog::Accepted) { homogeneityChartView->setYmax(changeAxisDialog.getMaxVal()); diff --git a/agrolib/interpolation/interpolation.cpp b/agrolib/interpolation/interpolation.cpp index 0b91d694..69691eb6 100644 --- a/agrolib/interpolation/interpolation.cpp +++ b/agrolib/interpolation/interpolation.cpp @@ -1015,6 +1015,67 @@ float gaussWeighted(vector &myPointList) } */ +void localSelection_new(std::vector &inputPoints, std::vector &selectedPoints, + float x, float y, float z, Crit3DInterpolationSettings& mySettings) +{ + std::vector tempPoints; + unsigned int i; + float radius; + unsigned int nrValid = 0; + unsigned int minPoints = unsigned(mySettings.getMinPointsLocalDetrending() * 1.2); + float shepardInitialRadius = computeShepardInitialRadius(mySettings.getPointsBoundingBoxArea(), unsigned(inputPoints.size()), minPoints); + + // define a first neighborhood inside initial radius + for (i=0; i < inputPoints.size(); i++) + { + inputPoints[i].distance = gis::computeDistance(x, y, float((inputPoints[i]).point->utm.x), float((inputPoints[i]).point->utm.y)); + if (inputPoints[i].distance <= shepardInitialRadius && + inputPoints[i].distance > 0 && + checkLapseRateCode(inputPoints[i].lapseRateCode, mySettings.getUseLapseRateCode(), true)) + { + tempPoints.push_back(inputPoints[i]); + nrValid++; + } + } + + if (tempPoints.size() <= minPoints) + { + nrValid = sortPointsByDistance(minPoints + 1, inputPoints, selectedPoints); + if (nrValid > minPoints) + { + radius = selectedPoints[minPoints].distance; + selectedPoints.pop_back(); + } + else + radius = selectedPoints[nrValid-1].distance + float(EPSILON); + } + else if (tempPoints.size() > minPoints) + { + nrValid = sortPointsByDistance(minPoints + 1, tempPoints, selectedPoints); + radius = selectedPoints[minPoints].distance; + selectedPoints.pop_back(); + } + else + { + selectedPoints = tempPoints; + radius = shepardInitialRadius; + } + + for (int i = 0; i < selectedPoints.size(); i++) + { + selectedPoints[i].regressionWeight = MAXVALUE((-(1/std::pow(radius,4)*(std::pow(selectedPoints[i].distance,4)))+1),EPSILON); + //selectedPoints[i].regressionWeight = 1; + //selectedPoints[i].heightWeight = 1./((2./maxHeightDelta)*selectedPoints[i].point->z+1); + selectedPoints[i].heightWeight = 1; + } + mySettings.setLocalRadius(float(radius)); + + return; +} + + + + // TODO elevation std dev? void localSelection(vector &inputPoints, vector &selectedPoints, float x, float y, float z, Crit3DInterpolationSettings& mySettings) @@ -1033,7 +1094,7 @@ void localSelection(vector &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 = 5000; // [m] + float stepRadius = 2500; // [m] float r0 = 0; // [m] float r1 = stepRadius; // [m] unsigned int i; @@ -1077,7 +1138,6 @@ void localSelection(vector &inputPoints, vector < mySettings.setLocalRadius(float(maxDistance)); } - bool checkPrecipitationZero(const std::vector &myPoints, float precThreshold, int &nrNotNull) { nrNotNull = 0; @@ -1100,6 +1160,7 @@ bool isThermal(meteoVariable myVar) myVar == dailyAirTemperatureAvg || myVar == dailyAirTemperatureMax || myVar == dailyAirTemperatureMin || + myVar == dailyReferenceEvapotranspirationHS || myVar == elaboration ) return true; else @@ -1114,6 +1175,7 @@ bool getUseDetrendingVar(meteoVariable myVar) myVar == dailyAirTemperatureAvg || myVar == dailyAirTemperatureMax || myVar == dailyAirTemperatureMin || + myVar == dailyReferenceEvapotranspirationHS || myVar == elaboration ) return true; @@ -1284,6 +1346,7 @@ bool regressionOrography(std::vector &myPoints, } } + void detrending(std::vector &myPoints, Crit3DProxyCombination myCombination, Crit3DInterpolationSettings* mySettings, Crit3DClimateParameters* myClimate, meteoVariable myVar, Crit3DTime myTime) @@ -1405,19 +1468,33 @@ bool proxyValidityWeighted(std::vector &myPoints, return true; } -bool setAllFittingRanges(Crit3DProxyCombination myCombination, Crit3DInterpolationSettings* mySettings) +bool setHeightFittingRange(Crit3DProxyCombination myCombination, Crit3DInterpolationSettings* mySettings) { if (mySettings->getMinMaxTemperature().empty()) return 0; + const double H0_MIN = -200; //height of inversion point (double piecewise) or first inversion point (triple piecewise) + const double H0_MAX = 5000; + const double DELTA_MIN = 300; //height difference between inversion points (for triple piecewise only) + const double DELTA_MAX = 1000; + const double SLOPE_MIN = 0.002; //ascending slope + const double SLOPE_MAX = 0.007; + const double INVSLOPE_MIN = -0.01; //inversion slope + const double INVSLOPE_MAX = -0.0015; + for (unsigned i=0; i < myCombination.getProxySize(); i++) if (myCombination.isProxyActive(i) == true) { if (getProxyPragaName(mySettings->getProxy(i)->getName()) == proxyHeight) { - double min = mySettings->getMinMaxTemperature()[0]; - double max = mySettings->getMinMaxTemperature()[1]; - + const double MIN_T = mySettings->getMinMaxTemperature()[0]; + const double MAX_T = mySettings->getMinMaxTemperature()[1]; + + /* + * following line allows to check if the function for elevation has been changed (GUI only) compared to the + * function read in the .ini file. if it hasn't been changed, only the minimum and maximum temperature get rewritten. + * otherwise appropriate parameters are loaded into the proxy (fittingParametersRange) + */ if (mySettings->getChosenElevationFunction() == mySettings->getProxy(i)->getFittingFunctionName()) { std::vector tempParam; @@ -1426,18 +1503,18 @@ bool setAllFittingRanges(Crit3DProxyCombination myCombination, Crit3DInterpolati { if (mySettings->getChosenElevationFunction() == piecewiseTwo) { - tempParam[1] = min-2; - tempParam[5] = max+2; + tempParam[1] = MIN_T-2; + tempParam[5] = MAX_T+2; } else if (mySettings->getChosenElevationFunction() == piecewiseThreeFree) { - tempParam[1] = min-2; - tempParam[7] = max+2; + tempParam[1] = MIN_T-2; + tempParam[7] = MAX_T+2; } else if (mySettings->getChosenElevationFunction() == piecewiseThree) { - tempParam[1] = min-2; - tempParam[6] = max+2; + tempParam[1] = MIN_T-2; + tempParam[6] = MAX_T+2; } mySettings->getProxy(i)->setFittingParametersRange(tempParam); } @@ -1448,22 +1525,20 @@ bool setAllFittingRanges(Crit3DProxyCombination myCombination, Crit3DInterpolati if (mySettings->getChosenElevationFunction() == piecewiseTwo) { mySettings->getProxy(i)->setFittingFunctionName(piecewiseTwo); - tempParam = {-200, min-2, 0.002, -0.01, 5000, max+2, 0.01, 0.0015}; + tempParam = {H0_MIN, MIN_T-2, SLOPE_MIN, INVSLOPE_MIN, + H0_MAX, MAX_T+2, SLOPE_MAX, INVSLOPE_MAX}; } else if (mySettings->getChosenElevationFunction() == piecewiseThreeFree) { mySettings->getProxy(i)->setFittingFunctionName(piecewiseThreeFree); - tempParam = {-200, min-2, 300, 0.002, -0.01, -0.01, 5000, max+2, 1000, 0.007, 0.0015, 0.0015}; + tempParam = {H0_MIN, MIN_T-2, DELTA_MIN, SLOPE_MIN, INVSLOPE_MIN, INVSLOPE_MIN, + H0_MAX, MAX_T+2, DELTA_MAX, SLOPE_MAX, INVSLOPE_MAX, INVSLOPE_MAX}; } else if (mySettings->getChosenElevationFunction() == piecewiseThree) { mySettings->getProxy(i)->setFittingFunctionName(piecewiseThree); - tempParam = {-200, min-2, 300, 0.002, -0.01, 5000, max+2, 1000, 0.007, 0.0015}; - } - else if (mySettings->getChosenElevationFunction() == freiFree) - { - mySettings->getProxy(i)->setFittingFunctionName(freiFree); - tempParam = {min, 0, -4, -200, 0.1, 0.002, max+10, 0.006, 4, 5000, 1000, 0.006}; + tempParam = {H0_MIN, MIN_T-2, DELTA_MIN, SLOPE_MIN, INVSLOPE_MIN, + H0_MAX, MAX_T+2, DELTA_MAX, SLOPE_MAX, INVSLOPE_MAX}; } mySettings->getProxy(i)->setFittingParametersRange(tempParam); } @@ -1504,11 +1579,9 @@ bool setAllFittingParameters_noRange(Crit3DProxyCombination myCombination, Crit3 if (mySettings->getChosenElevationFunction() == piecewiseTwo) myFunc[i] = lapseRatePiecewise_two; else if (mySettings->getChosenElevationFunction() == piecewiseThreeFree) - myFunc[i] = lapseRatePiecewiseFree; + myFunc[i] = lapseRatePiecewise_three_free; else if (mySettings->getChosenElevationFunction() == piecewiseThree) - myFunc[i] = lapseRatePiecewiseThree_withSlope; - else if (mySettings->getChosenElevationFunction() == freiFree) - myFunc[i] = lapseRateFreiFree; + myFunc[i] = lapseRatePiecewise_three; else { errorStr = "Missing or wrong fitting function for proxy: " + mySettings->getProxy(i)->getName(); @@ -1581,25 +1654,18 @@ bool setAllFittingParameters(Crit3DProxyCombination myCombination, Crit3DInterpo } else if (mySettings->getChosenElevationFunction() == piecewiseThreeFree) { - myFunc.push_back(lapseRatePiecewiseFree); + myFunc.push_back(lapseRatePiecewise_three_free); mySettings->getProxy(i)->setFittingFunctionName(piecewiseThreeFree); if (!(mySettings->getProxy(i)->getFittingParametersRange().empty())) tempParam = {-200, min-2, 100, 0.001, -0.006, -0.006, 1800, max+2, 1000, 0.01, 0, 0}; } else if (mySettings->getChosenElevationFunction() == piecewiseThree) { - myFunc.push_back(lapseRatePiecewiseThree_withSlope); + myFunc.push_back(lapseRatePiecewise_three); mySettings->getProxy(i)->setFittingFunctionName(piecewiseThree); if (!(mySettings->getProxy(i)->getFittingParametersRange().empty())) tempParam = {-200, min-2, 100, 0.002, -0.006, 1800, max+2, 1000, 0.01, 0}; } - else if (mySettings->getChosenElevationFunction() == freiFree) - { - myFunc.push_back(lapseRateFreiFree); - mySettings->getProxy(i)->setFittingFunctionName(freiFree); - if (!(mySettings->getProxy(i)->getFittingParametersRange().empty())) - tempParam = {min, 0, -4, -200, 0.1, 0, max+10, 0.006, 4, 1800, 1000, 0.006}; - } mySettings->getProxy(i)->setFittingParametersRange(tempParam); } else @@ -1882,12 +1948,12 @@ bool multipleDetrendingElevation(Crit3DProxyCombination elevationCombination, st mySettings->setSingleFittingFunction(detrendingLapseRatePiecewise_two, elevationPos); } else if (mySettings->getProxy(elevationPos)->getFittingFunctionName() == piecewiseThreeFree) { - myFunc[elevationPos] = detrendingLapseRatePiecewiseFree; - mySettings->setSingleFittingFunction(detrendingLapseRatePiecewiseFree, elevationPos); + myFunc[elevationPos] = detrendingLapseRatePiecewise_three_free; + mySettings->setSingleFittingFunction(detrendingLapseRatePiecewise_three_free, elevationPos); } else if (mySettings->getProxy(elevationPos)->getFittingFunctionName() == piecewiseThree) { - myFunc[elevationPos] = detrendingLapseRatePiecewiseThree_withSlope; - mySettings->setSingleFittingFunction(detrendingLapseRatePiecewiseThree_withSlope, elevationPos); + myFunc[elevationPos] = detrendingLapseRatePiecewise_three; + mySettings->setSingleFittingFunction(detrendingLapseRatePiecewise_three, elevationPos); } func = myFunc[elevationPos].target&)>();; @@ -2193,7 +2259,7 @@ void topographicDistanceOptimize(meteoVariable myVar, mySettings->setTopoDist_Kh(kh); if (computeResiduals(myVar, myMeteoPoints, nrMeteoPoints, interpolationPoints, mySettings, meteoSettings, true, true)) { - avgError = computeErrorCrossValidation(myVar, myMeteoPoints, nrMeteoPoints, myTime, meteoSettings); + avgError = computeErrorCrossValidation(myMeteoPoints, nrMeteoPoints); if (isEqual(bestError, NODATA) || avgError < bestError) { bestError = avgError; @@ -2240,7 +2306,7 @@ void optimalDetrending(meteoVariable myVar, Crit3DMeteoPoint* &myMeteoPoints, in if (computeResiduals(myVar, myMeteoPoints, nrMeteoPoints, interpolationPoints, mySettings, meteoSettings, true, true)) { - avgError = computeErrorCrossValidation(myVar, myMeteoPoints, nrMeteoPoints, myTime, meteoSettings); + avgError = computeErrorCrossValidation(myMeteoPoints, nrMeteoPoints); if (! isEqual(avgError, NODATA) && (isEqual(minError, NODATA) || avgError < minError)) { minError = avgError; @@ -2304,7 +2370,9 @@ bool preInterpolation(std::vector &myPoints, Crit } if (mySettings->getUseTD() && getUseTdVar(myVar)) + { topographicDistanceOptimize(myVar, myMeteoPoints, nrMeteoPoints, myPoints, mySettings, meteoSettings, myTime); + } return true; } diff --git a/agrolib/interpolation/interpolation.h b/agrolib/interpolation/interpolation.h index 013a12f0..88cbff1e 100644 --- a/agrolib/interpolation/interpolation.h +++ b/agrolib/interpolation/interpolation.h @@ -100,13 +100,16 @@ std::vector &selectedPoints, float x, float y, float z, Crit3DInterpolationSettings &mySettings); + void localSelection_new(std::vector &inputPoints, std::vector &selectedPoints, + float x, float y, float z, Crit3DInterpolationSettings& mySettings); + bool proxyValidity(std::vector &myPoints, int proxyPos, float stdDevThreshold, double &avg, double &stdDev); bool proxyValidityWeighted(std::vector &myPoints, int proxyPos, float stdDevThreshold, double &avg, double &stdDev); - bool setAllFittingRanges(Crit3DProxyCombination myCombination, Crit3DInterpolationSettings* mySettings); + bool setHeightFittingRange(Crit3DProxyCombination myCombination, Crit3DInterpolationSettings* mySettings); bool setAllFittingParameters_noRange(Crit3DProxyCombination myCombination, Crit3DInterpolationSettings* mySettings, std::vector&)>>& myFunc, diff --git a/agrolib/interpolation/interpolationConstants.h b/agrolib/interpolation/interpolationConstants.h index 5cb13de0..c6b3013a 100644 --- a/agrolib/interpolation/interpolationConstants.h +++ b/agrolib/interpolation/interpolationConstants.h @@ -38,13 +38,12 @@ { "water_index", proxyWaterIndex} }; - enum TFittingFunction { piecewiseTwo, piecewiseThreeFree, piecewiseThree, freiFree, linear, noFunction }; + enum TFittingFunction { piecewiseTwo, piecewiseThreeFree, piecewiseThree, linear, noFunction }; const std::map fittingFunctionNames = { { "piecewise_two", piecewiseTwo }, { "free_triple_piecewise", piecewiseThreeFree}, { "triple_piecewise", piecewiseThree}, - { "Nonlinear Frei function (6 parameters)", freiFree }, { "linear", linear } }; diff --git a/agrolib/interpolation/spatialControl.cpp b/agrolib/interpolation/spatialControl.cpp index f7834caf..932e0130 100644 --- a/agrolib/interpolation/spatialControl.cpp +++ b/agrolib/interpolation/spatialControl.cpp @@ -84,9 +84,6 @@ bool computeResiduals(meteoVariable myVar, Crit3DMeteoPoint* meteoPoints, int nr if (myVar == noMeteoVar) return false; - float myValue, interpolatedValue; - interpolatedValue = NODATA; - myValue = NODATA; std::vector myProxyValues; bool isValid; @@ -101,58 +98,67 @@ bool computeResiduals(meteoVariable myVar, Crit3DMeteoPoint* meteoPoints, int nr if (isValid && meteoPoints[i].quality == quality::accepted) { - myValue = meteoPoints[i].currentValue; + float myValue = meteoPoints[i].currentValue; - interpolatedValue = interpolate(interpolationPoints, settings, meteoSettings, myVar, + float interpolatedValue = interpolate(interpolationPoints, settings, meteoSettings, myVar, float(meteoPoints[i].point.utm.x), float(meteoPoints[i].point.utm.y), float(meteoPoints[i].point.z), myProxyValues, false); - if ( myVar == precipitation - || myVar == dailyPrecipitation) + if ( myVar == precipitation || myVar == dailyPrecipitation) { if (myValue != NODATA) - if (myValue < meteoSettings->getRainfallThreshold()) myValue=0.; + { + if (myValue < meteoSettings->getRainfallThreshold()) + myValue=0.; + } if (interpolatedValue != NODATA) - if (interpolatedValue < meteoSettings->getRainfallThreshold()) interpolatedValue=0.; + { + if (interpolatedValue < meteoSettings->getRainfallThreshold()) + interpolatedValue=0.; + } } // TODO derived var if ((interpolatedValue != NODATA) && (myValue != NODATA)) + { meteoPoints[i].residual = interpolatedValue - myValue; + } } } return true; } -float computeErrorCrossValidation(meteoVariable myVar, Crit3DMeteoPoint* myPoints, int nrMeteoPoints, const Crit3DTime& myTime, Crit3DMeteoSettings* meteoSettings) + +float computeErrorCrossValidation(Crit3DMeteoPoint* myPoints, int nrMeteoPoints) { std::vector obsValues, estValues; - float myValue, myEstimate, myResidual; for (int i=0; i < nrMeteoPoints; i++) { if (myPoints[i].active) { - myValue = myPoints[i].getMeteoPointValue(myTime, myVar, meteoSettings); - myResidual = myPoints[i].residual; + float value = myPoints[i].currentValue; + float residual = myPoints[i].residual; - if (myValue != NODATA && myResidual != NODATA) + if (value != NODATA && residual != NODATA) { - myEstimate = myValue + myResidual; - obsValues.push_back(myValue); - estValues.push_back(myEstimate); + obsValues.push_back(value); + estValues.push_back(value + residual); } } } if (obsValues.size() > 0) + { return statistics::meanAbsoluteError(obsValues, estValues); - else return NODATA; + } + else + return NODATA; } diff --git a/agrolib/interpolation/spatialControl.h b/agrolib/interpolation/spatialControl.h index 1cb2d7e9..347efdb8 100644 --- a/agrolib/interpolation/spatialControl.h +++ b/agrolib/interpolation/spatialControl.h @@ -31,7 +31,7 @@ bool computeResiduals(meteoVariable myVar, Crit3DMeteoPoint* meteoPoints, int nrMeteoPoints, std::vector &interpolationPoints, Crit3DInterpolationSettings* settings, Crit3DMeteoSettings* meteoSettings, bool excludeOutsideDem, bool excludeSupplemental); - float computeErrorCrossValidation(meteoVariable myVar, Crit3DMeteoPoint *myPoints, int nrMeteoPoints, const Crit3DTime& myTime, Crit3DMeteoSettings *meteoSettings); + float computeErrorCrossValidation(Crit3DMeteoPoint *myPoints, int nrMeteoPoints); bool spatialQualityControl(meteoVariable myVar, Crit3DMeteoPoint* meteoPoints, int nrMeteoPoints, Crit3DInterpolationSettings *settings, Crit3DMeteoSettings* meteoSettings, diff --git a/agrolib/mathFunctions/furtherMathFunctions.cpp b/agrolib/mathFunctions/furtherMathFunctions.cpp index d85a3e29..a85d650f 100644 --- a/agrolib/mathFunctions/furtherMathFunctions.cpp +++ b/agrolib/mathFunctions/furtherMathFunctions.cpp @@ -69,33 +69,8 @@ double lapseRateFrei(double x, std::vector & par) return y - 0.5*par[2]*(1 + cos(PI*(x-par[3])/(par[4]))); } -double lapseRateFreiFree(double x, std::vector & par) -{ - /* - par[0] = T0; - par[1] = gamma1; - par[2] = a; - par[3] = h0; - par[4] = h1-h0; - par[5] = gamma2 - */ - - if (par.size() < 6) return NODATA; - double h1 = par[3]+par[4]; - if (x <= par[3]) - { - return par[0] - par[1]*x - par[2]; - } - else if (x >= (par[4]+par[3])) - { - return par[0] - par[5]*x; - } - return par[0] - ((par[5]*par[3]+par[1]*h1)/(par[3]+h1))*x - 0.5*par[2]*(1 + cos(PI*(x-par[3])/(par[4]))); -} - - -double lapseRatePiecewise_three(double x, std::vector & par) +double lapseRatePiecewise_three_noSlope(double x, std::vector & par) { // the piecewise line is parameterized as follows // the line passes through A(par[0];par[1])and B(par[0]+par[2];par[3]). par[4] is the slope of the 2 externals pieces @@ -124,7 +99,7 @@ double lapseRatePiecewise_three(double x, std::vector & par) } -double lapseRatePiecewiseThree_withSlope(double x, std::vector & par) +double lapseRatePiecewise_three(double x, std::vector & par) { //xa (par[0],par[1]), xb-xa = par[2], par[3] is the slope of the middle piece, //par[4] the slope of the first and last piece @@ -139,7 +114,7 @@ double lapseRatePiecewiseThree_withSlope(double x, std::vector & par) return par[3]*x - par[3]*par[0]+par[1]; } -double detrendingLapseRatePiecewiseThree_withSlope(double x, std::vector & par) +double detrendingLapseRatePiecewise_three(double x, std::vector & par) { //xa (par[0],par[1]), xb-xa = par[2], par[3] is the slope of the middle piece, //par[4] the slope of the first and last piece @@ -154,7 +129,7 @@ double detrendingLapseRatePiecewiseThree_withSlope(double x, std::vector & par) +double lapseRatePiecewise_three_free(double x, std::vector & par) { // the piecewise line is parameterized as follows // the line passes through A(par[0];par[1])and B(par[0]+par[2];...) @@ -186,7 +161,7 @@ double lapseRatePiecewiseFree(double x, std::vector & par) } } -double detrendingLapseRatePiecewiseFree(double x, std::vector & par) +double detrendingLapseRatePiecewise_three_free(double x, std::vector & par) { // the piecewise line is parameterized as follows // the line passes through A(par[0];par[1])and B(par[0]+par[2];...) @@ -1346,7 +1321,7 @@ namespace interpolation //grigliato const int numSteps = 20; - std::vector steps = {2*(parametersMax[0][0]-parametersMin[0][0])/numSteps, 2*(parametersMax[0][1]-parametersMin[0][1])/numSteps}; + std::vector stepSize = {2*(parametersMax[0][0]-parametersMin[0][0])/numSteps, 2*(parametersMax[0][1]-parametersMin[0][1])/numSteps}; int directions[] = {1, -1}; std::vector> firstGuessParam = parameters; @@ -1405,9 +1380,9 @@ namespace interpolation if (parameters.size() > proxyIndex && !parameters[proxyIndex].empty()) { if (dir == 0) - parameters[proxyIndex][paramIndex] = MINVALUE(firstGuessParam[proxyIndex][paramIndex] + directions[dir] * step * steps[paramIndex], parametersMax[proxyIndex][paramIndex]); + parameters[proxyIndex][paramIndex] = MINVALUE(firstGuessParam[proxyIndex][paramIndex] + directions[dir] * step * stepSize[paramIndex], parametersMax[proxyIndex][paramIndex]); else - parameters[proxyIndex][paramIndex] = MAXVALUE(firstGuessParam[proxyIndex][paramIndex] + directions[dir] * step * steps[paramIndex], parametersMin[proxyIndex][paramIndex]); + parameters[proxyIndex][paramIndex] = MAXVALUE(firstGuessParam[proxyIndex][paramIndex] + directions[dir] * step * stepSize[paramIndex], parametersMin[proxyIndex][paramIndex]); } if ((counter > nrTrials) || ((R2Previous[0] != NODATA) && fabs(R2Previous[0]-R2Previous[nrMinima-1]) < deltaR2 )) @@ -2068,12 +2043,12 @@ namespace interpolation int counter = 0; //grigliato - std::vector steps; + std::vector stepSize; const int numSteps = 30; if (parameters.size() == 4) - steps = {2*(parametersMax[0]-parametersMin[0])/numSteps, 2*(parametersMax[1]-parametersMin[1])/numSteps, 2*(parametersMax[2]-parametersMin[2])/numSteps, 2*(parametersMax[3]-parametersMin[3])/numSteps}; + stepSize = {2*(parametersMax[0]-parametersMin[0])/numSteps, 2*(parametersMax[1]-parametersMin[1])/numSteps, 20*(parametersMax[2]-parametersMin[2])/numSteps, 20*(parametersMax[3]-parametersMin[3])/numSteps}; else if (parameters.size() == 6) - steps = {2*(parametersMax[0]-parametersMin[0])/numSteps, 2*(parametersMax[1]-parametersMin[1])/numSteps, 4*(parametersMax[2]-parametersMin[2])/numSteps, 20*(parametersMax[3]-parametersMin[3])/numSteps,20*(parametersMax[3]-parametersMin[3])/numSteps,20*(parametersMax[3]-parametersMin[3])/numSteps }; + stepSize = {2*(parametersMax[0]-parametersMin[0])/numSteps, 2*(parametersMax[1]-parametersMin[1])/numSteps, 4*(parametersMax[2]-parametersMin[2])/numSteps, 20*(parametersMax[3]-parametersMin[3])/numSteps,20*(parametersMax[3]-parametersMin[3])/numSteps,20*(parametersMax[3]-parametersMin[3])/numSteps }; else return false; int directions[] = {1, -1}; @@ -2122,9 +2097,9 @@ namespace interpolation counter++; if (dir == 0) - parameters[paramIndex] = MINVALUE(firstGuessParam[paramIndex] + directions[dir] * step * steps[paramIndex], parametersMax[paramIndex]); + parameters[paramIndex] = MINVALUE(firstGuessParam[paramIndex] + directions[dir] * step * stepSize[paramIndex], parametersMax[paramIndex]); else - parameters[paramIndex] = MAXVALUE(firstGuessParam[paramIndex] + directions[dir] * step * steps[paramIndex], parametersMin[paramIndex]); + parameters[paramIndex] = MAXVALUE(firstGuessParam[paramIndex] + directions[dir] * step * stepSize[paramIndex], parametersMin[paramIndex]); } } diff --git a/agrolib/mathFunctions/furtherMathFunctions.h b/agrolib/mathFunctions/furtherMathFunctions.h index c5a3fc5b..0b9aa663 100644 --- a/agrolib/mathFunctions/furtherMathFunctions.h +++ b/agrolib/mathFunctions/furtherMathFunctions.h @@ -50,18 +50,20 @@ enum estimatedFunction {FUNCTION_CODE_SPHERICAL, FUNCTION_CODE_LINEAR, FUNCTION_ double functionLinear(double x, std::vector & par); double functionLinear_intercept(double x, std::vector & par); double multilinear(std::vector &x, std::vector &par); - double lapseRatePiecewise_three(double x, std::vector & par); - double lapseRatePiecewiseForInterpolation(double x, std::vector & par); - double lapseRatePiecewiseFree(double x, std::vector & par); - double detrendingLapseRatePiecewiseFree(double x, std::vector & par); - double lapseRatePiecewiseThree_withSlope(double x, std::vector & par); - double detrendingLapseRatePiecewiseThree_withSlope(double x, std::vector & par); - double lapseRatePiecewise_two(double x, std::vector & par); - double detrendingLapseRatePiecewise_two(double x, std::vector & par); double lapseRateFrei(double x, std::vector & par); double lapseRateFreiFree(double x, std::vector & par); double lapseRateRotatedSigmoid(double x, std::vector par); + double lapseRatePiecewise_two(double x, std::vector & par); + double lapseRatePiecewise_three_noSlope(double x, std::vector & par); + double lapseRatePiecewise_three(double x, std::vector & par); + double lapseRatePiecewise_three_free(double x, std::vector & par); + + double detrendingLapseRatePiecewise_two(double x, std::vector & par); + double detrendingLapseRatePiecewise_three(double x, std::vector & par); + double detrendingLapseRatePiecewise_three_free(double x, std::vector & par); + + namespace integration { float trapzdParametric(float (*func)(TfunctionInput), int nrPar, float *par , float a , float b , int n); diff --git a/agrolib/meteo/meteo.cpp b/agrolib/meteo/meteo.cpp index 21709b86..df096947 100644 --- a/agrolib/meteo/meteo.cpp +++ b/agrolib/meteo/meteo.cpp @@ -927,6 +927,10 @@ std::string getVariableString(meteoVariable myVar) else if (myVar == leafAreaIndex) return "Leaf area index (m2 m-2)"; + else if (myVar == elaboration) + return "Elaboration"; + else if (myVar == anomaly) + return "Anomaly"; else if (myVar == noMeteoTerrain) return "Elevation (m)"; else diff --git a/agrolib/meteoWidget/meteoWidget.cpp b/agrolib/meteoWidget/meteoWidget.cpp index 20e0e43a..5324411c 100644 --- a/agrolib/meteoWidget/meteoWidget.cpp +++ b/agrolib/meteoWidget/meteoWidget.cpp @@ -64,7 +64,7 @@ Crit3DMeteoWidget::Crit3DMeteoWidget(bool isGrid_, QString projectPath, Crit3DMe this->resize(1240, 700); this->setAttribute(Qt::WA_DeleteOnClose); - currentFreq = noFrequency; + _currentFrequency = noFrequency; QDate noDate = QDate(1800,1,1); currentDate = noDate; @@ -140,11 +140,11 @@ Crit3DMeteoWidget::Crit3DMeteoWidget(bool isGrid_, QString projectPath, Crit3DMe } if (key.contains("DAILY")) { - currentFreq = daily; + _currentFrequency = daily; } else { - currentFreq = hourly; + _currentFrequency = hourly; } MapCSVDefault.insert(key,items); zeroLine = new QLineSeries(); @@ -207,7 +207,7 @@ Crit3DMeteoWidget::Crit3DMeteoWidget(bool isGrid_, QString projectPath, Crit3DMe if (currentVariables.isEmpty() || (dailyVar != 0 && hourlyVar != 0)) { QMessageBox::information(nullptr, "Warning", "Wrong variables in Crit3DPlotDefault.csv"); - currentFreq = noFrequency; + _currentFrequency = noFrequency; currentVariables.clear(); nameLines.clear(); nameBar.clear(); @@ -293,19 +293,19 @@ Crit3DMeteoWidget::Crit3DMeteoWidget(bool isGrid_, QString projectPath, Crit3DMe firstDate->setMinimumWidth(firstDate->width()-firstDate->width()*0.3); lastDate->setMinimumWidth(lastDate->width()-lastDate->width()*0.3); - if (currentFreq == daily || currentFreq == noFrequency) + if (_currentFrequency == daily || _currentFrequency == noFrequency) { dailyButton->setEnabled(false); hourlyButton->setEnabled(true); monthlyButton->setEnabled(true); } - else if (currentFreq == hourly) + else if (_currentFrequency == hourly) { hourlyButton->setEnabled(false); dailyButton->setEnabled(true); monthlyButton->setEnabled(true); } - else if (currentFreq == monthly) + else if (_currentFrequency == monthly) { monthlyButton->setEnabled(false); dailyButton->setEnabled(true); @@ -421,24 +421,24 @@ Crit3DMeteoWidget::~Crit3DMeteoWidget() } -void Crit3DMeteoWidget::setCurrentFrequency(frequencyType frequency) +void Crit3DMeteoWidget::setFrequency(frequencyType frequency) { - currentFreq = frequency; + _currentFrequency = frequency; // update gui - if (currentFreq == daily || currentFreq == noFrequency) + if (_currentFrequency == daily || _currentFrequency == noFrequency) { dailyButton->setEnabled(false); hourlyButton->setEnabled(true); monthlyButton->setEnabled(true); } - else if (currentFreq == hourly) + else if (_currentFrequency == hourly) { hourlyButton->setEnabled(false); dailyButton->setEnabled(true); monthlyButton->setEnabled(true); } - else if (currentFreq == monthly) + else if (_currentFrequency == monthly) { monthlyButton->setEnabled(false); dailyButton->setEnabled(true); @@ -582,11 +582,11 @@ void Crit3DMeteoWidget::drawMeteoPoint(Crit3DMeteoPoint mp, bool isAppend) lastDate->setDate(currentDate); // draw period (31 days for daily, 3 days for hourly) - if (currentFreq == daily) + if (_currentFrequency == daily) { firstDate->setDate(currentDate.addDays(-30)); } - else if (currentFreq == hourly) + else if (_currentFrequency == hourly) { firstDate->setDate(currentDate.addDays(-2)); } @@ -680,9 +680,9 @@ void Crit3DMeteoWidget::drawEnsemble() lastDate->setDate(currentDate); // draw period (31 days for daily, 3 days for hourly) - if (currentFreq == daily) + if (_currentFrequency == daily) firstDate->setDate(currentDate.addDays(-30)); - else if (currentFreq == hourly) + else if (_currentFrequency == hourly) firstDate->setDate(currentDate.addDays(-2)); redraw(); @@ -2052,19 +2052,19 @@ void Crit3DMeteoWidget::showVar() { if (! isInitialized) return; - if (currentFreq == noFrequency) + if (_currentFrequency == noFrequency) { if (dailyButton->isChecked()) // dailyButton is pressed { - currentFreq = daily; + _currentFrequency = daily; } else if (hourlyButton->isChecked()) { - currentFreq = hourly; + _currentFrequency = hourly; } else if (monthlyButton->isChecked()) { - currentFreq = monthly; + _currentFrequency = monthly; } } QList allKeys = MapCSVStyles.keys(); @@ -2072,21 +2072,21 @@ void Crit3DMeteoWidget::showVar() QList allVar; for (int i = 0; idate().daysTo(lastDate->date()); QDate firstValidDate; - if (currentFreq == daily) + if (_currentFrequency == daily) { if (! firstDailyDate.isValid() || firstDailyDate.year() == 1800) return; firstValidDate = firstDailyDate; } - else if (currentFreq == hourly) + else if (_currentFrequency == hourly) { if (! firstHourlyDate.isValid() || firstHourlyDate.year() == 1800) return; firstValidDate = firstHourlyDate; } - else if (currentFreq == monthly) + else if (_currentFrequency == monthly) { if (! firstMonthlyDate.isValid() || firstMonthlyDate.year() == 1800) return; @@ -2477,7 +2477,7 @@ void Crit3DMeteoWidget::shiftFollowing() int nDays = firstDate->date().daysTo(lastDate->date()); QDate lastValidDate; - if (currentFreq == daily) + if (_currentFrequency == daily) { if (! lastDailyDate.isValid() || lastDailyDate.year() == 1800) return; @@ -2509,7 +2509,7 @@ void Crit3DMeteoWidget::shiftFollowing() void Crit3DMeteoWidget::showTable() { - DialogMeteoTable meteoTable(meteoSettings, meteoPoints, firstDate->date(), lastDate->date(), currentFreq, currentVariables); + DialogMeteoTable meteoTable(meteoSettings, meteoPoints, firstDate->date(), lastDate->date(), _currentFrequency, currentVariables); } void Crit3DMeteoWidget::tooltipLineSeries(QPointF point, bool state) @@ -2608,18 +2608,18 @@ bool Crit3DMeteoWidget::computeTooltipLineSeries(QLineSeries *series, QPointF po if (!valueExist) { // missing data - if (currentFreq == daily) + if (_currentFrequency == daily) { QDate xDate = firstDate->date().addDays(doy); m_tooltip->setText(QString("%1 \n%2 nan ").arg(series->name()).arg(xDate.toString("MMM dd yyyy"))); } - else if (currentFreq == hourly) + else if (_currentFrequency == hourly) { QDateTime xDate(firstDate->date(), QTime(0,0,0), Qt::UTC); xDate = xDate.addSecs(3600*doy); m_tooltip->setText(QString("%1 \n%2 nan ").arg(series->name()).arg(xDate.toString("MMM dd yyyy hh:mm"))); } - else if (currentFreq == monthly) + else if (_currentFrequency == monthly) { QDate xDate = firstDate->date().addMonths(doy); m_tooltip->setText(QString("%1 \n%2 nan ").arg(series->name()).arg(xDate.toString("MMM yyyy"))); @@ -2739,7 +2739,7 @@ bool Crit3DMeteoWidget::computeTooltipLineSeries(QLineSeries *series, QPointF po } - if (currentFreq == daily) + if (_currentFrequency == daily) { QDate xDate = firstDate->date().addDays(doy); for(int i = 0; i < series->count(); i++) @@ -2753,7 +2753,7 @@ bool Crit3DMeteoWidget::computeTooltipLineSeries(QLineSeries *series, QPointF po double value = series->at(doyRelative).y(); m_tooltip->setText(QString("%1 \n%2 %3 ").arg(series->name()).arg(xDate.toString("MMM dd yyyy")).arg(value, 0, 'f', 1)); } - else if (currentFreq == hourly) + else if (_currentFrequency == hourly) { QDateTime xDate(firstDate->date(), QTime(0,0,0), Qt::UTC); xDate = xDate.addSecs(3600*doy); @@ -2768,7 +2768,7 @@ bool Crit3DMeteoWidget::computeTooltipLineSeries(QLineSeries *series, QPointF po double value = series->at(doyRelative).y(); m_tooltip->setText(QString("%1 \n%2 %3 ").arg(series->name()).arg(xDate.toString("MMM dd yyyy hh:mm")).arg(value, 0, 'f', 1)); } - else if (currentFreq == monthly) + else if (_currentFrequency == monthly) { QDate xDate = firstDate->date().addMonths(doy); for(int i = 0; i < series->count(); i++) @@ -2829,19 +2829,19 @@ void Crit3DMeteoWidget::tooltipBar(bool state, int index, QBarSet *barset) } QString valueStr; - if (currentFreq == daily) + if (_currentFrequency == daily) { QDate xDate = firstDate->date().addDays(index); valueStr = QString("%1 \n%2 %3 ").arg(xDate.toString("MMM dd yyyy")).arg(barset->label()).arg(barset->at(index), 0, 'f', 1); } - else if (currentFreq == hourly) + else if (_currentFrequency == hourly) { QDateTime xDate(firstDate->date(), QTime(0,0,0), Qt::UTC); xDate = xDate.addSecs(3600*index); valueStr = QString("%1 \n%2 %3 ").arg(xDate.toString("MMM dd yyyy hh:mm")).arg(barset->label()).arg(barset->at(index), 0, 'f', 1); } - else if (currentFreq == monthly) + else if (_currentFrequency == monthly) { QDate xDate = firstDate->date().addMonths(index); valueStr = QString("%1 \n%2 %3 ").arg(xDate.toString("MMM yyyy")).arg(barset->label()).arg(barset->at(index), 0, 'f', 1); @@ -2999,7 +2999,7 @@ void Crit3DMeteoWidget::setIsEnsemble(bool value) void Crit3DMeteoWidget::on_actionChangeLeftAxis() { - DialogChangeAxis changeAxisDialog(true); + DialogChangeAxis changeAxisDialog(1, false); if (changeAxisDialog.result() == QDialog::Accepted) { axisY_sx->setMax(changeAxisDialog.getMaxVal()); @@ -3010,7 +3010,7 @@ void Crit3DMeteoWidget::on_actionChangeLeftAxis() void Crit3DMeteoWidget::on_actionChangeRightAxis() { - DialogChangeAxis changeAxisDialog(false); + DialogChangeAxis changeAxisDialog(2, false); if (changeAxisDialog.result() == QDialog::Accepted) { axisY_dx->setMax(changeAxisDialog.getMaxVal()); diff --git a/agrolib/meteoWidget/meteoWidget.h b/agrolib/meteoWidget/meteoWidget.h index 2186dfea..d27e6c62 100644 --- a/agrolib/meteoWidget/meteoWidget.h +++ b/agrolib/meteoWidget/meteoWidget.h @@ -24,7 +24,7 @@ bool getIsEnsemble() { return isEnsemble; } void setNrMembers(int value) { nrMembers = value; } - void setCurrentFrequency(frequencyType frequency); + void setFrequency(frequencyType frequency); void setDailyRange(QDate firstDate, QDate lastDate); void setHourlyRange(QDate firstDate, QDate lastDate); @@ -46,7 +46,7 @@ QVector meteoPointsEnsemble; Crit3DMeteoSettings* meteoSettings; - frequencyType currentFreq; + frequencyType _currentFrequency; QDate firstDailyDate; QDate lastDailyDate; QDate firstHourlyDate; diff --git a/agrolib/outputPoints/dbOutputPointsHandler.cpp b/agrolib/outputPoints/dbOutputPointsHandler.cpp index 3e17f7b1..636add04 100644 --- a/agrolib/outputPoints/dbOutputPointsHandler.cpp +++ b/agrolib/outputPoints/dbOutputPointsHandler.cpp @@ -13,14 +13,10 @@ Crit3DOutputPointsDbHandler::Crit3DOutputPointsDbHandler(QString dbname_) { _db.close(); } - errorString = ""; _db = QSqlDatabase::addDatabase("QSQLITE", QUuid::createUuid().toString()); _db.setDatabaseName(dbname_); - if (! _db.open()) - { - errorString = _db.lastError().text(); - } + _db.open(); } @@ -97,23 +93,26 @@ bool Crit3DOutputPointsDbHandler::addCriteria3DColumn(const QString &tableName, } // column name - QString newField = variableString + "_" + QString::number(depth); + if (depth != NODATA) + { + variableString += "_" + QString::number(depth); + } // column exists already QList fieldList = getFields(&_db, tableName); - if ( fieldList.contains(newField) ) + if ( fieldList.contains(variableString) ) { return true; } // add column QString queryString = "ALTER TABLE '" + tableName + "'"; - queryString += " ADD " + newField + " REAL"; + queryString += " ADD " + variableString + " REAL"; QSqlQuery myQuery = _db.exec(queryString); if (myQuery.lastError().isValid()) { - errorStr = "Error in add column: " + newField + "\n" + myQuery.lastError().text(); + errorStr = "Error in add column: " + variableString + "\n" + myQuery.lastError().text(); return false; } @@ -123,26 +122,78 @@ bool Crit3DOutputPointsDbHandler::addCriteria3DColumn(const QString &tableName, bool Crit3DOutputPointsDbHandler::saveHourlyMeteoData(const QString &tableName, const QDateTime &myTime, const std::vector &varList, - const std::vector &values, + const std::vector &valuesList, QString &errorStr) { - if (varList.size() != values.size()) + if (varList.size() != valuesList.size()) { errorStr = "Error saving values: number of variables is different from values"; return false; } - QSqlQuery qry(_db); QString timeStr = myTime.toString("yyyy-MM-dd HH:mm:ss"); - QString queryString = QString("DELETE FROM '%1' WHERE DATE_TIME ='%2'").arg(tableName, timeStr); + QString queryString = QString("SELECT * FROM '%1' WHERE DATE_TIME='%2'").arg(tableName, timeStr); + QSqlQuery query(_db); + if (! query.exec(queryString)) + { + errorStr = QString("Error in reading table: %1 \nTime: %2 \n%3") + .arg(tableName, timeStr, query.lastError().text()); + return false; + } + + query.last(); + int querySize = query.at() + 1; // SQLITE doesn't support SIZE + if (querySize > 0) + { + return saveHourlyMeteoData_update(tableName, timeStr, varList, valuesList, errorStr); + } + else + { + return saveHourlyMeteoData_insert(tableName, timeStr, varList, valuesList, errorStr); + } + +} + + +bool Crit3DOutputPointsDbHandler::saveHourlyMeteoData_update(const QString &tableName, const QString timeStr, + const std::vector &varList, + const std::vector &values, + QString &errorStr) +{ + QList assignList; + for (unsigned int i = 0; i < varList.size(); i++) + { + QString fieldStr = QString::fromStdString(getMeteoVarName(varList[i])); + if (fieldStr.isEmpty()) + { + errorStr = QString("Error saving values in table:%1 \nMissing variable name").arg(tableName); + return false; + } + + QString valueStr = QString::number(values[i], 'f', 2); + + assignList.push_back("'" + fieldStr + "'=" + valueStr); + } + + QSqlQuery qry(_db); + QString queryString = QString("UPDATE '%1' SET %2 WHERE DATE_TIME = DATETIME('%3')").arg(tableName, assignList.join(','), timeStr); if (! qry.exec(queryString)) { - errorStr = QString("Error deleting values in table:%1 Time:%2\n%3") - .arg(tableName, timeStr, qry.lastError().text()); + errorStr = QString("Error saving values in table:%1 Time:%2\n%3") + .arg(tableName, timeStr, qry.lastError().text()); return false; } + return true; +} + + +bool Crit3DOutputPointsDbHandler::saveHourlyMeteoData_insert(const QString &tableName, const QString timeStr, + const std::vector &varList, + const std::vector &values, + QString &errorStr) +{ // field list QString fieldList = "'DATE_TIME'"; for (unsigned int i = 0; i < varList.size(); i++) @@ -154,23 +205,24 @@ bool Crit3DOutputPointsDbHandler::saveHourlyMeteoData(const QString &tableName, } else { - errorStr = "Error saving values: missing variable name."; + errorStr = QString("Error saving values in table:%1 \nMissing variable name").arg(tableName); return false; } } // values list QString valuesList = "'" + timeStr + "'"; - for (unsigned int i = 0; i < varList.size(); i++) + for (unsigned int i = 0; i < values.size(); i++) { valuesList += "," + QString::number(values[i], 'f', 2); } - queryString = QString("INSERT INTO '%1' (%2) VALUES (%3)").arg(tableName, fieldList, valuesList); + QSqlQuery qry(_db); + QString queryString = QString("INSERT INTO '%1' (%2) VALUES (%3)").arg(tableName, fieldList, valuesList); if (! qry.exec(queryString)) { errorStr = QString("Error saving values in table:%1 Time:%2\n%3") - .arg(tableName, timeStr, qry.lastError().text()); + .arg(tableName, timeStr, qry.lastError().text()); return false; } @@ -178,22 +230,20 @@ bool Crit3DOutputPointsDbHandler::saveHourlyMeteoData(const QString &tableName, } -// layerDepth [m] + +// variableDepth [cm] bool Crit3DOutputPointsDbHandler::saveHourlyCriteria3D_Data(const QString &tableName, const QDateTime& myTime, - const std::vector& varList, - const std::vector& values, - const std::vector & layerDepth, - QString &errorStr) + const std::vector& values, + const std::vector& waterContentDepth, + const std::vector& waterPotentialDepth, + const std::vector& degreeOfSaturationDepth, + const std::vector& factorOfSafetyDepth, + QString &errorStr) { - int nrSoilLayers = int(layerDepth.size()) - 1; - if (nrSoilLayers <= 0) - { - errorStr = "Error saving values: missing soil layers."; - return false; - } + int nrValues = int(waterContentDepth.size() + waterPotentialDepth.size() + + degreeOfSaturationDepth.size() + factorOfSafetyDepth.size()); - int nrValues = int(varList.size()) * nrSoilLayers; - if (nrValues != values.size()) + if (nrValues != int(values.size())) { errorStr = "Error saving values: number of values is not as expected."; return false; @@ -201,35 +251,16 @@ bool Crit3DOutputPointsDbHandler::saveHourlyCriteria3D_Data(const QString &table QString timeStr = myTime.toString("yyyy-MM-dd HH:mm:ss"); - // field list - QString setList = ""; - for (unsigned int i = 0; i < varList.size(); i++) - { - QString variableString = QString::fromStdString(getCriteria3DVarName(varList[i])); - if (variableString == "") - { - errorStr = "Missing variable name."; - return false; - } - - for (int layer = 1; layer <= nrSoilLayers; layer++) - { - int depth = round(layerDepth[layer] * 100); // [cm] + QList valueList; + int firstIndex = 0; + appendCriteria3DOutputValue(volumetricWaterContent, waterContentDepth, values, firstIndex, valueList); + appendCriteria3DOutputValue(waterMatricPotential, waterPotentialDepth, values, firstIndex, valueList); + appendCriteria3DOutputValue(degreeOfSaturation, degreeOfSaturationDepth, values, firstIndex, valueList); + appendCriteria3DOutputValue(factorOfSafety, factorOfSafetyDepth, values, firstIndex, valueList); - QString fieldName = variableString + "_" + QString::number(depth); - setList += "'" + fieldName + "'="; - - int index = i * nrSoilLayers + layer - 1; - QString valueStr = QString::number(values[index], 'f', 3); - setList += valueStr; - - if (index < (nrValues - 1)) - setList += ","; - } - } QSqlQuery qry(_db); - QString queryString = QString("UPDATE '%1' SET %2 WHERE DATE_TIME ='%3'").arg(tableName, setList, timeStr); + QString queryString = QString("UPDATE '%1' SET %2 WHERE DATE_TIME = DATETIME('%3')").arg(tableName, valueList.join(','), timeStr); if (! qry.exec(queryString)) { errorStr = QString("Error in query: " + queryString + "\n" + qry.lastError().text()); @@ -239,3 +270,26 @@ bool Crit3DOutputPointsDbHandler::saveHourlyCriteria3D_Data(const QString &table return true; } + +void Crit3DOutputPointsDbHandler::appendCriteria3DOutputValue(criteria3DVariable myVar, const std::vector &depthList, + const std::vector& values, int &firstIndex, + QList &outputList) +{ + QString variableString = QString::fromStdString(getCriteria3DVarName(myVar)); + + for (int l = 0; l < depthList.size(); l++) + { + float depth_cm = depthList[l]; + QString fieldName = variableString + "_" + QString::number(depth_cm); + + int index = firstIndex + l; + QString valueStr = QString::number(values[index], 'f', 3); + + QString assignStr = "'" + fieldName + "'=" + valueStr; + + outputList.push_back(assignStr); + } + + firstIndex += int(depthList.size()); +} + diff --git a/agrolib/outputPoints/dbOutputPointsHandler.h b/agrolib/outputPoints/dbOutputPointsHandler.h index 00742d1f..7007fbab 100644 --- a/agrolib/outputPoints/dbOutputPointsHandler.h +++ b/agrolib/outputPoints/dbOutputPointsHandler.h @@ -6,6 +6,7 @@ #endif #include + #include #include class Crit3DOutputPointsDbHandler @@ -14,14 +15,14 @@ explicit Crit3DOutputPointsDbHandler(QString dbname_); ~Crit3DOutputPointsDbHandler(); - QString getDbName() { - return _db.databaseName(); } + QString getDbName() + { return _db.databaseName(); } - QString getErrorString() { - return errorString; } + QString getErrorString() + { return _db.lastError().text(); } - bool isOpen() { - return _db.isOpen(); } + bool isOpen() + { return _db.isOpen(); } bool createTable(const QString &tableName, QString &errorStr); @@ -31,16 +32,33 @@ bool saveHourlyMeteoData(const QString &tableName, const QDateTime &myTime, const std::vector &varList, - const std::vector &values, QString &errorStr); + const std::vector &valuesList, QString &errorStr); - bool saveHourlyCriteria3D_Data(const QString &tableName, const QDateTime &myTime, - const std::vector &varList, - const std::vector &values, - const std::vector &layerDepth, QString &errorStr); + bool saveHourlyCriteria3D_Data(const QString &tableName, const QDateTime& myTime, + const std::vector& values, + const std::vector& waterContentDepth, + const std::vector& waterPotentialDepth, + const std::vector& degreeOfSaturationDepth, + const std::vector& factorOfSafetyDepth, + QString &errorStr); + + void appendCriteria3DOutputValue(criteria3DVariable myVar, const std::vector &depthList, + const std::vector& values, int &firstIndex, + QList &outputList); private: + QSqlDatabase _db; - QString errorString; + + bool saveHourlyMeteoData_insert(const QString &tableName, const QString timeStr, + const std::vector &varList, + const std::vector &values, + QString &errorStr); + + bool saveHourlyMeteoData_update(const QString &tableName, const QString timeStr, + const std::vector &varList, + const std::vector &values, + QString &errorStr); }; diff --git a/agrolib/pointStatisticsWidget/pointStatisticsWidget.cpp b/agrolib/pointStatisticsWidget/pointStatisticsWidget.cpp index 3a7fdbc5..2e7f70ae 100644 --- a/agrolib/pointStatisticsWidget/pointStatisticsWidget.cpp +++ b/agrolib/pointStatisticsWidget/pointStatisticsWidget.cpp @@ -1689,7 +1689,7 @@ void Crit3DPointStatisticsWidget::computePlot() void Crit3DPointStatisticsWidget::on_actionChangeLeftAxis() { - DialogChangeAxis changeAxisDialog(true); + DialogChangeAxis changeAxisDialog(1, false); if (changeAxisDialog.result() == QDialog::Accepted) { chartView->setYmax(changeAxisDialog.getMaxVal()); diff --git a/agrolib/project/dialogInterpolation.cpp b/agrolib/project/dialogInterpolation.cpp index 11ae05a7..b78bd227 100644 --- a/agrolib/project/dialogInterpolation.cpp +++ b/agrolib/project/dialogInterpolation.cpp @@ -142,7 +142,7 @@ DialogInterpolation::DialogInterpolation(Project *myProject) std::map::const_iterator itElFunc; for (itElFunc = fittingFunctionNames.begin(); itElFunc != fittingFunctionNames.end(); ++itElFunc) { - if (itElFunc->first == "linear" || itElFunc->first == "Nonlinear Frei function (6 parameters)") + if (itElFunc->first == "linear") continue; elevationFunctionEdit.addItem(QString::fromStdString(itElFunc->first), QString::fromStdString(itElFunc->first)); } diff --git a/agrolib/project/interpolationCmd.cpp b/agrolib/project/interpolationCmd.cpp index d7adce14..97620915 100644 --- a/agrolib/project/interpolationCmd.cpp +++ b/agrolib/project/interpolationCmd.cpp @@ -264,7 +264,10 @@ bool interpolationRaster(std::vector &myPoints, C return false; } - raster.initializeParameters(*raster.header); + if (mySettings->getUseLocalDetrending()) + { + raster.initializeParameters(*raster.header); + } float myX, myY; std::vector proxyValues; diff --git a/agrolib/project/project.cpp b/agrolib/project/project.cpp index 1f311929..dadd9e02 100644 --- a/agrolib/project/project.cpp +++ b/agrolib/project/project.cpp @@ -198,9 +198,10 @@ void Project::setProxyDEM() } -bool Project::checkProxy(const Crit3DProxy &myProxy, QString* error) +bool Project::checkProxy(Crit3DProxy &myProxy, QString* error) { std::string name_ = myProxy.getName(); + QList myList; if (name_ == "") { @@ -216,6 +217,79 @@ bool Project::checkProxy(const Crit3DProxy &myProxy, QString* error) return false; } + if (isHeight) + { + if (parameters->contains("fitting_function")) + { + std::string elevationFuction = parameters->value("fitting_function").toString().toStdString(); + if (fittingFunctionNames.find(elevationFuction) == fittingFunctionNames.end()) + { + errorString = "Unknown function for elevation. Remove the field from the .ini file or choose between: piecewise_two, triple_piecewise, free_triple_piecewise."; + return false; + } + else + myProxy.setFittingFunctionName(fittingFunctionNames.at(elevationFuction)); + + if (parameters->contains("fitting_parameters")) + { + unsigned int nrParameters = NODATA; + + if (myProxy.getFittingFunctionName() == piecewiseTwo) + nrParameters = 4; + else if (myProxy.getFittingFunctionName() == piecewiseThree) + nrParameters = 5; + else if (myProxy.getFittingFunctionName()== piecewiseThreeFree) + nrParameters = 6; + + myList = parameters->value("fitting_parameters").toStringList(); + if (myList.size() != nrParameters*2) + { + *error = "Wrong number of fitting parameters for proxy: " + QString::fromStdString(name_); + return false; + } + + myProxy.setFittingParametersRange(StringListToDouble(myList)); + } + } + else + { + if (parameters->contains("fitting_parameters")) + { + myList = parameters->value("fitting_parameters").toStringList(); + + if (myList.size() == 8) + myProxy.setFittingFunctionName(piecewiseTwo); + else if (myList.size() == 10) + myProxy.setFittingFunctionName(piecewiseThree); + else if (myList.size() == 12) + myProxy.setFittingFunctionName(piecewiseThreeFree); + else + { + *error = "Wrong number of fitting parameters for proxy: " + QString::fromStdString(name_); + return false; + } + myProxy.setFittingParametersRange(StringListToDouble(myList)); + } + } + } + else + { + myProxy.setFittingFunctionName(linear); + if(parameters->contains("fitting_parameters")) + { + unsigned int nrParameters = 2; + + myList = parameters->value("fitting_parameters").toStringList(); + if (myList.size() != nrParameters*2) + { + *error = "Wrong number of fitting parameters for proxy: " + QString::fromStdString(name_); + return false; + } + + myProxy.setFittingParametersRange(StringListToDouble(myList)); + } + } + return true; } @@ -688,98 +762,6 @@ bool Project::loadParameters(QString parametersFileName) if (parameters->contains("stddev_threshold")) myProxy->setStdDevThreshold(parameters->value("stddev_threshold").toFloat()); - /*if (parameters->contains("fitting_parameters")) - { - unsigned int nrParameters; - - if (getProxyPragaName(name_.toStdString()) == proxyHeight) - nrParameters = 5; - else - nrParameters = 2; - - myList = parameters->value("fitting_parameters").toStringList(); - if (myList.size() != nrParameters*2 && myList.size() != (nrParameters-1)*2 && myList.size() != (nrParameters+1)*2) //TODO: change - { - errorString = "Wrong number of fitting parameters for proxy: " + name_; - return false; - } - - myProxy->setFittingParametersRange(StringListToDouble(myList)); - }*/ - - if (getProxyPragaName(name_.toStdString()) == proxyHeight) - { - if (parameters->contains("fitting_function")) - { - std::string elevationFuction = parameters->value("fitting_function").toString().toStdString(); - if (fittingFunctionNames.find(elevationFuction) == fittingFunctionNames.end()) - { - errorString = "Unknown function for elevation. Remove the field from the .ini file or choose between: piecewise_two, triple_piecewise, free_triple_piecewise."; - return false; - } - else - myProxy->setFittingFunctionName(fittingFunctionNames.at(elevationFuction)); - - if (parameters->contains("fitting_parameters")) - { - unsigned int nrParameters; - - if (myProxy->getFittingFunctionName() == piecewiseTwo) - nrParameters = 4; - else if (myProxy->getFittingFunctionName() == piecewiseThree) - nrParameters = 5; - else if (myProxy->getFittingFunctionName()== piecewiseThreeFree) - nrParameters = 6; - - myList = parameters->value("fitting_parameters").toStringList(); - if (myList.size() != nrParameters*2) - { - errorString = "Wrong number of fitting parameters for proxy: " + name_; - return false; - } - - myProxy->setFittingParametersRange(StringListToDouble(myList)); - } - } - else - { - if (parameters->contains("fitting_parameters")) - { - myList = parameters->value("fitting_parameters").toStringList(); - - if (myList.size() == 8) - myProxy->setFittingFunctionName(piecewiseTwo); - else if (myList.size() == 10) - myProxy->setFittingFunctionName(piecewiseThree); - else if (myList.size() == 12) - myProxy->setFittingFunctionName(piecewiseThreeFree); - else - { - errorString = "Wrong number of fitting parameters for proxy: " + name_; - return false; - } - myProxy->setFittingParametersRange(StringListToDouble(myList)); - } - } - } - else - { - myProxy->setFittingFunctionName(linear); - if(parameters->contains("fitting_parameters")) - { - unsigned int nrParameters = 2; - - myList = parameters->value("fitting_parameters").toStringList(); - if (myList.size() != nrParameters*2) - { - errorString = "Wrong number of fitting parameters for proxy: " + name_; - return false; - } - - myProxy->setFittingParametersRange(StringListToDouble(myList)); - } - } - if (! parameters->contains("active")) { errorString = "active not specified for proxy " + QString::fromStdString(myProxy->getName()); @@ -934,8 +916,13 @@ Crit3DTime Project::getCrit3DCurrentTime() QDateTime Project::getCurrentTime() { QDateTime myDateTime; - myDateTime.setDate(this->currentDate); - return myDateTime.addSecs(this->currentHour * HOUR_SECONDS); + if (gisSettings.isUTC) + { + myDateTime.setTimeSpec(Qt::UTC); + } + + myDateTime.setDate(currentDate); + return myDateTime.addSecs(currentHour * HOUR_SECONDS); } @@ -1102,7 +1089,8 @@ bool Project::loadDEM(QString myFileName) setProxyDEM(); interpolationSettings.setProxyLoaded(false); - if (! updateProxy()) return false; + if (! updateProxy()) + return false; //set interpolation settings DEM interpolationSettings.setCurrentDEM(&DEM); @@ -1899,7 +1887,7 @@ bool Project::loadProxyGrids() } else { - errorString = "Error loading proxy grid " + fileName; + errorString = "Error loading raster proxy:\n" + fileName + "\nHow to fix it: check the proxy section in the parameters.ini"; return false; } @@ -1972,13 +1960,11 @@ bool Project::updateProxy() { if (! interpolationSettings.getProxyLoaded()) { - if (loadProxyGrids()) - { - interpolationSettings.setProxyLoaded(true); - } - else + if (! loadProxyGrids()) return false; + interpolationSettings.setProxyLoaded(true); + if (meteoPointsDbHandler != nullptr) { if (! readProxyValues()) @@ -2227,13 +2213,11 @@ bool Project::computeStatisticsCrossValidation(Crit3DTime myTime, meteoVariable std::vector obs; std::vector pre; - float value; - for (int i = 0; i < nrMeteoPoints; i++) { if (meteoPoints[i].active) { - value = meteoPoints[i].getMeteoPointValue(myTime, myVar, meteoSettings); + float value = meteoPoints[i].currentValue; if (! isEqual(value, NODATA) && ! isEqual(meteoPoints[i].residual, NODATA)) { @@ -2278,11 +2262,11 @@ bool Project::interpolationCv(meteoVariable myVar, const Crit3DTime& myTime, cro } if (myVar == dailyGlobalRadiation || + myVar == globalIrradiance || myVar == dailyLeafWetness || myVar == dailyWindVectorDirectionPrevailing || myVar == dailyWindVectorIntensityAvg || - myVar == dailyWindVectorIntensityMax || - myVar == globalIrradiance) + myVar == dailyWindVectorIntensityMax ) { logError("Cross validation is not available for " + QString::fromStdString(getVariableString(myVar))); return false; @@ -2292,7 +2276,7 @@ bool Project::interpolationCv(meteoVariable myVar, const Crit3DTime& myTime, cro std::string errorStdStr; // check quality and pass data to interpolation - if (!checkAndPassDataToInterpolation(quality, myVar, meteoPoints, nrMeteoPoints, myTime, + if (! checkAndPassDataToInterpolation(quality, myVar, meteoPoints, nrMeteoPoints, myTime, &qualityInterpolationSettings, &interpolationSettings, meteoSettings, &climateParameters, interpolationPoints, checkSpatialQuality, errorStdStr)) @@ -2427,7 +2411,7 @@ bool Project::interpolationDemLocalDetrending(meteoVariable myVar, const Crit3DT myRaster->initializeGrid(myHeader); myRaster->initializeParameters(myHeader); - if(!setAllFittingRanges(myCombination, &interpolationSettings)) + if(!setHeightFittingRange(myCombination, &interpolationSettings)) { errorString = "Error in function preInterpolation: \n couldn't set fitting ranges."; return false; @@ -2705,7 +2689,7 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) proxyValues.resize(unsigned(interpolationSettings.getProxyNr())); if (interpolationSettings.getUseLocalDetrending()) - if(!setAllFittingRanges(myCombination, &interpolationSettings)) + if(!setHeightFittingRange(myCombination, &interpolationSettings)) { errorString = "Error in function preInterpolation: \n couldn't set fitting ranges."; return false; @@ -2800,7 +2784,9 @@ float Project::meteoDataConsistency(meteoVariable myVar, const Crit3DTime& timeI { float dataConsistency = 0.0; for (int i = 0; i < nrMeteoPoints; i++) + { dataConsistency = MAXVALUE(dataConsistency, meteoPoints[i].obsDataConsistencyH(myVar, timeIni, timeFin)); + } return dataConsistency; } @@ -2963,6 +2949,7 @@ frequencyType Project::getCurrentFrequency() const return currentFrequency; } + void Project::setCurrentFrequency(const frequencyType &value) { currentFrequency = value; @@ -3220,7 +3207,6 @@ bool Project::loadProject() { errorType = ERROR_SETTINGS; errorString = "Load parameters failed.\n" + errorString; - logError(); return false; } @@ -3362,7 +3348,7 @@ void Project::showMeteoWidgetPoint(std::string idMeteoPoint, std::string namePoi QDateTime lastHourly = meteoPointsDbHandler->getLastDate(hourly, idMeteoPoint); bool hasHourlyData = !(firstHourly.isNull() || lastHourly.isNull()); - if (!hasDailyData && !hasHourlyData) + if (! hasDailyData && ! hasHourlyData) { logInfoGUI("No data."); return; @@ -3430,6 +3416,10 @@ void Project::showMeteoWidgetPoint(std::string idMeteoPoint, std::string namePoi if (hasHourlyData) { meteoWidgetPoint->setHourlyRange(firstHourly.date(), lastHourly.date()); + if (! hasDailyData) + { + meteoWidgetPoint->setFrequency(hourly); + } } meteoWidgetPoint->setCurrentDate(this->currentDate); @@ -3702,12 +3692,12 @@ void Project::showLocalProxyGraph(gis::Crit3DGeoPoint myPoint, gis::Crit3DRaster int row, col; std::vector> parameters; - if (myDataRaster->isLoaded && !myDataRaster->parametersCell.empty()) + if (myDataRaster->isLoaded && !myDataRaster->singleCell.empty()) { gis::getRowColFromXY(*(myDataRaster->header), myUtm, &row, &col); parameters = myDataRaster->getParametersFromRowCol(row, col); } - if (this->meteoGridLoaded && !this->meteoGridDbHandler->meteoGrid()->dataMeteoGrid.parametersCell.empty()) + if (this->meteoGridLoaded && !this->meteoGridDbHandler->meteoGrid()->dataMeteoGrid.singleCell.empty()) { gis::getGridRowColFromLonLat(meteoGridDbHandler->meteoGrid()->gridStructure().header(), myPoint.longitude, myPoint.latitude, &row, &col); parameters = meteoGridDbHandler->meteoGrid()->dataMeteoGrid.getParametersFromRowCol(row, col); @@ -4775,9 +4765,7 @@ void Project::waterTableShowSingleWell(WaterTable &waterTable, const QString &id { DialogSummary* dialogResult = new DialogSummary(waterTable); // show results dialogResult->show(); - WaterTableWidget* chartResult = new WaterTableWidget(idWell, waterTable.getMyDates(), waterTable.getMyHindcastSeries(), - waterTable.getMyInterpolateSeries(), waterTable.getObsDepths(), - quality->getWaterTableMaximumDepth()); + WaterTableWidget* chartResult = new WaterTableWidget(idWell, waterTable, quality->getWaterTableMaximumDepth()); chartResult->show(); return; } diff --git a/agrolib/project/project.h b/agrolib/project/project.h index c8ef74d5..03f37e3c 100644 --- a/agrolib/project/project.h +++ b/agrolib/project/project.h @@ -181,7 +181,7 @@ void setProxyDEM(); void clearProxyDEM(); - bool checkProxy(const Crit3DProxy &myProxy, QString *error); + bool checkProxy(Crit3DProxy &myProxy, QString *error); bool addProxyToProject(std::vector proxyList, std::deque proxyActive, std::vector proxyOrder); void addProxyGridSeries(QString name_, std::vector gridNames, std::vector gridYears); void setCurrentDate(QDate myDate); diff --git a/agrolib/proxyWidget/localProxyWidget.cpp b/agrolib/proxyWidget/localProxyWidget.cpp index 77af178b..181a0180 100644 --- a/agrolib/proxyWidget/localProxyWidget.cpp +++ b/agrolib/proxyWidget/localProxyWidget.cpp @@ -1,34 +1,9 @@ -/*! - CRITERIA3D - \copyright 2016 Fausto Tomei, Gabriele Antolini, Laura Costantini - Alberto Pistocchi, Marco Bittelli, Antonio Volta - You should have received a copy of the GNU General Public License - along with Nome-Programma. If not, see . - This file is part of CRITERIA3D. - CRITERIA3D has been developed under contract issued by A.R.P.A. Emilia-Romagna - CRITERIA3D is free software: you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - CRITERIA3D is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - You should have received a copy of the /NU Lesser General Public License - along with CRITERIA3D. If not, see . - contacts: - fausto.tomei@gmail.com - ftomei@arpae.it -*/ - #include "meteo.h" #include "localProxyWidget.h" -#include "proxyWidget.h" #include "utilities.h" #include "interpolation.h" #include "spatialControl.h" #include "commonConstants.h" -#include "formInfo.h" #include "math.h" #include "furtherMathFunctions.h" @@ -62,7 +37,12 @@ Crit3DLocalProxyWidget::Crit3DLocalProxyWidget(double x, double y, std::vectoraddWidget(&detrended); selectionOptionBoxLayout->addWidget(&modelLR); selectionOptionBoxLayout->addWidget(&climatologicalLR); + selectionOptionBoxLayout->addWidget(&stationWeights); selectionOptionEditLayout->addWidget(r2Label); selectionOptionEditLayout->addWidget(&r2); @@ -140,25 +121,7 @@ Crit3DLocalProxyWidget::Crit3DLocalProxyWidget(double x, double y, std::vectoraddStretch(30); selectionLayout->addLayout(selectionOptionLayout); - if (!parameters.empty() && interpolationSettings->getProxy(proxyPos)->getFittingFunctionName() == piecewiseThree && parameters[proxyPos].size() == 5) - { - QVBoxLayout *parametriLayout = new QVBoxLayout(); - - QLabel *H0Lab = new QLabel(QString("H0: %1").arg(parameters[proxyPos][0])); - QLabel *T0Lab = new QLabel(QString("T0: %1").arg(parameters[proxyPos][1])); - QLabel *H1Lab = new QLabel(QString("H1: %1").arg(parameters[proxyPos][0]+parameters[proxyPos][2])); - QLabel *T1Lab = new QLabel(QString("T1: %1").arg(parameters[proxyPos][1]+parameters[proxyPos][3])); - QLabel *slopeLab = new QLabel(QString("slope: %1").arg(parameters[proxyPos][4])); - - parametriLayout->addWidget(H0Lab); - parametriLayout->addWidget(T0Lab); - parametriLayout->addWidget(H1Lab); - parametriLayout->addWidget(T1Lab); - parametriLayout->addWidget(slopeLab); - - selectionLayout->addLayout(parametriLayout); - } - else if (!parameters.empty() && interpolationSettings->getProxy(proxyPos)->getFittingFunctionName() == piecewiseTwo && parameters[proxyPos].size() == 4) + if (!parameters.empty() && interpolationSettings->getProxy(proxyPos)->getFittingFunctionName() == piecewiseTwo && parameters[proxyPos].size() == 4) { QVBoxLayout *parametriLayout = new QVBoxLayout(); @@ -242,6 +205,7 @@ Crit3DLocalProxyWidget::Crit3DLocalProxyWidget(double x, double y, std::vectorclimatologicalLRClicked(toggled); }); connect(&modelLR, &QCheckBox::toggled, [=](int toggled){ this->modelLRClicked(toggled); }); connect(&detrended, &QCheckBox::toggled, [=](){ this->plot(); }); + connect(&stationWeights, &QCheckBox::toggled, [=] () {this->plot();}); connect(updateStations, &QAction::triggered, this, [=](){ this->plot(); }); if (currentFrequency != noFrequency) @@ -371,6 +335,14 @@ void Crit3DLocalProxyWidget::plot() outInterpolationPoints.clear(); subsetInterpolationPoints.clear(); std::string errorStdStr; + + for (QGraphicsTextItem* label : weightLabels) + { + chartView->scene()->removeItem(label); + delete label; + } + weightLabels.clear(); + if (detrended.isChecked()) { outInterpolationPoints.clear(); @@ -407,7 +379,7 @@ void Crit3DLocalProxyWidget::plot() point.setY(varValue); QString text = "id: " + QString::fromStdString(meteoPoints[subsetInterpolationPoints[i].index].id) + "\n" + "name: " + QString::fromStdString(meteoPoints[subsetInterpolationPoints[i].index].name) + "\n" - + "weight: " + QString::number(subsetInterpolationPoints[i].regressionWeight); + + "weight: " + QString::number(subsetInterpolationPoints[i].regressionWeight, 'f', 5); if (subsetInterpolationPoints[i].isMarked) { pointListMarked.append(point); @@ -485,6 +457,37 @@ void Crit3DLocalProxyWidget::plot() { modelLRClicked(1); } + + if (stationWeights.isChecked()) + { + QChart* chart = chartView->chart(); + QRectF chartRect = chart->plotArea(); + double xMin = chartView->axisX->min(); + double xMax = chartView->axisX->max(); + double yMin = chartView->axisY->min(); + double yMax = chartView->axisY->max(); + + for (int i = 0; i < int(subsetInterpolationPoints.size()); i++) + { + float proxyVal = subsetInterpolationPoints[i].getProxyValue(proxyPos); + float varValue = subsetInterpolationPoints[i].value; + + if (proxyVal != NODATA && varValue != NODATA) + { + double xRatio = (proxyVal - xMin) / (xMax - xMin); + double yRatio = (varValue - yMin) / (yMax - yMin); + + QPointF scenePos; + scenePos.setX(chartRect.left() + xRatio * chartRect.width()); + scenePos.setY(chartRect.bottom() - yRatio * chartRect.height()); + + QGraphicsTextItem* weightLabel = new QGraphicsTextItem(QString::number(subsetInterpolationPoints[i].regressionWeight, 'f', 3)); + weightLabel->setPos(scenePos); + chartView->scene()->addItem(weightLabel); + weightLabels.push_back(weightLabel); + } + } + } } @@ -544,7 +547,7 @@ void Crit3DLocalProxyWidget::modelLRClicked(int toggled) for (int p = 0; p < int(xVector.size()); p++) { point.setX(xVector[p]); - point.setY(lapseRatePiecewiseThree_withSlope(xVector[p], parameters[proxyPos])); + point.setY(lapseRatePiecewise_three(xVector[p], parameters[proxyPos])); point_vector.append(point); } } @@ -581,7 +584,7 @@ void Crit3DLocalProxyWidget::modelLRClicked(int toggled) for (int p = 0; p < int(xVector.size()); p++) { point.setX(xVector[p]); - point.setY(lapseRatePiecewiseFree(xVector[p], parameters[proxyPos])); + point.setY(lapseRatePiecewise_three_free(xVector[p], parameters[proxyPos])); point_vector.append(point); } diff --git a/agrolib/proxyWidget/localProxyWidget.h b/agrolib/proxyWidget/localProxyWidget.h index cb62df4f..63a291ba 100644 --- a/agrolib/proxyWidget/localProxyWidget.h +++ b/agrolib/proxyWidget/localProxyWidget.h @@ -23,7 +23,6 @@ class Crit3DLocalProxyWidget : public QWidget void plot(); void climatologicalLRClicked(int toggled); void modelLRClicked(int toggled); - private: double x; double y; @@ -49,11 +48,13 @@ class Crit3DLocalProxyWidget : public QWidget QCheckBox detrended; QCheckBox climatologicalLR; QCheckBox modelLR; + QCheckBox stationWeights; QTextEdit r2; QTextEdit lapseRate; ChartView *chartView; meteoVariable myVar; int proxyPos; + std::vector weightLabels; Crit3DTime getCurrentTime(); diff --git a/agrolib/proxyWidget/proxyWidget.cpp b/agrolib/proxyWidget/proxyWidget.cpp index 1961b967..8f0d0e55 100644 --- a/agrolib/proxyWidget/proxyWidget.cpp +++ b/agrolib/proxyWidget/proxyWidget.cpp @@ -1,33 +1,9 @@ -/*! - CRITERIA3D - \copyright 2016 Fausto Tomei, Gabriele Antolini, Laura Costantini - Alberto Pistocchi, Marco Bittelli, Antonio Volta - You should have received a copy of the GNU General Public License - along with Nome-Programma. If not, see . - This file is part of CRITERIA3D. - CRITERIA3D has been developed under contract issued by A.R.P.A. Emilia-Romagna - CRITERIA3D is free software: you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - CRITERIA3D is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - You should have received a copy of the /NU Lesser General Public License - along with CRITERIA3D. If not, see . - contacts: - fausto.tomei@gmail.com - ftomei@arpae.it -*/ - #include "meteo.h" #include "proxyWidget.h" #include "utilities.h" #include "interpolation.h" #include "spatialControl.h" #include "commonConstants.h" -#include "formInfo.h" #include "math.h" #include diff --git a/agrolib/soilFluxes3D/header/soilFluxes3D.h b/agrolib/soilFluxes3D/header/soilFluxes3D.h index 7341310c..9c8ebd07 100644 --- a/agrolib/soilFluxes3D/header/soilFluxes3D.h +++ b/agrolib/soilFluxes3D/header/soilFluxes3D.h @@ -49,8 +49,9 @@ __EXTERN int DLL_EXPORT __STDCALL setNodeSurface(long nodeIndex, int surfaceIndex); // WATER - __EXTERN int DLL_EXPORT __STDCALL setHydraulicProperties(int waterRetentionCurve, int conductivityMeanType, float horizVertRatioConductivity); + __EXTERN int DLL_EXPORT __STDCALL setHydraulicProperties(int waterRetentionCurve, int conductivityMeanType, float conductivityHorizVertRatio); __EXTERN int DLL_EXPORT __STDCALL setWaterContent(long index, double myWaterContent); + __EXTERN int DLL_EXPORT __STDCALL setDegreeOfSaturation(long nodeIndex, double degreeOfSaturation); __EXTERN int DLL_EXPORT __STDCALL setMatricPotential(long index, double psi); __EXTERN int DLL_EXPORT __STDCALL setTotalPotential(long index, double totalPotential); __EXTERN int DLL_EXPORT __STDCALL setPrescribedTotalPotential(long index, double prescribedTotalPotential); diff --git a/agrolib/soilFluxes3D/header/solver.h b/agrolib/soilFluxes3D/header/solver.h index ba8739c2..52195724 100644 --- a/agrolib/soilFluxes3D/header/solver.h +++ b/agrolib/soilFluxes3D/header/solver.h @@ -1,7 +1,7 @@ #ifndef SOLVER_H #define SOLVER_H - inline double square(double x) {return ((x)*(x));} + inline double square(double x) { return ((x)*(x)); } double distance(unsigned long index1, unsigned long index2); diff --git a/agrolib/soilFluxes3D/soilFluxes3D.cpp b/agrolib/soilFluxes3D/soilFluxes3D.cpp index 33d9fc19..4cb717a9 100644 --- a/agrolib/soilFluxes3D/soilFluxes3D.cpp +++ b/agrolib/soilFluxes3D/soilFluxes3D.cpp @@ -182,18 +182,18 @@ int DLL_EXPORT __STDCALL setNumericalParameters(float minDeltaT, float maxDeltaT * k_lateral_vertical_ratio = 10 * \param waterRetentionCurve * \param conductivityMeanType - * \param horizVertRatioConductivity + * \param conductivityHorizVertRatio * \return OK or PARAMETER_ERROR */ int DLL_EXPORT __STDCALL setHydraulicProperties(int waterRetentionCurve, - int conductivityMeanType, float horizVertRatioConductivity) + int conductivityMeanType, float conductivityHorizVertRatio) { myParameters.waterRetentionCurve = waterRetentionCurve; myParameters.meanType = conductivityMeanType; - if ((horizVertRatioConductivity >= 0.1) && (horizVertRatioConductivity <= 100)) + if ((conductivityHorizVertRatio >= 0.1) && (conductivityHorizVertRatio <= 100)) { - myParameters.k_lateral_vertical_ratio = horizVertRatioConductivity; + myParameters.k_lateral_vertical_ratio = conductivityHorizVertRatio; return CRIT3D_OK; } else @@ -499,7 +499,7 @@ int DLL_EXPORT __STDCALL setHydraulicProperties(int waterRetentionCurve, /*! * \brief setWaterContent * \param nodeIndex - * \param waterContent [m^3 m^-3] + * \param waterContent [m] surface - [m3 m-3] sub-surface * \return OK/ERROR */ int DLL_EXPORT __STDCALL setWaterContent(long nodeIndex, double waterContent) @@ -531,6 +531,31 @@ int DLL_EXPORT __STDCALL setHydraulicProperties(int waterRetentionCurve, } + /*! + * \brief setDegreeOfSaturation + * \param nodeIndex + * \param degreeOfSaturation [-] (only sub-surface) + * \return OK/ERROR + */ + int DLL_EXPORT __STDCALL setDegreeOfSaturation(long nodeIndex, double degreeOfSaturation) + { + if (nodeListPtr == nullptr) return MEMORY_ERROR; + + if ((nodeIndex < 0) || (nodeIndex >= myStructure.nrNodes)) return INDEX_ERROR; + + if (nodeListPtr[nodeIndex].isSurface) return INDEX_ERROR; + + if ((degreeOfSaturation < 0.) || (degreeOfSaturation > 1.)) return PARAMETER_ERROR; + + nodeListPtr[nodeIndex].Se = degreeOfSaturation; + nodeListPtr[nodeIndex].H = nodeListPtr[nodeIndex].z - psi_from_Se(nodeIndex); + nodeListPtr[nodeIndex].oldH = nodeListPtr[nodeIndex].H; + nodeListPtr[nodeIndex].k = computeK(nodeIndex); + + return CRIT3D_OK; + } + + /*! * \brief setWaterSinkSource * \param nodeIndex diff --git a/agrolib/synchronicityWidget/synchronicityWidget.cpp b/agrolib/synchronicityWidget/synchronicityWidget.cpp index 277f9d47..663ed37f 100644 --- a/agrolib/synchronicityWidget/synchronicityWidget.cpp +++ b/agrolib/synchronicityWidget/synchronicityWidget.cpp @@ -398,7 +398,7 @@ void Crit3DSynchronicityWidget::changeInterpolationDate() void Crit3DSynchronicityWidget::on_actionChangeLeftSynchAxis() { - DialogChangeAxis changeAxisDialog(true); + DialogChangeAxis changeAxisDialog(1, false); if (changeAxisDialog.result() == QDialog::Accepted) { synchronicityChartView->setYmax(changeAxisDialog.getMaxVal()); @@ -408,7 +408,7 @@ void Crit3DSynchronicityWidget::on_actionChangeLeftSynchAxis() void Crit3DSynchronicityWidget::on_actionChangeLeftInterpolationAxis() { - DialogChangeAxis changeAxisDialog(true); + DialogChangeAxis changeAxisDialog(1, false); if (changeAxisDialog.result() == QDialog::Accepted) { interpolationChartView->setYmax(changeAxisDialog.getMaxVal()); diff --git a/agrolib/waterTable/waterTable.cpp b/agrolib/waterTable/waterTable.cpp index 12eafbf0..c96e1001 100644 --- a/agrolib/waterTable/waterTable.cpp +++ b/agrolib/waterTable/waterTable.cpp @@ -108,8 +108,7 @@ bool WaterTable::computeWTClimate() H_num.push_back(0); } - QMap myDepths = well.getObsDepths(); - QMapIterator it(myDepths); + QMapIterator it(well.depths); while (it.hasNext()) { it.next(); @@ -216,7 +215,6 @@ bool WaterTable::computeETP_allSeries(bool isUpdateAvgCWB) // Ricerca del periodo di correlazione migliore bool WaterTable::computeCWBCorrelation(int stepDays) { - QMap myDepths = well.getObsDepths(); std::vector myCWBSum; std::vector myObsWT; float a = NODATA; @@ -232,7 +230,7 @@ bool WaterTable::computeCWBCorrelation(int stepDays) { myCWBSum.clear(); myObsWT.clear(); - QMapIterator it(myDepths); + QMapIterator it(well.depths); while (it.hasNext()) { @@ -314,8 +312,7 @@ double WaterTable::computeCWB(QDate myDate, int nrDays) // function to compute several statistical indices for watertable depth bool WaterTable::computeWaterTableIndices() { - QMap myDepths = well.getObsDepths(); - QMapIterator it(myDepths); + QMapIterator it(well.depths); std::vector myObs; std::vector myComputed; std::vector myClimate; @@ -491,8 +488,7 @@ bool WaterTable::getWaterTableInterpolation(QDate myDate, float* myValue, float* QDate nextDate; int dT; - QMap myDepths = well.getObsDepths(); - QList keys = myDepths.keys(); + QList keys = well.depths.keys(); // check previuos and next observed data int lastIndex = keys.size() - 1; @@ -501,10 +497,10 @@ bool WaterTable::getWaterTableInterpolation(QDate myDate, float* myValue, float* { indexPrev = i; previousDate = keys[indexPrev]; - previosValue = myDepths[previousDate]; + previosValue = well.depths[previousDate]; indexNext = i; nextDate = keys[indexNext]; - nextValue = myDepths[nextDate]; + nextValue = well.depths[nextDate]; } else { @@ -512,13 +508,13 @@ bool WaterTable::getWaterTableInterpolation(QDate myDate, float* myValue, float* { indexNext = 0; nextDate = keys[indexNext]; - nextValue = myDepths[nextDate]; + nextValue = well.depths[nextDate]; } else if (keys[lastIndex] < myDate) { indexPrev = lastIndex; previousDate = keys[indexPrev]; - previosValue = myDepths[previousDate]; + previosValue = well.depths[previousDate]; } else { @@ -527,10 +523,10 @@ bool WaterTable::getWaterTableInterpolation(QDate myDate, float* myValue, float* { indexPrev = i; previousDate = keys[indexPrev]; - previosValue = myDepths[previousDate]; + previosValue = well.depths[previousDate]; indexNext = i + 1; nextDate = keys[indexNext]; - nextValue = myDepths[nextDate]; + nextValue = well.depths[nextDate]; break; } } @@ -612,30 +608,29 @@ bool WaterTable::getWaterTableInterpolation(QDate myDate, float* myValue, float* void WaterTable::computeWaterTableSeries() { - myDates.clear(); - myHindcastSeries.clear(); - myInterpolateSeries.clear(); + hindcastSeries.clear(); + interpolationSeries.clear(); + float myDelta; int myDeltaDays; - QDate firstObsDate = std::min(well.getFirstDate(), firstMeteoDate); - int numValues = firstObsDate.daysTo(lastMeteoDate) + 1; + firstDate = std::min(well.getFirstDate(), firstMeteoDate); + int numValues = firstDate.daysTo(lastMeteoDate) + 1; for (int i = 0; i < numValues; i++) { - QDate currentDate = firstObsDate.addDays(i); - myDates.push_back(currentDate); + QDate currentDate = firstDate.addDays(i); float currentDepth = getWaterTableDaily(currentDate); - myHindcastSeries.push_back(currentDepth); + hindcastSeries.push_back(currentDepth); if (getWaterTableInterpolation(currentDate, ¤tDepth, &myDelta, &myDeltaDays)) { - myInterpolateSeries.push_back(currentDepth); + interpolationSeries.push_back(currentDepth); } else { - myInterpolateSeries.push_back(NODATA); + interpolationSeries.push_back(NODATA); } } } diff --git a/agrolib/waterTable/waterTable.h b/agrolib/waterTable/waterTable.h index b31075ee..f6c9ece7 100644 --- a/agrolib/waterTable/waterTable.h +++ b/agrolib/waterTable/waterTable.h @@ -18,6 +18,13 @@ class WaterTable { public: + float WTClimateMonthly[12]; + float WTClimateDaily[366]; + + QDate firstDate; + std::vector hindcastSeries; + std::vector interpolationSeries; + WaterTable(std::vector &inputTMin, std::vector &inputTMax, std::vector &inputPrec, QDate firstMeteoDate, QDate lastMeteoDate, Crit3DMeteoSettings meteoSettings); @@ -59,13 +66,7 @@ class WaterTable QDate getFirstDateWell() { return well.getFirstDate(); } QDate getLastDateWell() { return well.getLastDate(); } - QMap getObsDepths() { return well.getObsDepths(); } - - std::vector getMyDates() { return myDates; } - - std::vector getMyHindcastSeries() { return myHindcastSeries; } - - std::vector getMyInterpolateSeries() { return myInterpolateSeries; } + Well* getWell() { return &well; } void cleanAllMeteoVector(); @@ -92,17 +93,10 @@ class WaterTable float EF; bool isClimateReady; - float WTClimateMonthly[12]; - float WTClimateDaily[366]; int nrObsData; bool isCWBEquationReady; double avgDailyCWB; //[mm] - - // graph - std::vector myDates; - std::vector myHindcastSeries; - std::vector myInterpolateSeries; }; #endif // WATERTABLE_H diff --git a/agrolib/waterTable/waterTable.pro b/agrolib/waterTable/waterTable.pro index 55a69b5a..069634f6 100644 --- a/agrolib/waterTable/waterTable.pro +++ b/agrolib/waterTable/waterTable.pro @@ -4,7 +4,7 @@ # #------------------------------------------------- -QT += widgets charts core xml +QT += core widgets charts QT -= gui TEMPLATE = lib @@ -24,7 +24,7 @@ win32:{ TARGET = waterTable } -INCLUDEPATH += ../crit3dDate ../mathFunctions ../meteo ../interpolation ../gis ../commonChartElements +INCLUDEPATH += ../mathFunctions ../crit3dDate ../meteo ../gis ../commonChartElements SOURCES += \ dialogSelectWell.cpp \ diff --git a/agrolib/waterTable/waterTableChartView.cpp b/agrolib/waterTable/waterTableChartView.cpp index 30891b63..ffd5d83c 100644 --- a/agrolib/waterTable/waterTableChartView.cpp +++ b/agrolib/waterTable/waterTableChartView.cpp @@ -1,12 +1,14 @@ #include "waterTableChartView.h" + WaterTableChartView::WaterTableChartView(QWidget *parent) : QChartView(new QChart(), parent) { obsDepthSeries = new QScatterSeries(); obsDepthSeries->setName("Observed"); - obsDepthSeries->setColor(Qt::green); - obsDepthSeries->setMarkerSize(8.0); + obsDepthSeries->setColor(QColor(0, 192, 255)); + obsDepthSeries->setBorderColor(QColor(0,0,1)); + obsDepthSeries->setMarkerSize(4.0); hindcastSeries = new QLineSeries(); hindcastSeries->setName("hindcast"); @@ -16,6 +18,13 @@ WaterTableChartView::WaterTableChartView(QWidget *parent) : interpolationSeries->setName("interpolation"); interpolationSeries->setColor(QColor(0,0,1)); + QPen pen; + pen.setWidth(2); + climateSeries = new QLineSeries(); + climateSeries->setName("climate"); + climateSeries->setPen(pen); + climateSeries->setColor(QColor(0, 200, 0, 255)); + axisX = new QDateTimeAxis(); axisX->setFormat("yyyy/MM"); axisY = new QValueAxis(); @@ -32,40 +41,47 @@ WaterTableChartView::WaterTableChartView(QWidget *parent) : } -void WaterTableChartView::draw(std::vector &myDates, std::vector &myHindcastSeries, std::vector &myInterpolateSeries, - QMap obsDepths, float maximumObservedDepth) +void WaterTableChartView::drawWaterTable(WaterTable &waterTable, float maximumObservedDepth) { axisY->setMax(maximumObservedDepth); // unit of observed watertable data, usually [cm] axisY->setMin(0); axisY->setLabelFormat("%d"); axisY->setTickCount(16); - QDateTime firstDate; - firstDate.setDate(myDates[0]); - QDateTime lastDate; - lastDate.setDate(myDates[myDates.size()-1]); + QDateTime firstDateTime, lastDateTime; + int nrDays = int(waterTable.interpolationSeries.size()); + firstDateTime.setDate(waterTable.firstDate); + lastDateTime.setDate(waterTable.firstDate); + lastDateTime = lastDateTime.addDays(nrDays-1); - axisX->setRange(firstDate, lastDate); + axisX->setRange(firstDateTime, lastDateTime); axisX->setTickCount(15); - int nDays = int(myDates.size()); - QDateTime currentDateTime; - for (int day = 0; day < nDays; day++) + QDateTime currentDateTime = firstDateTime; + for (int day = 0; day < nrDays; day++) { - currentDateTime.setDate(myDates[day]); - hindcastSeries->append(currentDateTime.toMSecsSinceEpoch(), myHindcastSeries[day]); - interpolationSeries->append(currentDateTime.toMSecsSinceEpoch(), myInterpolateSeries[day]); + QDate firstJanuary; + firstJanuary.setDate(currentDateTime.date().year(), 1, 1); + int doyIndex = firstJanuary.daysTo(currentDateTime.date()); // from 0 to 365 - if(obsDepths.contains(myDates[day])) + hindcastSeries->append(currentDateTime.toMSecsSinceEpoch(), waterTable.hindcastSeries[day]); + interpolationSeries->append(currentDateTime.toMSecsSinceEpoch(), waterTable.interpolationSeries[day]); + climateSeries->append(currentDateTime.toMSecsSinceEpoch(), waterTable.WTClimateDaily[doyIndex]); + + if(waterTable.getWell()->depths.contains(currentDateTime.date())) { - int myDepth = obsDepths[myDates[day]]; + int myDepth = waterTable.getWell()->depths[currentDateTime.date()]; obsDepthSeries->append(currentDateTime.toMSecsSinceEpoch(), myDepth); } + + currentDateTime = currentDateTime.addDays(1); } - chart()->addSeries(obsDepthSeries); + chart()->addSeries(hindcastSeries); + chart()->addSeries(climateSeries); chart()->addSeries(interpolationSeries); + chart()->addSeries(obsDepthSeries); obsDepthSeries->attachAxis(axisX); obsDepthSeries->attachAxis(axisY); @@ -73,10 +89,13 @@ void WaterTableChartView::draw(std::vector &myDates, std::vector & hindcastSeries->attachAxis(axisY); interpolationSeries->attachAxis(axisX); interpolationSeries->attachAxis(axisY); + climateSeries->attachAxis(axisX); + climateSeries->attachAxis(axisY); connect(obsDepthSeries, &QScatterSeries::hovered, this, &WaterTableChartView::tooltipObsDepthSeries); connect(hindcastSeries, &QLineSeries::hovered, this, &WaterTableChartView::tooltipLineSeries); connect(interpolationSeries, &QLineSeries::hovered, this, &WaterTableChartView::tooltipLineSeries); + foreach(QLegendMarker* marker, chart()->legend()->markers()) { marker->setVisible(true); diff --git a/agrolib/waterTable/waterTableChartView.h b/agrolib/waterTable/waterTableChartView.h index 879bc998..278f3b23 100644 --- a/agrolib/waterTable/waterTableChartView.h +++ b/agrolib/waterTable/waterTableChartView.h @@ -11,20 +11,21 @@ public: WaterTableChartView(QWidget *parent = 0); - void draw(std::vector &myDates, std::vector &myHindcastSeries, std::vector &myInterpolateSeries, - QMap obsDepths, float maximumObservedDepth); + void drawWaterTable(WaterTable &waterTable, float maximumObservedDepth); void tooltipObsDepthSeries(QPointF point, bool state); void tooltipLineSeries(QPointF point, bool state); void handleMarkerClicked(); QList exportInterpolationValues(); + QDateTimeAxis* axisX; + private: QScatterSeries* obsDepthSeries; QLineSeries* hindcastSeries; QLineSeries* interpolationSeries; + QLineSeries* climateSeries; - QDateTimeAxis* axisX; QValueAxis* axisY; Callout *m_tooltip; }; diff --git a/agrolib/waterTable/waterTableWidget.cpp b/agrolib/waterTable/waterTableWidget.cpp index b8e4bec8..de70e5dc 100644 --- a/agrolib/waterTable/waterTableWidget.cpp +++ b/agrolib/waterTable/waterTableWidget.cpp @@ -1,7 +1,7 @@ #include "waterTableWidget.h" +#include "dialogChangeAxis.h" -WaterTableWidget::WaterTableWidget(const QString &id, std::vector myDates, std::vector myHindcastSeries, - std::vector myInterpolateSeries, QMap obsDepths, float maxObservedDepth) +WaterTableWidget::WaterTableWidget(const QString &id, WaterTable &waterTable, float maxObservedDepth) { this->setWindowTitle("Graph Id well: "+ id); this->resize(1240, 700); @@ -21,20 +21,38 @@ WaterTableWidget::WaterTableWidget(const QString &id, std::vector myDates menuBar->addMenu(editMenu); mainLayout->setMenuBar(menuBar); - QAction* exportInterpolation = new QAction(tr("&Export interpolation as csv"), this); + QAction* exportInterpolation = new QAction(tr("&Export interpolation as csv.."), this); + QAction* changeXAxis = new QAction(tr("&Change period (X axis).."), this); editMenu->addAction(exportInterpolation); + editMenu->addAction(changeXAxis); mainLayout->addLayout(plotLayout); setLayout(mainLayout); connect(exportInterpolation, &QAction::triggered, this, &WaterTableWidget::on_actionExportInterpolationData); + connect(changeXAxis, &QAction::triggered, this, &WaterTableWidget::on_actionChangeXAxis); - waterTableChartView->draw(myDates, myHindcastSeries, myInterpolateSeries, obsDepths, maxObservedDepth); + waterTableChartView->drawWaterTable(waterTable, maxObservedDepth); } -void WaterTableWidget::on_actionExportInterpolationData() + +void WaterTableWidget::on_actionChangeXAxis() { + DialogChangeAxis changeAxisDialog(0, true); + + if (changeAxisDialog.result() == QDialog::Accepted) + { + QDateTime firstDate, lastDate; + firstDate.setDate(changeAxisDialog.getMinDate()); + lastDate.setDate(changeAxisDialog.getMaxDate()); + waterTableChartView->axisX->setRange(firstDate, lastDate); + waterTableChartView->update(); + } +} + +void WaterTableWidget::on_actionExportInterpolationData() +{ QString csvFileName = QFileDialog::getSaveFileName(this, tr("Save current data"), "", tr("csv files (*.csv)")); if (csvFileName != "") { @@ -66,12 +84,3 @@ void WaterTableWidget::on_actionExportInterpolationData() } } -WaterTableWidget::~WaterTableWidget() -{ - -} - -void WaterTableWidget::closeEvent(QCloseEvent *event) -{ - event->accept(); -} diff --git a/agrolib/waterTable/waterTableWidget.h b/agrolib/waterTable/waterTableWidget.h index cb43494d..ac8f440d 100644 --- a/agrolib/waterTable/waterTableWidget.h +++ b/agrolib/waterTable/waterTableWidget.h @@ -10,15 +10,19 @@ { Q_OBJECT public: - WaterTableWidget(const QString &id, std::vector myDates, std::vector myHindcastSeries, - std::vector myInterpolateSeries, QMap obsDepths, float maxObservedDepth); - ~WaterTableWidget(); - void closeEvent(QCloseEvent *event); - void on_actionExportInterpolationData(); + WaterTableWidget(const QString &id, WaterTable &waterTable, float maxObservedDepth); + + ~WaterTableWidget() { ; } + + void closeEvent(QCloseEvent *event) + { event->accept(); } private: WaterTableChartView *waterTableChartView; + void on_actionExportInterpolationData(); + void on_actionChangeXAxis(); + }; #endif // WATERTABLEWIDGET_H diff --git a/agrolib/waterTable/well.cpp b/agrolib/waterTable/well.cpp index be4b33c9..963a303b 100644 --- a/agrolib/waterTable/well.cpp +++ b/agrolib/waterTable/well.cpp @@ -25,11 +25,6 @@ int Well::getObsDepthNr() } -QMap Well::getObsDepths() -{ - return depths; -} - QDate Well::getFirstDate() { diff --git a/agrolib/waterTable/well.h b/agrolib/waterTable/well.h index aeb72800..9d26ca77 100644 --- a/agrolib/waterTable/well.h +++ b/agrolib/waterTable/well.h @@ -9,6 +9,8 @@ class Well { public: + QMap depths; + Well(); QString getId() const { return id; } @@ -32,9 +34,6 @@ class Well QDate getLastDate(); int getObsDepthNr(); - - QMap getObsDepths(); - int minValuesPerMonth(); private: @@ -43,8 +42,6 @@ class Well double utmX, utmY; double lat, lon; - QMap depths; - QDate firstDate; QDate lastDate; };