From d09b84771e027e0784030155555dc3fcd6e2118c Mon Sep 17 00:00:00 2001 From: avolta Date: Fri, 18 Oct 2024 16:28:39 +0200 Subject: [PATCH 01/34] aggiunta commento in r2 generalizzato --- mathFunctions/furtherMathFunctions.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mathFunctions/furtherMathFunctions.cpp b/mathFunctions/furtherMathFunctions.cpp index e2c6447ab..e004c151e 100644 --- a/mathFunctions/furtherMathFunctions.cpp +++ b/mathFunctions/furtherMathFunctions.cpp @@ -1272,7 +1272,7 @@ namespace interpolation double computeWeighted_R2_generalized(const std::vector& observed, const std::vector& predicted, const std::vector>& weights) { - + // !!! still to be checked // This function computes the generalized weighted R-squared (coefficient of determination) int n = int(observed.size()); double sum_weighted_squared_total = 0.; From 0d0a20c92a6c25e2d111f041477c1aba36d741f9 Mon Sep 17 00:00:00 2001 From: Caterina Toscano Date: Tue, 22 Oct 2024 11:25:20 +0200 Subject: [PATCH 02/34] added glocal detrending option for DEM interpolation --- interpolation/interpolation.cpp | 89 +++++- interpolation/interpolation.h | 3 +- interpolation/interpolationPoint.cpp | 1 + interpolation/interpolationPoint.h | 1 + interpolation/interpolationSettings.cpp | 40 +++ interpolation/interpolationSettings.h | 12 + interpolation/spatialControl.cpp | 1 + meteo/meteoPoint.cpp | 1 + meteo/meteoPoint.h | 1 + pragaProject/pragaProject.cpp | 7 + project/dialogInterpolation.cpp | 24 ++ project/dialogInterpolation.h | 2 + project/project.cpp | 357 +++++++++++++++++++++++- project/project.h | 4 + proxyWidget/localProxyWidget.cpp | 40 +++ 15 files changed, 578 insertions(+), 5 deletions(-) diff --git a/interpolation/interpolation.cpp b/interpolation/interpolation.cpp index e479828c6..f5e7522af 100644 --- a/interpolation/interpolation.cpp +++ b/interpolation/interpolation.cpp @@ -1727,7 +1727,7 @@ bool multipleDetrendingMain(std::vector &myPoints elevationPos = pos; } - if (mySettings->getCurrentCombination().isProxyActive(elevationPos)) + if (elevationPos != NODATA && mySettings->getCurrentCombination().isProxyActive(elevationPos)) { if (!multipleDetrendingElevationFitting(elevationPos, myPoints, mySettings, myVar, errorStr)) return false; @@ -2061,6 +2061,82 @@ void detrendingOtherProxies(int elevationPos, std::vector &myPoints, Crit3DInterpolationSettings* mySettings, meteoVariable myVar, std::string& errorStr) +{ + //per n aree, si creano gruppi di stazioni in base al codice (sorta di selection) + //poi si chiama fitting per ogni area, il fitting va salvato in una matrice + //matrice avrà, per ogni area, i parametri di fitting + la combination + il valore delle stazioni detrendate secondo il fitting corrispondente + std::vector subsetPoints; + int i = 0; + int j = 0; + int areasNr = int(myPoints.front().macroAreaCode); + int start = areasNr; + + //look for highest and lowest "macroAreaCode", that's the areas we'll have to iterate the fitting for + for (i = 0; i < myPoints.size(); i++) + { + if (myPoints[i].macroAreaCode > areasNr) + areasNr = int(myPoints[i].macroAreaCode); + if (myPoints[i].macroAreaCode < start && !(myPoints[i].macroAreaCode == NODATA)) + start = int(myPoints[i].macroAreaCode); + } + + std::vector>> allAreaParameters; + std::vector allAreaCombinations; + //now create the subset of points, including only the right points, and then call multipleDetrending + for (i = start-1; i < areasNr+1; i++) //ATTENZIONE INDICI + { + + for (j = 0; j < myPoints.size(); j++) + { + if (myPoints[j].macroAreaCode == i) + subsetPoints.push_back(myPoints[j]); + } + + if (! subsetPoints.empty()) + { + mySettings->setCurrentCombination(mySettings->getSelectedCombination()); + mySettings->clearFitting(); + + int elevationPos = NODATA; + for (unsigned int pos=0; pos < mySettings->getCurrentCombination().getProxySize(); pos++) + { + if (getProxyPragaName(mySettings->getProxy(pos)->getName()) == proxyHeight) + elevationPos = pos; //SPOSTA + } + + if (elevationPos != NODATA && mySettings->getCurrentCombination().isProxyActive(elevationPos)) + { + if (!multipleDetrendingElevationFitting(elevationPos, subsetPoints, mySettings, myVar, errorStr)) return false; + + detrendingElevation(elevationPos, subsetPoints, mySettings); + + } + if (!multipleDetrendingOtherProxiesFitting(elevationPos, subsetPoints, mySettings, myVar, errorStr)) + return false; + + detrendingOtherProxies(elevationPos, subsetPoints, mySettings); + + //salvataggio parametri e valori staz(?) parlane con gabri pk dovresti salvare tutti i crit3dintpoint(?) + allAreaParameters.push_back(mySettings->getFittingParameters()); + allAreaCombinations.push_back(mySettings->getCurrentCombination()); + } + else + { + std::vector> emptyParameters; + Crit3DProxyCombination emptyCombination; + allAreaParameters.push_back(emptyParameters); + allAreaCombinations.push_back(emptyCombination); + } + subsetPoints.clear(); + } + + mySettings->setAreaParameters(allAreaParameters); + mySettings->setAreaCombination(allAreaCombinations); + + return true; +} + void topographicDistanceOptimize(meteoVariable myVar, Crit3DMeteoPoint* &myMeteoPoints, @@ -2169,8 +2245,15 @@ bool preInterpolation(std::vector &myPoints, Crit { if (!mySettings->getUseLocalDetrending()) setMultipleDetrendingHeightTemperatureRange(mySettings); - - if (! multipleDetrendingMain(myPoints, mySettings, myVar, errorStr)) return false; + + if (mySettings->getUseGlocalDetrending()) + { + glocalDetrendingFitting(myPoints, mySettings, myVar, errorStr); + } + else + { + if (! multipleDetrendingMain(myPoints, mySettings, myVar, errorStr)) return false; + } } else { diff --git a/interpolation/interpolation.h b/interpolation/interpolation.h index e48348e18..63b6c6a22 100644 --- a/interpolation/interpolation.h +++ b/interpolation/interpolation.h @@ -79,7 +79,8 @@ void detrendingElevation(int elevationPos, std::vector &myPoints, Crit3DInterpolationSettings* mySettings); void detrendingOtherProxies(int elevationPos, std::vector &myPoints, Crit3DInterpolationSettings* mySettings); - + + bool glocalDetrendingFitting(std::vector &myPoints, Crit3DInterpolationSettings* mySettings, meteoVariable myVar, std::string& errorStr); bool getUseDetrendingVar(meteoVariable myVar); bool isThermal(meteoVariable myVar); bool getUseTdVar(meteoVariable myVar); diff --git a/interpolation/interpolationPoint.cpp b/interpolation/interpolationPoint.cpp index c456f72fd..5e9aaff39 100644 --- a/interpolation/interpolationPoint.cpp +++ b/interpolation/interpolationPoint.cpp @@ -39,6 +39,7 @@ Crit3DInterpolationDataPoint::Crit3DInterpolationDataPoint() regressionWeight = NODATA; lapseRateCode = primary; + macroAreaCode = NODATA; topographicDistance = nullptr; diff --git a/interpolation/interpolationPoint.h b/interpolation/interpolationPoint.h index 82caaef55..cd641a803 100644 --- a/interpolation/interpolationPoint.h +++ b/interpolation/interpolationPoint.h @@ -25,6 +25,7 @@ lapseRateCodeType lapseRateCode; gis::Crit3DRasterGrid* topographicDistance; std::vector proxyValues; + int macroAreaCode; Crit3DInterpolationDataPoint(); diff --git a/interpolation/interpolationSettings.cpp b/interpolation/interpolationSettings.cpp index c7cafa4df..51ffb12c5 100644 --- a/interpolation/interpolationSettings.cpp +++ b/interpolation/interpolationSettings.cpp @@ -439,10 +439,12 @@ void Crit3DInterpolationSettings::initializeProxy() void Crit3DInterpolationSettings::initialize() { currentDEM = nullptr; + macroAreasMap = nullptr; interpolationMethod = idw; useThermalInversion = true; useTD = false; useLocalDetrending = false; + useGlocalDetrending = false; topoDist_maxKh = 128; useDewPoint = true; useInterpolatedTForRH = true; @@ -504,6 +506,38 @@ std::string getKeyStringElevationFunction(TFittingFunction value) return key; } +void Crit3DInterpolationSettings::setMacroAreasMap(gis::Crit3DRasterGrid *value) +{ + macroAreasMap = value; +} + +gis::Crit3DRasterGrid* Crit3DInterpolationSettings::getMacroAreasMap() +{ + return macroAreasMap; +} + +std::vector>> Crit3DInterpolationSettings::getAreaParameters() +{ + return areaParameters; +} + +void Crit3DInterpolationSettings::setAreaParameters(std::vector > > myAreaParameters) +{ + areaParameters.clear(); + areaParameters = myAreaParameters; +} + +std::vector Crit3DInterpolationSettings::getAreaCombination() +{ + return areaCombination; +} + +void Crit3DInterpolationSettings::setAreaCombination(std::vector myAreaCombination) +{ + areaCombination.clear(); + areaCombination = myAreaCombination; +} + TInterpolationMethod Crit3DInterpolationSettings::getInterpolationMethod() { return interpolationMethod;} @@ -513,6 +547,9 @@ bool Crit3DInterpolationSettings::getUseTD() bool Crit3DInterpolationSettings::getUseLocalDetrending() { return useLocalDetrending;} +bool Crit3DInterpolationSettings::getUseGlocalDetrending() +{ return useGlocalDetrending;} + bool Crit3DInterpolationSettings::getUseDoNotRetrend() { return useDoNotRetrend;} @@ -537,6 +574,9 @@ void Crit3DInterpolationSettings::setUseTD(bool myValue) void Crit3DInterpolationSettings::setUseLocalDetrending(bool myValue) { useLocalDetrending = myValue;} +void Crit3DInterpolationSettings::setUseGlocalDetrending(bool myValue) +{ useGlocalDetrending = myValue;} + void Crit3DInterpolationSettings::setUseDoNotRetrend(bool myValue) { useDoNotRetrend = myValue;} diff --git a/interpolation/interpolationSettings.h b/interpolation/interpolationSettings.h index c21801eb9..42357ab54 100644 --- a/interpolation/interpolationSettings.h +++ b/interpolation/interpolationSettings.h @@ -137,6 +137,7 @@ { private: gis::Crit3DRasterGrid* currentDEM; //for TD + gis::Crit3DRasterGrid* macroAreasMap; //for glocal detrending TInterpolationMethod interpolationMethod; @@ -144,6 +145,7 @@ bool useThermalInversion; bool useTD; bool useLocalDetrending; + bool useGlocalDetrending; int maxTdMultiplier; bool useLapseRateCode; bool useBestDetrending; @@ -175,6 +177,8 @@ std::vector > fittingParameters; std::vector&)>> fittingFunction; std::vector pointsRange; + std::vector>> areaParameters; + std::vector areaCombination; public: @@ -203,6 +207,14 @@ void setUseLocalDetrending(bool myValue); bool getUseLocalDetrending(); + void setUseGlocalDetrending(bool myValue); + bool getUseGlocalDetrending(); + void setMacroAreasMap(gis::Crit3DRasterGrid *value); + gis::Crit3DRasterGrid *getMacroAreasMap(); + std::vector > > getAreaParameters(); + void setAreaParameters(std::vector>> myAreaParameters); + std::vector getAreaCombination(); + void setAreaCombination(std::vector myAreaCombination); void setUseDoNotRetrend(bool myValue); bool getUseDoNotRetrend(); diff --git a/interpolation/spatialControl.cpp b/interpolation/spatialControl.cpp index c450e51da..5f9a9d532 100644 --- a/interpolation/spatialControl.cpp +++ b/interpolation/spatialControl.cpp @@ -429,6 +429,7 @@ bool passDataToInterpolation(Crit3DMeteoPoint* meteoPoints, int nrMeteoPoints, myPoint.topographicDistance = meteoPoints[i].topographicDistance; myPoint.isActive = true; myPoint.isMarked = meteoPoints[i].marked; + myPoint.macroAreaCode = meteoPoints[i].macroAreaCode; if (isEqual(xMin, NODATA)) { diff --git a/meteo/meteoPoint.cpp b/meteo/meteoPoint.cpp index e32f7f411..97ae29fc4 100644 --- a/meteo/meteoPoint.cpp +++ b/meteo/meteoPoint.cpp @@ -61,6 +61,7 @@ void Crit3DMeteoPoint::clear() this->latInt = NODATA; this->lonInt = NODATA; this->isInsideDem = false; + this->macroAreaCode = NODATA; this->nrObsDataDaysH = 0; this->nrObsDataDaysD = 0; diff --git a/meteo/meteoPoint.h b/meteo/meteoPoint.h index 2c86a6da4..dfbde836f 100644 --- a/meteo/meteoPoint.h +++ b/meteo/meteoPoint.h @@ -91,6 +91,7 @@ int latInt; int lonInt; bool isInsideDem; + int macroAreaCode; bool isUTC; bool isForecast; diff --git a/pragaProject/pragaProject.cpp b/pragaProject/pragaProject.cpp index 27b4aeecc..9472a24fe 100644 --- a/pragaProject/pragaProject.cpp +++ b/pragaProject/pragaProject.cpp @@ -2636,6 +2636,13 @@ bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QL return false; } + if (interpolationSettings.getUseGlocalDetrending()) + { + logInfoGUI("Loading macro areas map for glocal detrending..."); + if (! loadMacroAreaGlocalMap(false)) + return false; + } + // save also time aggregated variables foreach (myVar, aggrVariables) varToSave.push_back(myVar); diff --git a/project/dialogInterpolation.cpp b/project/dialogInterpolation.cpp index e95a194f1..8b4151a93 100644 --- a/project/dialogInterpolation.cpp +++ b/project/dialogInterpolation.cpp @@ -128,6 +128,11 @@ DialogInterpolation::DialogInterpolation(Project *myProject) connect(localDetrendingEdit, SIGNAL(stateChanged(int)), this, SLOT(localDetrendingChanged(int))); + glocalDetrendingEdit = new QCheckBox(tr("glocal detrending")); + glocalDetrendingEdit->setChecked(_interpolationSettings->getUseGlocalDetrending()); + + connect(glocalDetrendingEdit, SIGNAL(stateChanged(int)), this, SLOT(glocalDetrendingChanged(int))); + doNotRetrendEdit = new QCheckBox(tr("do not retrend")); doNotRetrendEdit->setChecked(_interpolationSettings->getUseDoNotRetrend()); @@ -143,6 +148,7 @@ DialogInterpolation::DialogInterpolation(Project *myProject) layoutDetrending->addWidget(&minPointsLocalDetrendingEdit); layoutDetrending->addWidget(localDetrendingEdit); + layoutDetrending->addWidget(glocalDetrendingEdit); layoutDetrending->addWidget(doNotRetrendEdit); layoutDetrending->addWidget(retrendOnlyEdit); @@ -237,12 +243,29 @@ void DialogInterpolation::localDetrendingChanged(int active) topographicDistanceEdit->setEnabled(active == Qt::Unchecked); maxTdMultiplierEdit.setEnabled(active == Qt::Unchecked); minPointsLocalDetrendingEdit.setEnabled(active == Qt::Checked); + if (active == Qt::Checked) optimalDetrendingEdit->setChecked(Qt::Unchecked); + optimalDetrendingEdit->setEnabled(active == Qt::Unchecked); + if (active == Qt::Checked) glocalDetrendingEdit->setChecked(Qt::Unchecked); + glocalDetrendingEdit->setEnabled(active == Qt::Unchecked); +} + +void DialogInterpolation::glocalDetrendingChanged(int active) +{ + if (active == Qt::Checked) topographicDistanceEdit->setChecked(Qt::Unchecked); + topographicDistanceEdit->setEnabled(active == Qt::Unchecked); + maxTdMultiplierEdit.setEnabled(active == Qt::Unchecked); + if (active == Qt::Checked) optimalDetrendingEdit->setChecked(Qt::Unchecked); + optimalDetrendingEdit->setEnabled(active == Qt::Unchecked); + if (active == Qt::Checked) localDetrendingEdit->setChecked(Qt::Unchecked); + localDetrendingEdit->setEnabled(active == Qt::Unchecked); } void DialogInterpolation::optimalDetrendingChanged(int active) { if (active == Qt::Checked) localDetrendingEdit->setChecked(Qt::Unchecked); localDetrendingEdit->setEnabled(active == Qt::Unchecked); + if (active == Qt::Checked) glocalDetrendingEdit->setChecked(Qt::Unchecked); + glocalDetrendingEdit->setEnabled(active == Qt::Unchecked); } void DialogInterpolation::redrawProxies() @@ -313,6 +336,7 @@ void DialogInterpolation::accept() _interpolationSettings->setMeteoGridUpscaleFromDem(upscaleFromDemEdit->isChecked()); _interpolationSettings->setUseTD(topographicDistanceEdit->isChecked()); _interpolationSettings->setUseLocalDetrending(localDetrendingEdit->isChecked()); + _interpolationSettings->setUseGlocalDetrending(glocalDetrendingEdit->isChecked()); _interpolationSettings->setUseLapseRateCode(lapseRateCodeEdit->isChecked()); _interpolationSettings->setUseBestDetrending(optimalDetrendingEdit->isChecked()); _interpolationSettings->setUseMultipleDetrending(multipleDetrendingEdit->isChecked()); diff --git a/project/dialogInterpolation.h b/project/dialogInterpolation.h index f708ecfcf..4a26141da 100644 --- a/project/dialogInterpolation.h +++ b/project/dialogInterpolation.h @@ -29,6 +29,7 @@ QCheckBox* multipleDetrendingEdit; QCheckBox* topographicDistanceEdit; QCheckBox* localDetrendingEdit; + QCheckBox* glocalDetrendingEdit; QCheckBox* doNotRetrendEdit; QCheckBox* retrendOnlyEdit; QCheckBox* upscaleFromDemEdit; @@ -52,6 +53,7 @@ void upscaleFromDemChanged(int active); void multipleDetrendingChanged(int active); void localDetrendingChanged(int active); + void glocalDetrendingChanged(int active); void optimalDetrendingChanged(int active); }; diff --git a/project/project.cpp b/project/project.cpp index 9a4672cd7..65f01c26b 100644 --- a/project/project.cpp +++ b/project/project.cpp @@ -18,6 +18,7 @@ #include "dialogSummary.h" #include "waterTableWidget.h" #include "utilities.h" +#include "furtherMathFunctions.h" #include @@ -100,6 +101,7 @@ void Project::initializeProject() dbGridXMLFileName = ""; outputPointsFileName = ""; currentDbOutputFileName = ""; + glocalMapName = ""; meteoPointsLoaded = false; meteoGridLoaded = false; @@ -626,6 +628,12 @@ bool Project::loadParameters(QString parametersFileName) if (parametersSettings->contains("localDetrending")) interpolationSettings.setUseLocalDetrending(parametersSettings->value("localDetrending").toBool()); + + if (parametersSettings->contains("glocalDetrending")) + interpolationSettings.setUseGlocalDetrending(parametersSettings->value("glocalDetrending").toBool()); + + if (parametersSettings->contains("glocalMapName")) + this->glocalMapName = parametersSettings->value("glocalMapName").toString(); if (parametersSettings->contains("meteogrid_upscalefromdem")) interpolationSettings.setMeteoGridUpscaleFromDem(parametersSettings->value("meteogrid_upscalefromdem").toBool()); @@ -2134,6 +2142,89 @@ bool Project::loadTopographicDistanceMaps(bool onlyWithData, bool showInfo) return true; } +bool Project::loadMacroAreaGlocalMap(bool showInfo) +{ + //TODO: add a check for the code values? + + if (nrMeteoPoints == 0) + { + logError("Open a meteo points DB before."); + return false; + } + + QString mapsFolder = defaultPath + PATH_GEO; + if (! QDir(mapsFolder).exists()) + { + //errore? + } + + if (showInfo) + { + //QString infoStr = "Loading map of macro areas for glocal detrending..."; + //logInfo(infoStr); + } + + std::string myError; + gis::Crit3DRasterGrid* macroAreasGrid = new gis::Crit3DRasterGrid(); + std::string fileName = mapsFolder.toStdString() + glocalMapName.toStdString(); + + if (!QFile::exists(QString::fromStdString(fileName + ".flt"))) + { + myError = "Couldn't find " + fileName + " file in " + mapsFolder.toStdString(); + logError(QString::fromStdString(myError)); + return false; + } + + if (gis::readEsriGrid(fileName, macroAreasGrid, myError)) + interpolationSettings.setMacroAreasMap(macroAreasGrid); + else + { + myError = "Error loading:\n" + fileName; + return false; + } + + + //associa il codice a ogni stazione + for (int i=0; i < nrMeteoPoints; i++) + { + int myValue = gis::getValueFromXY(*macroAreasGrid, meteoPoints[i].point.utm.x, meteoPoints[i].point.utm.y); + + if (myValue != (*macroAreasGrid).header->flag) + meteoPoints[i].macroAreaCode = myValue; + } + + //if (showInfo) closeLogInfo(); + + return true; +} + +bool Project::groupCellsInArea(std::vector &areaPoints, unsigned int index) +{ + int zoneNr; + double myX, myY; + gis::Crit3DRasterGrid* macroAreas = interpolationSettings.getMacroAreasMap(); + int nrCols = DEM.header->nrCols; + + //DEM. GRID ANCORA DA SCRIVERE + for (int i = 0; i < DEM.header->nrRows; i++) + { + for (int j = 0; j < DEM.header->nrCols; j++) + { + if (!isEqual(DEM.getValueFromRowCol(i,j),DEM.header->flag)) + { + gis::getUtmXYFromRowCol(DEM.header, i, j, &myX, &myY); + + zoneNr = int(macroAreas->getValueFromXY(myX, myY)); + + if (zoneNr == index) + areaPoints.push_back(i*nrCols+j); + } + } + } + + return true; +} + void Project::passInterpolatedTemperatureToHumidityPoints(Crit3DTime myTime, Crit3DMeteoSettings* meteoSettings) { if (! hourlyMeteoMaps->mapHourlyTair->isLoaded) return; @@ -2474,6 +2565,137 @@ bool Project::interpolationDemLocalDetrending(meteoVariable myVar, const Crit3DT return true; } +bool Project::interpolationDemGlocalDetrending(meteoVariable myVar, const Crit3DTime& myTime, gis::Crit3DRasterGrid *myRaster) +{ + if (!getUseDetrendingVar(myVar) || !interpolationSettings.getUseGlocalDetrending()) + return false; + + // pass data to interpolation + std::vector interpolationPoints; + std::string errorStdStr; + + if (!checkAndPassDataToInterpolation(quality, myVar, meteoPoints, nrMeteoPoints, myTime, + &qualityInterpolationSettings, &interpolationSettings, meteoSettings, &climateParameters, interpolationPoints, + checkSpatialQuality, errorStdStr)) + { + errorString = "No data available: " + QString::fromStdString(getVariableString(myVar)) + + "\n" + QString::fromStdString(errorStdStr); + return false; + } + + std::vector proxyValues; + proxyValues.resize(unsigned(interpolationSettings.getProxyNr())); + double x, y; + + Crit3DProxyCombination myCombination = interpolationSettings.getSelectedCombination(); + interpolationSettings.setCurrentCombination(myCombination); + + int elevationPos = NODATA; + for (unsigned int pos=0; pos < interpolationSettings.getCurrentCombination().getProxySize(); pos++) + { + if (getProxyPragaName(interpolationSettings.getProxy(pos)->getName()) == proxyHeight) + elevationPos = pos; + } + + //obtain fitting parameters for every macro area + if (! preInterpolation(interpolationPoints, &interpolationSettings, meteoSettings, &climateParameters, + meteoPoints, nrMeteoPoints, myVar, myTime, errorStdStr)) + { + errorString = "Error in function preInterpolation:\n" + QString::fromStdString(errorStdStr); + return false; + } + + if (getComputeOnlyPoints()) + { + //TODO + } + else + { + gis::Crit3DRasterHeader myHeader = *(DEM.header); + myRaster->initializeGrid(myHeader); + + if(!setMultipleDetrendingHeightTemperatureRange(&interpolationSettings)) + { + errorString = "Error in function preInterpolation: \n couldn't set temperature ranges for height proxy."; + return false; + } + + //preparing to start a "global" interpolation for every single macro area + unsigned row, col; + std::vector>> allAreaFittingParameters = interpolationSettings.getAreaParameters(); + std::vector allAreaCombinations = interpolationSettings.getAreaCombination(); + std::vector subsetInterpolationPoints; + std::vector areaCells; + unsigned numAreas = allAreaFittingParameters.size(); + + if (allAreaFittingParameters.size() < numAreas || allAreaCombinations.size() < numAreas) + { + errorString = "Error in function interpolationGrid: there is something wrong with the macro areas raster file."; + return false; + } + + for (unsigned areaIndex = 0; areaIndex < numAreas; areaIndex++) + { + //save all the cells in that specific area in a vector (areaCells) + areaCells.clear(); + groupCellsInArea(areaCells, areaIndex); + + //take the parameters+combination for that area + std::vector &)>> fittingFunction; + interpolationSettings.setFittingParameters(allAreaFittingParameters[areaIndex]); + interpolationSettings.setCurrentCombination(allAreaCombinations[areaIndex]); + + //find the fitting functions vector based on the length of the parameters vector for every proxy + for (int l = 0; l < allAreaFittingParameters[areaIndex].size(); l++) + { + if (allAreaFittingParameters[areaIndex][l].size() == 2) + fittingFunction.push_back(functionLinear); + else if (allAreaFittingParameters[areaIndex][l].size() == 4) + fittingFunction.push_back(lapseRatePiecewise_two); + else if (allAreaFittingParameters[areaIndex][l].size() == 5) + fittingFunction.push_back(lapseRatePiecewise_three); + else if (allAreaFittingParameters[areaIndex][l].size() == 6) + fittingFunction.push_back(lapseRatePiecewise_three_free); + } + + //create group of macro area interpolation points + interpolationSettings.setFittingFunction(fittingFunction); + subsetInterpolationPoints.clear(); + for (int l = 0; l < interpolationPoints.size(); l++) + { + if (interpolationPoints[l].macroAreaCode == areaIndex) + subsetInterpolationPoints.push_back(interpolationPoints[l]); + } + + //detrending (done once for every macro area) + if (elevationPos != NODATA && myCombination.isProxyActive(elevationPos)) + detrendingElevation(elevationPos, subsetInterpolationPoints, &interpolationSettings); + + detrendingOtherProxies(elevationPos, subsetInterpolationPoints, &interpolationSettings); + + for (unsigned cellIndex = 0; cellIndex < areaCells.size(); cellIndex++) + { + + row = unsigned(areaCells[cellIndex]/DEM.header->nrCols); + col = areaCells[cellIndex]%DEM.header->nrCols; + + float z = DEM.value[row][col]; + if (! isEqual(z, myHeader.flag)) + { + gis::getUtmXYFromRowCol(myHeader, row, col, &x, &y); + + getProxyValuesXY(x, y, &interpolationSettings, proxyValues); + + myRaster->value[row][col] = interpolate(subsetInterpolationPoints, &interpolationSettings, meteoSettings, + myVar, x, y, z, proxyValues, true); + } + } + } + } + + return true; +} + bool Project::interpolateDemRadiation(const Crit3DTime& myTime, gis::Crit3DRasterGrid *myRaster) { @@ -2620,10 +2842,18 @@ bool Project::interpolationDemMain(meteoVariable myVar, const Crit3DTime& myTime if (interpolationSettings.getUseMultipleDetrending()) interpolationSettings.clearFitting(); + if (interpolationSettings.getUseGlocalDetrending()) + { + if (!loadMacroAreaGlocalMap(true)) return false; + } // dynamic lapserate if (getUseDetrendingVar(myVar) && interpolationSettings.getUseLocalDetrending()) { return interpolationDemLocalDetrending(myVar, myTime, myRaster); + } + else if (getUseDetrendingVar(myVar) && interpolationSettings.getUseGlocalDetrending()) + { + return interpolationDemGlocalDetrending(myVar, myTime, myRaster); } else { @@ -2670,6 +2900,13 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) if (interpolationSettings.getUseMultipleDetrending()) interpolationSettings.clearFitting(); + unsigned numAreas; + if (interpolationSettings.getUseGlocalDetrending()) + { + if (!loadMacroAreaGlocalMap(true)) return false; + numAreas = interpolationSettings.getAreaParameters().size(); + } + // check quality and pass data to interpolation if (! checkAndPassDataToInterpolation(quality, myVar, meteoPoints, nrMeteoPoints, myTime, &qualityInterpolationSettings, &interpolationSettings, meteoSettings, &climateParameters, interpolationPoints, @@ -2716,6 +2953,8 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) float interpolatedValue = NODATA; unsigned int i, proxyIndex; + if (!interpolationSettings.getUseGlocalDetrending()) + { for (unsigned col = 0; col < unsigned(meteoGridDbHandler->meteoGrid()->gridStructure().header().nrCols); col++) { for (unsigned row = 0; row < unsigned(meteoGridDbHandler->meteoGrid()->gridStructure().header().nrRows); row++) @@ -2788,6 +3027,120 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) } } } + } + else + { + unsigned row, col; + std::vector>> allAreaFittingParameters = interpolationSettings.getAreaParameters(); + std::vector allAreaCombinations = interpolationSettings.getAreaCombination(); + std::vector areaCells; + + if (allAreaCombinations.size() != numAreas) + { + errorString = "Error in function interpolationGrid: there is something wrong with the macro areas raster file."; + return false; + } + + for (unsigned areaIndex = 0; areaIndex < numAreas; areaIndex++) + { + areaCells.clear(); + groupCellsInArea(areaCells, areaIndex); + std::vector &)> > fittingFunction; + interpolationSettings.setFittingParameters(allAreaFittingParameters[areaIndex]); + interpolationSettings.setCurrentCombination(allAreaCombinations[areaIndex]); + + int elevationPos = NODATA; + for (unsigned int pos=0; pos < interpolationSettings.getCurrentCombination().getProxySize(); pos++) + { + if (getProxyPragaName(interpolationSettings.getProxy(pos)->getName()) == proxyHeight) + elevationPos = pos; + } + + for (int l = 0; l < allAreaFittingParameters[areaIndex].size(); l++) + { + if (allAreaFittingParameters[areaIndex][l].size() == 2) + fittingFunction.push_back(functionLinear); + else if (allAreaFittingParameters[areaIndex][l].size() == 4) + fittingFunction.push_back(lapseRatePiecewise_two); + else if (allAreaFittingParameters[areaIndex][l].size() == 5) + fittingFunction.push_back(lapseRatePiecewise_three); + else if (allAreaFittingParameters[areaIndex][l].size() == 6) + fittingFunction.push_back(lapseRatePiecewise_three_free); + } + + interpolationSettings.setFittingFunction(fittingFunction); + + std::vector subsetInterpolationPoints; + + for (int l = 0; l < interpolationPoints.size(); l++) + { + if (interpolationPoints[l].macroAreaCode == areaIndex) + subsetInterpolationPoints.push_back(interpolationPoints[l]); + } + + if (elevationPos != NODATA && myCombination.isProxyActive(elevationPos)) + detrendingElevation(elevationPos, subsetInterpolationPoints, &interpolationSettings); + + detrendingOtherProxies(elevationPos, subsetInterpolationPoints, &interpolationSettings); + + unsigned nrCols = meteoGridDbHandler->meteoGrid()->gridStructure().header().nrCols; + + for (unsigned cellIndex = 0; cellIndex < areaCells.size(); cellIndex++) + { + row = unsigned(areaCells[cellIndex]/nrCols); + col = areaCells[cellIndex]%nrCols; + + + if (meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->active) + { + myX = meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->point.utm.x; + myY = meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->point.utm.y; + myZ = meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->point.z; + + if (getUseDetrendingVar(myVar)) + { + proxyIndex = 0; + + for (i=0; i < interpolationSettings.getProxyNr(); i++) + { + proxyValues[i] = NODATA; + + if (myCombination.isProxyActive(i)) + { + if (proxyIndex < meteoGridProxies.size()) + { + float proxyValue = gis::getValueFromXY(*meteoGridProxies[proxyIndex], myX, myY); + if (proxyValue != meteoGridProxies[proxyIndex]->header->flag) + proxyValues[i] = double(proxyValue); + } + + proxyIndex++; + } + } + interpolatedValue = interpolate(subsetInterpolationPoints, &interpolationSettings, meteoSettings, myVar, myX, myY, myZ, proxyValues, true); + + } + if (freq == hourly) + { + if (meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->nrObsDataDaysH == 0) + meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->initializeObsDataH(1, 1, myTime.date); + + meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->setMeteoPointValueH(myTime.date, myTime.getHour(), myTime.getMinutes(), myVar, float(interpolatedValue)); + meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->currentValue = float(interpolatedValue); + } + else if (freq == daily) + { + if (meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->nrObsDataDaysD == 0) + meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->initializeObsDataD(1, myTime.date); + + meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->setMeteoPointValueD(myTime.date, myVar, float(interpolatedValue)); + meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->currentValue = float(interpolatedValue); + + } + } + } + } + } return true; } @@ -3052,6 +3405,7 @@ void Project::saveInterpolationParameters() parametersSettings->setValue("topographicDistanceMaxMultiplier", QString::number(interpolationSettings.getTopoDist_maxKh())); parametersSettings->setValue("optimalDetrending", interpolationSettings.getUseBestDetrending()); parametersSettings->setValue("multipleDetrending", interpolationSettings.getUseMultipleDetrending()); + parametersSettings->setValue("glocalDetrending", interpolationSettings.getUseGlocalDetrending()); parametersSettings->setValue("doNotRetrend", interpolationSettings.getUseDoNotRetrend()); parametersSettings->setValue("retrendOnly", interpolationSettings.getUseRetrendOnly()); parametersSettings->setValue("useDewPoint", interpolationSettings.getUseDewPoint()); @@ -3726,7 +4080,8 @@ void Project::showLocalProxyGraph(gis::Crit3DGeoPoint myPoint) myZDEM = DEM.value[row][col]; } - + if (interpolationSettings.getUseGlocalDetrending()) + loadMacroAreaGlocalMap(false); localProxyWidget = new Crit3DLocalProxyWidget(myUtm.x, myUtm.y, myZDEM, myZGrid, this->gisSettings, &interpolationSettings, meteoPoints, nrMeteoPoints, currentVariable, currentFrequency, currentDate, currentHour, quality, &qualityInterpolationSettings, meteoSettings, &climateParameters, checkSpatialQuality); return; diff --git a/project/project.h b/project/project.h index 858f99ab2..7b217b0c9 100644 --- a/project/project.h +++ b/project/project.h @@ -104,6 +104,7 @@ QString aggregationPath; QString dbGridXMLFileName; QString parametersFileName; + QString glocalMapName; std::ofstream logFile; // output points @@ -259,6 +260,8 @@ bool writeTopographicDistanceMap(int pointIndex, const gis::Crit3DRasterGrid& demMap, QString pathTd); bool loadTopographicDistanceMaps(bool onlyWithData, bool showInfo); void passInterpolatedTemperatureToHumidityPoints(Crit3DTime myTime, Crit3DMeteoSettings *meteoSettings); + bool loadMacroAreaGlocalMap(bool showInfo); + bool groupCellsInArea(std::vector &areaPoints, unsigned index); bool checkInterpolation(meteoVariable myVar); bool checkInterpolationGrid(meteoVariable myVar); @@ -266,6 +269,7 @@ bool interpolationDemMain(meteoVariable myVar, const Crit3DTime& myTime, gis::Crit3DRasterGrid *myRaster); bool interpolationDem(meteoVariable myVar, const Crit3DTime& myTime, gis::Crit3DRasterGrid *myRaster); bool interpolationDemLocalDetrending(meteoVariable myVar, const Crit3DTime& myTime, gis::Crit3DRasterGrid *myRaster); + bool interpolationDemGlocalDetrending(meteoVariable myVar, const Crit3DTime& myTime, gis::Crit3DRasterGrid *myRaster); bool interpolateDemRadiation(const Crit3DTime& myTime, gis::Crit3DRasterGrid *myRaster); bool interpolationOutputPoints(std::vector &interpolationPoints, gis::Crit3DRasterGrid *outputGrid, meteoVariable myVar); diff --git a/proxyWidget/localProxyWidget.cpp b/proxyWidget/localProxyWidget.cpp index 6c5b5b831..9870cf168 100644 --- a/proxyWidget/localProxyWidget.cpp +++ b/proxyWidget/localProxyWidget.cpp @@ -328,7 +328,10 @@ void Crit3DLocalProxyWidget::plot() delete label; } weightLabels.clear(); + int areaCode = NODATA; + if (interpolationSettings->getUseLocalDetrending()) + { if (detrended.isChecked()) { outInterpolationPoints.clear(); @@ -347,6 +350,43 @@ void Crit3DLocalProxyWidget::plot() outInterpolationPoints, checkSpatialQuality, errorStdStr); localSelection(outInterpolationPoints, subsetInterpolationPoints, x, y, *interpolationSettings); } + } + else if (interpolationSettings->getUseGlocalDetrending()) + { + areaCode = gis::getValueFromXY(*interpolationSettings->getMacroAreasMap(), x, y); + if (detrended.isChecked()) + { + outInterpolationPoints.clear(); + + checkAndPassDataToInterpolation(quality, myVar, meteoPoints, nrMeteoPoints, getCurrentTime(), SQinterpolationSettings, + interpolationSettings, meteoSettings, climateParam, + outInterpolationPoints, checkSpatialQuality, errorStdStr); + + for (int k = 0; k < outInterpolationPoints.size(); k++) + { + if (outInterpolationPoints[k].macroAreaCode == areaCode) + { + subsetInterpolationPoints.push_back(outInterpolationPoints[k]); + } + } + + detrending(subsetInterpolationPoints, interpolationSettings->getSelectedCombination(), interpolationSettings, climateParam, myVar, getCurrentTime()); + } + else + { + checkAndPassDataToInterpolation(quality, myVar, meteoPoints, nrMeteoPoints, getCurrentTime(), SQinterpolationSettings, + interpolationSettings, meteoSettings, climateParam, + outInterpolationPoints, checkSpatialQuality, errorStdStr); + + for (int k = 0; k < outInterpolationPoints.size(); k++) + { + if (outInterpolationPoints[k].macroAreaCode == areaCode) + { + subsetInterpolationPoints.push_back(outInterpolationPoints[k]); + } + } + } + } QList pointListPrimary, pointListSecondary, pointListSupplemental, pointListMarked; QMap< QString, QPointF > idPointMap1; QMap< QString, QPointF > idPointMap2; From 4df780726c4df45ca9d37d7b2b700e19859f3564 Mon Sep 17 00:00:00 2001 From: cate-to Date: Wed, 23 Oct 2024 09:44:47 +0200 Subject: [PATCH 03/34] changes in glocal detrending --- pragaProject/pragaProject.cpp | 2 +- project/project.cpp | 84 +++++++++++++++++++++-------------- project/project.h | 4 +- 3 files changed, 53 insertions(+), 37 deletions(-) diff --git a/pragaProject/pragaProject.cpp b/pragaProject/pragaProject.cpp index 9472a24fe..78ed1b01a 100644 --- a/pragaProject/pragaProject.cpp +++ b/pragaProject/pragaProject.cpp @@ -2639,7 +2639,7 @@ bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QL if (interpolationSettings.getUseGlocalDetrending()) { logInfoGUI("Loading macro areas map for glocal detrending..."); - if (! loadMacroAreaGlocalMap(false)) + if (loadMacroAreaGlocalMap() < 1) return false; } diff --git a/project/project.cpp b/project/project.cpp index 65f01c26b..4517beaeb 100644 --- a/project/project.cpp +++ b/project/project.cpp @@ -2142,28 +2142,16 @@ bool Project::loadTopographicDistanceMaps(bool onlyWithData, bool showInfo) return true; } -bool Project::loadMacroAreaGlocalMap(bool showInfo) +int Project::loadMacroAreaGlocalMap() { //TODO: add a check for the code values? - if (nrMeteoPoints == 0) - { - logError("Open a meteo points DB before."); - return false; - } - QString mapsFolder = defaultPath + PATH_GEO; if (! QDir(mapsFolder).exists()) { //errore? } - if (showInfo) - { - //QString infoStr = "Loading map of macro areas for glocal detrending..."; - //logInfo(infoStr); - } - std::string myError; gis::Crit3DRasterGrid* macroAreasGrid = new gis::Crit3DRasterGrid(); std::string fileName = mapsFolder.toStdString() + glocalMapName.toStdString(); @@ -2184,40 +2172,68 @@ bool Project::loadMacroAreaGlocalMap(bool showInfo) } - //associa il codice a ogni stazione + //associa il codice a ogni stazione e trova il numero di aree tot + int numAreas = 1; for (int i=0; i < nrMeteoPoints; i++) { int myValue = gis::getValueFromXY(*macroAreasGrid, meteoPoints[i].point.utm.x, meteoPoints[i].point.utm.y); + if (myValue > numAreas) + numAreas = myValue; + if (myValue != (*macroAreasGrid).header->flag) meteoPoints[i].macroAreaCode = myValue; } - //if (showInfo) closeLogInfo(); - - return true; + return numAreas+1; } -bool Project::groupCellsInArea(std::vector &areaPoints, unsigned int index) +bool Project::groupCellsInArea(std::vector &areaPoints, unsigned int index, bool isGrid) { int zoneNr; double myX, myY; gis::Crit3DRasterGrid* macroAreas = interpolationSettings.getMacroAreasMap(); - int nrCols = DEM.header->nrCols; + int nrCols; + areaPoints.clear(); + + if (!isGrid) + { + nrCols = DEM.header->nrCols; - //DEM. GRID ANCORA DA SCRIVERE - for (int i = 0; i < DEM.header->nrRows; i++) + for (int i = 0; i < DEM.header->nrRows; i++) + { + for (int j = 0; j < DEM.header->nrCols; j++) + { + if (!isEqual(DEM.getValueFromRowCol(i,j),DEM.header->flag)) + { + gis::getUtmXYFromRowCol(DEM.header, i, j, &myX, &myY); + + zoneNr = int(macroAreas->getValueFromXY(myX, myY)); + + if (zoneNr == index) + areaPoints.push_back(i*nrCols+j); + } + } + } + } + else { - for (int j = 0; j < DEM.header->nrCols; j++) + nrCols = meteoGridDbHandler->meteoGrid()->gridStructure().header().nrCols; + + for (int i = 0; i < meteoGridDbHandler->meteoGrid()->gridStructure().header().nrRows; i++) { - if (!isEqual(DEM.getValueFromRowCol(i,j),DEM.header->flag)) + for (int j = 0; j < meteoGridDbHandler->meteoGrid()->gridStructure().header().nrCols; j++) { - gis::getUtmXYFromRowCol(DEM.header, i, j, &myX, &myY); + if (meteoGridDbHandler->meteoGrid()->meteoPoints()[i][j]->active) + { + myX = meteoGridDbHandler->meteoGrid()->meteoPoints()[i][j]->point.utm.x; + myY = meteoGridDbHandler->meteoGrid()->meteoPoints()[i][j]->point.utm.y; - zoneNr = int(macroAreas->getValueFromXY(myX, myY)); + zoneNr = int(macroAreas->getValueFromXY(myX, myY)); - if (zoneNr == index) - areaPoints.push_back(i*nrCols+j); + if (zoneNr == index) + areaPoints.push_back(i*nrCols+j); + } } } } @@ -2638,7 +2654,7 @@ bool Project::interpolationDemGlocalDetrending(meteoVariable myVar, const Crit3D { //save all the cells in that specific area in a vector (areaCells) areaCells.clear(); - groupCellsInArea(areaCells, areaIndex); + groupCellsInArea(areaCells, areaIndex, false); //mettere qui calcolo buffer, riutilizzo funzione distanza da area //take the parameters+combination for that area std::vector &)>> fittingFunction; @@ -2844,7 +2860,7 @@ bool Project::interpolationDemMain(meteoVariable myVar, const Crit3DTime& myTime if (interpolationSettings.getUseGlocalDetrending()) { - if (!loadMacroAreaGlocalMap(true)) return false; + if (loadMacroAreaGlocalMap() < 1) return false; } // dynamic lapserate if (getUseDetrendingVar(myVar) && interpolationSettings.getUseLocalDetrending()) @@ -2900,11 +2916,11 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) if (interpolationSettings.getUseMultipleDetrending()) interpolationSettings.clearFitting(); - unsigned numAreas; + unsigned numAreas = 1; if (interpolationSettings.getUseGlocalDetrending()) { - if (!loadMacroAreaGlocalMap(true)) return false; - numAreas = interpolationSettings.getAreaParameters().size(); + numAreas = loadMacroAreaGlocalMap(); + if (numAreas < 1) return false; } // check quality and pass data to interpolation @@ -3044,7 +3060,7 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) for (unsigned areaIndex = 0; areaIndex < numAreas; areaIndex++) { areaCells.clear(); - groupCellsInArea(areaCells, areaIndex); + groupCellsInArea(areaCells, areaIndex, true); std::vector &)> > fittingFunction; interpolationSettings.setFittingParameters(allAreaFittingParameters[areaIndex]); interpolationSettings.setCurrentCombination(allAreaCombinations[areaIndex]); @@ -4081,7 +4097,7 @@ void Project::showLocalProxyGraph(gis::Crit3DGeoPoint myPoint) } if (interpolationSettings.getUseGlocalDetrending()) - loadMacroAreaGlocalMap(false); + loadMacroAreaGlocalMap(); localProxyWidget = new Crit3DLocalProxyWidget(myUtm.x, myUtm.y, myZDEM, myZGrid, this->gisSettings, &interpolationSettings, meteoPoints, nrMeteoPoints, currentVariable, currentFrequency, currentDate, currentHour, quality, &qualityInterpolationSettings, meteoSettings, &climateParameters, checkSpatialQuality); return; diff --git a/project/project.h b/project/project.h index 7b217b0c9..2a9790021 100644 --- a/project/project.h +++ b/project/project.h @@ -260,8 +260,8 @@ bool writeTopographicDistanceMap(int pointIndex, const gis::Crit3DRasterGrid& demMap, QString pathTd); bool loadTopographicDistanceMaps(bool onlyWithData, bool showInfo); void passInterpolatedTemperatureToHumidityPoints(Crit3DTime myTime, Crit3DMeteoSettings *meteoSettings); - bool loadMacroAreaGlocalMap(bool showInfo); - bool groupCellsInArea(std::vector &areaPoints, unsigned index); + int loadMacroAreaGlocalMap(); + bool groupCellsInArea(std::vector &areaPoints, unsigned index, bool isGrid); bool checkInterpolation(meteoVariable myVar); bool checkInterpolationGrid(meteoVariable myVar); From 56c1aed168bcca5f5cc3a6ea027fed23a4d17388 Mon Sep 17 00:00:00 2001 From: avolta Date: Wed, 23 Oct 2024 16:24:33 +0200 Subject: [PATCH 04/34] aggiunte funzioni di marquardt senza pesi --- mathFunctions/furtherMathFunctions.cpp | 504 ++++++++++++++++++++++++- mathFunctions/furtherMathFunctions.h | 30 ++ 2 files changed, 529 insertions(+), 5 deletions(-) diff --git a/mathFunctions/furtherMathFunctions.cpp b/mathFunctions/furtherMathFunctions.cpp index e004c151e..bb1456450 100644 --- a/mathFunctions/furtherMathFunctions.cpp +++ b/mathFunctions/furtherMathFunctions.cpp @@ -1413,6 +1413,23 @@ namespace interpolation return sqrt(weighted_MSE); } + double computeRMSE(const std::vector& observed, const std::vector& predicted) + { + // This function computes the root mean squared error + double sum_residuals = 0.0; + + // Calculate the sum of the squared errors + for (int i = 0; i < int(observed.size()); i++) + { + double residual = (observed[i] - predicted[i]) * (observed[i] - predicted[i]); + sum_residuals += residual; + } + + // Calculate MSE + double MSE = sum_residuals / (observed.size()); + + return sqrt(MSE); + } int bestFittingMarquardt_nDimension(double (*func)(std::vector&)>>&, std::vector& , std::vector >&), std::vector&)>>& myFunc, @@ -1602,6 +1619,51 @@ namespace interpolation } + int bestFittingMarquardt_nDimension_clean(double (*func)(std::vector&)>>&, std::vector& , std::vector >&), + std::vector&)>>& myFunc, + std::vector >& parametersMin, std::vector >& parametersMax, + std::vector >& parameters, std::vector >& parametersDelta, + int maxIterationsNr, double myEpsilon, + std::vector >& x ,std::vector& y) + { + int i,j; + //int nrData = int(y.size()); + std::vector nrParameters(parameters.size()); + int nrParametersTotal = 0; + for (i=0; i> bestParameters(parameters.size()); + std::vector > correspondenceTag(2,std::vector(nrParametersTotal)); + int counterTag = 0; + for (i=0; i ySim(nrData); + + fittingMarquardt_nDimension_noSquares(func,myFunc,parametersMin, parametersMax, + parameters, parametersDelta,correspondenceTag, maxIterationsNr, + myEpsilon, x, y); + + //for (i=0;i&)>>&, std::vector& , std::vector >&), std::vector &)> >& myFunc, @@ -1723,11 +1785,6 @@ namespace interpolation paramChange[i].resize(nrParameters[i]); newParameters[i].resize(nrParameters[i]); lambda[i].resize(nrParameters[i],0.01); - /*for (j=0; j myEpsilon && iterationNr <= maxIterationsNr); return (fabs(diffSSE) <= myEpsilon); } + bool fittingMarquardt_nDimension_noSquares(double (*func)(std::vector&)>>&, std::vector& , std::vector >&), + std::vector &)> >& myFunc, + std::vector > ¶metersMin, std::vector > ¶metersMax, + std::vector > ¶meters, std::vector > ¶metersDelta, std::vector > &correspondenceParametersTag, + int maxIterationsNr, double myEpsilon, + std::vector >& x, std::vector& y) + { + int i; + int nrPredictors = 0; + for (int k = 0; k < parameters.size(); k++) + if (parameters[k].size() == 2) nrPredictors++; + int nrData = int(y.size()); + double mySSE, diffSSE, newSSE; + static double VFACTOR = 10; + std::vector nrParameters(nrPredictors); + std::vector > paramChange(nrPredictors); + std::vector > newParameters(nrPredictors); + std::vector > lambda(nrPredictors); + + for (i=0; i nrParameters(nrPredictors); + for (i=0; i g(nrParametersTotal,0); + //std::vector z(nrParametersTotal); + std::vector firstEst(nrData); + std::vector> a(nrParametersTotal, std::vector(nrParametersTotal,0)); + std::vector> P(nrParametersTotal, std::vector(nrData)); + + // matrix P corresponds to the Jacobian + // first set of estimates + for (i = 0; i < nrData; i++) + { + firstEst[i] = func(myFunc,x[i], parameters); + } + + // change parameters and compute derivatives + int counterDim = 0; + //double newEst; + for (i = 0; i < nrPredictors; i++) + { + for (k=0;k= 0; i--) + { + top = g[i]; + for (k = i + 1; k < nrParametersTotal; k++) + { + top -= a[i][k] * paramChange[correspondenceParametersTag[0][k]][correspondenceParametersTag[1][k]]; + } + paramChange[correspondenceParametersTag[0][i]][correspondenceParametersTag[1][i]] = top / a[i][i]; + } + // change parameters + for (i = 0; i < nrPredictors; i++) + { + for (j=0; j parametersMax[i][j]) && (lambda[i][j] < 1000)) + { + newParameters[i][j] = parametersMax[i][j]; + if (lambda[i][j] < 1000) + lambda[i][j] *= VFACTOR; + } + if (newParameters[i][j] < parametersMin[i][j]) + { + newParameters[i][j] = parametersMin[i][j]; + if (lambda[i][j] < 1000) + lambda[i][j] *= VFACTOR; + } + } + } + newSSE = normGeneric_nDimension(func, myFunc, newParameters, x, y); + + if (newSSE == NODATA) + return false; + + diffSSE = mySSE - newSSE ; + + if (diffSSE > 0) + { + mySSE = newSSE; + for (i=0; i myEpsilon && iterationNr <= maxIterationsNr); + return (fabs(diffSSE) <= myEpsilon); + } void leastSquares_nDimension(double (*func)(std::vector &)> > &, std::vector &, std::vector >&), std::vector &)> > myFunc, @@ -2209,7 +2461,21 @@ namespace interpolation } return norm; } + double normGeneric_nDimension(double (*func)(std::vector&)>>&, std::vector&, std::vector >&), + std::vector&)>> myFunc, + std::vector > ¶meters,std::vector >& x, + std::vector& y) + { + double error; + double norm = 0; + for (int i = 0; i < int(y.size()); i++) + { + error = y[i] - func(myFunc,x[i], parameters); + norm += error * error; + } + return norm; + } double bestFittingMarquardt_nDimension_singleFunction(double (*func)(double, std::vector&), int nrTrials, int nrMinima, std::vector & parametersMin, std::vector & parametersMax, @@ -2280,6 +2546,77 @@ namespace interpolation return bestR2; } + double bestFittingMarquardt_nDimension_singleFunction(double (*func)(double, std::vector&), + int nrTrials, int nrMinima, + std::vector & parametersMin, std::vector & parametersMax, + std::vector & parameters, std::vector & parametersDelta, + std::vector & stepSize, int numSteps, + int maxIterationsNr, double myEpsilon, double deltaR2, + std::vector & x ,std::vector& y, + std::vector> firstGuessCombinations) + { + int i,j; + int nrData = int(y.size()); + int nrParameters = int(parameters.size()); + + std::vector bestParameters(nrParameters); + + double bestR2 = NODATA; + double bestRMSE = NODATA; + int RMSEindex = NODATA; + double R2, RMSE; + std::vector R2Previous(nrMinima,NODATA); + std::vector ySim(nrData); + + //grigliato + for (int k = 0; k < firstGuessCombinations.size(); k++) + { + parameters = firstGuessCombinations[k]; + fittingMarquardt_nDimension_noSquares_singleFunction(func,parametersMin, + parametersMax,parameters, + parametersDelta,maxIterationsNr, + myEpsilon,x,y); + + for (i=0;i (bestR2 + deltaR2)) + { + for (j=0; j&), std::vector ¶metersMin, std::vector ¶metersMax, std::vector ¶meters, std::vector ¶metersDelta, @@ -2446,6 +2783,163 @@ namespace interpolation return (fabs(diffSSE) <= myEpsilon); } + bool fittingMarquardt_nDimension_noSquares_singleFunction(double (*func) (double, std::vector&), + std::vector ¶metersMin, std::vector ¶metersMax, + std::vector ¶meters, std::vector ¶metersDelta, + int maxIterationsNr, double myEpsilon, + std::vector & x, std::vector& y) + { + int i,j,k; + double mySSE, diffSSE, newSSE; + static double VFACTOR = 10; + int nrParameters = int(parameters.size()); + std::vector paramChange(nrParameters,0); + std::vector newParameters(nrParameters,0); + std::vector lambda (nrParameters,0.01); + int nrData = int(y.size()); + + double error = 0; + mySSE = 0; + for (i = 0; i < nrData; i++) + { + error = y[i] - func(x[i], parameters); + mySSE += error * error; + } + int iterationNr = 0; + do + { + double pivot, mult, top; + std::vector g(nrParameters,0); + std::vector firstEst(nrData); + std::vector> a(nrParameters, std::vector(nrParameters,0)); + std::vector> P(nrParameters, std::vector(nrData)); + //std::vector> weightsP(nrParameters, std::vector(nrData)); + + // matrix P corresponds to the Jacobian + // first set of estimates + for (i = 0; i < nrData; i++) + { + firstEst[i] = func(x[i], parameters); + } + + // change parameters and compute derivatives + for (k=0;k= 0; i--) + { + top = g[i]; + for (k = i + 1; k < nrParameters; k++) + { + top -= a[i][k] * paramChange[k]; + } + paramChange[i] = top / std::max(a[i][i], EPSILON); + } + + // change parameters + for (j=0; j parametersMax[j]) && (lambda[j] < 1000)) + { + newParameters[j] = parametersMax[j]; + if (lambda[j] < 1000) + lambda[j] *= VFACTOR; + } + if (newParameters[j] < parametersMin[j]) + { + newParameters[j] = parametersMin[j]; + if (lambda[j] < 1000) + lambda[j] *= VFACTOR; + } + } + + newSSE = 0; + for (i = 0; i < nrData; i++) + { + error = y[i] - func(x[i], newParameters); + newSSE += error * error; + } + + if (newSSE == NODATA) + return false; + + diffSSE = mySSE - newSSE ; + + if (diffSSE > 0) + { + mySSE = newSSE; + for (j=0; j myEpsilon && iterationNr <= maxIterationsNr); + return (fabs(diffSSE) <= myEpsilon); + } + } diff --git a/mathFunctions/furtherMathFunctions.h b/mathFunctions/furtherMathFunctions.h index 3f4dbb30f..9092fa2e1 100644 --- a/mathFunctions/furtherMathFunctions.h +++ b/mathFunctions/furtherMathFunctions.h @@ -102,6 +102,7 @@ enum estimatedFunction {FUNCTION_CODE_SPHERICAL, FUNCTION_CODE_LINEAR, FUNCTION_ void tridiagonalThomasAlgorithm (int n, double *subDiagonal, double *mainDiagonal, double *superDiagonal, double *constantTerm, double* output); // not working to be checked double computeR2(const std::vector& obs, const std::vector& sim); + double computeWeightedRMSE(const std::vector& observed, const std::vector& predicted); double computeWeighted_R2(const std::vector& observed, const std::vector& predicted, const std::vector& weights); double computeWeighted_R2_secondFormulation(const std::vector& observed, const std::vector& predicted, const std::vector& weights); double computeWeighted_R2_generalized_independentObservedData(const std::vector& observed, const std::vector& predicted, const std::vector& weights); @@ -125,6 +126,12 @@ enum estimatedFunction {FUNCTION_CODE_SPHERICAL, FUNCTION_CODE_LINEAR, FUNCTION_ int maxIterationsNr, double myEpsilon, std::vector >& x ,std::vector& y, std::vector& weights); + int bestFittingMarquardt_nDimension_clean(double (*func)(std::vector&)>>&, std::vector& , std::vector >&), + std::vector&)>>& myFunc, + std::vector >& parametersMin, std::vector >& parametersMax, + std::vector >& parameters, std::vector >& parametersDelta, + int maxIterationsNr, double myEpsilon, + std::vector >& x ,std::vector& y); bool fittingMarquardt_nDimension(double (*func)(std::vector &)> > &, std::vector &, std::vector >&), std::vector &)> > &myFunc, @@ -140,12 +147,24 @@ enum estimatedFunction {FUNCTION_CODE_SPHERICAL, FUNCTION_CODE_LINEAR, FUNCTION_ int maxIterationsNr, double myEpsilon, std::vector >& x, std::vector& y, std::vector& weights); + bool fittingMarquardt_nDimension_noSquares(double (*func)(std::vector&)>>&, std::vector& , std::vector >&), + std::vector &)> >& myFunc, + std::vector > ¶metersMin, std::vector > ¶metersMax, + std::vector > ¶meters, std::vector > ¶metersDelta, std::vector > &correspondenceParametersTag, + int maxIterationsNr, double myEpsilon, + std::vector >& x, std::vector& y); + bool fittingMarquardt_nDimension_noSquares_singleFunction(double (*func) (double, std::vector&), std::vector ¶metersMin, std::vector ¶metersMax, std::vector ¶meters, std::vector ¶metersDelta, int maxIterationsNr, double myEpsilon, std::vector & x, std::vector& y, std::vector& weights); + bool fittingMarquardt_nDimension_noSquares_singleFunction(double (*func) (double, std::vector&), + std::vector ¶metersMin, std::vector ¶metersMax, + std::vector ¶meters, std::vector ¶metersDelta, + int maxIterationsNr, double myEpsilon, + std::vector & x, std::vector& y); double bestFittingMarquardt_nDimension_singleFunction(double (*func)(double, std::vector&), int nrTrials, int nrMinima, std::vector & parametersMin, std::vector & parametersMax, @@ -154,10 +173,21 @@ enum estimatedFunction {FUNCTION_CODE_SPHERICAL, FUNCTION_CODE_LINEAR, FUNCTION_ int maxIterationsNr, double myEpsilon, double deltaR2, std::vector & x , std::vector& y, std::vector& weights, std::vector > firstGuessCombinations); + double bestFittingMarquardt_nDimension_singleFunction(double (*func)(double, std::vector&), + int nrTrials, int nrMinima, + std::vector & parametersMin, std::vector & parametersMax, + std::vector & parameters, std::vector & parametersDelta, + std::vector &stepSize, int numSteps, + int maxIterationsNr, double myEpsilon, double deltaR2, + std::vector & x , std::vector& y, + std::vector > firstGuessCombinations); double normGeneric_nDimension(double (*func)(std::vector &)>> &, std::vector &, std::vector >&), std::vector &)> > myFunc, std::vector > ¶meters, std::vector >& x, std::vector& y, std::vector& weights); + double normGeneric_nDimension(double (*func)(std::vector &)>> &, std::vector &, std::vector >&), + std::vector &)> > myFunc, + std::vector > ¶meters, std::vector >& x, std::vector& y); void leastSquares_nDimension(double (*func)(std::vector&)>>&, std::vector& , std::vector >&), std::vector &)> > myFunc, From 11a50f61d57b0e99b6f7844552b63ccbbbf5215a Mon Sep 17 00:00:00 2001 From: Gabriele Antolini Date: Wed, 23 Oct 2024 17:31:32 +0200 Subject: [PATCH 05/34] glocal weight maps --- mathFunctions/commonConstants.h | 1 + project/interpolationCmd.cpp | 1 + project/project.cpp | 92 +++++++++++++++++++++++++++++++++ project/project.h | 1 + 4 files changed, 95 insertions(+) diff --git a/mathFunctions/commonConstants.h b/mathFunctions/commonConstants.h index 1342b4961..f507f7b06 100644 --- a/mathFunctions/commonConstants.h +++ b/mathFunctions/commonConstants.h @@ -51,6 +51,7 @@ #define PATH_LOG "LOG/" #define PATH_OUTPUT "OUTPUT/" #define PATH_TD "TD/" + #define PATH_GLOCAL "GLOCAL/" #define PATH_STATES "STATES/" #define PATH_NETCDF "NETCDF/" diff --git a/project/interpolationCmd.cpp b/project/interpolationCmd.cpp index 045b0adcb..cd452a287 100644 --- a/project/interpolationCmd.cpp +++ b/project/interpolationCmd.cpp @@ -380,3 +380,4 @@ bool topographicIndex(const gis::Crit3DRasterGrid& DEM, std::vector wind return true; } + diff --git a/project/project.cpp b/project/project.cpp index 4517beaeb..6f89205e1 100644 --- a/project/project.cpp +++ b/project/project.cpp @@ -2142,6 +2142,98 @@ bool Project::loadTopographicDistanceMaps(bool onlyWithData, bool showInfo) return true; } +bool Project::glocalWeightsMaps(float windowWidth) +{ + //controlla se c'è mappa aree + + gis::Crit3DRasterGrid* zoneGrid = interpolationSettings.getMacroAreasMap(); + if (! zoneGrid->isLoaded) + { + errorString = "Load a zone grid before"; + return false; + } + + QString mapsFolder = projectPath + PATH_GLOCAL; + if (! QDir(mapsFolder).exists()) + { + QDir().mkdir(mapsFolder); + } + + std::vector zoneCount((unsigned int)(zoneGrid->maximum)); + + int i, j; + + //create raster for each zone + std::vector zoneWeightMaps((unsigned int)(zoneGrid->maximum)); + for (i=0; i < zoneGrid->maximum; i++) + zoneWeightMaps[i].initializeGrid(*zoneGrid); + + int cellNr = round(windowWidth / zoneGrid->header->cellSize); + float cellDelta; + + int r1, r2, c1, c2, windowRow, windowCol, cellCount; + + for (int zoneRow = 0; zoneRow < zoneGrid->header->nrRows; zoneRow++) + { + for (int zoneCol = 0; zoneCol < zoneGrid->header->nrCols; zoneCol++) + { + if (! isEqual(zoneGrid->value[zoneRow][zoneCol], zoneGrid->header->flag)) + { + cellCount = 0; + + for (j=0; j < zoneGrid->maximum; j++) + zoneCount[j] = 0; + + r1 = zoneRow - cellNr; + r2 = zoneRow + cellNr; + c1 = zoneCol - cellNr; + c2 = zoneCol + cellNr; + + for (windowRow = r1; windowRow <= r2; windowRow++) + { + for (windowCol = c1; windowCol <= c2; windowCol++) + { + if (! gis::isOutOfGridRowCol(windowRow, windowCol, *zoneGrid)) + { + i = zoneGrid->value[windowRow][windowCol]; + + if (! isEqual(i, zoneGrid->header->flag)) + { + cellDelta = gis::computeDistance(windowRow, windowCol, zoneRow, zoneCol); + + if (cellDelta <= cellNr) + { + cellCount++; + zoneCount[i-1]++; + } + } + } + } + } + + for (i=0; i < zoneGrid->maximum; i++) + if (cellCount > 0) + zoneWeightMaps[i].value[zoneRow][zoneCol] = float(zoneCount[i]) / float(cellCount); + + } + } + } + + QString fileName; + std::string myError; + for (int i=0; i < zoneGrid->maximum; i++) + { + fileName = mapsFolder + "glocalWeight_" + QString::number(i+1); + if (! gis::writeEsriGrid(fileName.toStdString(), &zoneWeightMaps[i], myError)) + { + errorString = QString::fromStdString(myError); + return false; + } + } + + return true; +} + int Project::loadMacroAreaGlocalMap() { //TODO: add a check for the code values? diff --git a/project/project.h b/project/project.h index 2a9790021..d531b02a0 100644 --- a/project/project.h +++ b/project/project.h @@ -262,6 +262,7 @@ void passInterpolatedTemperatureToHumidityPoints(Crit3DTime myTime, Crit3DMeteoSettings *meteoSettings); int loadMacroAreaGlocalMap(); bool groupCellsInArea(std::vector &areaPoints, unsigned index, bool isGrid); + bool glocalWeightsMaps(float windowWidth); bool checkInterpolation(meteoVariable myVar); bool checkInterpolationGrid(meteoVariable myVar); From b6b39dbe5989b504bd693dbae4f2145661618b52 Mon Sep 17 00:00:00 2001 From: cate-to Date: Thu, 24 Oct 2024 15:33:34 +0200 Subject: [PATCH 06/34] changes in glocal detrending --- gis/gis.cpp | 31 +++++ gis/gis.h | 2 + interpolation/interpolation.cpp | 59 ++++---- interpolation/interpolationPoint.cpp | 1 - interpolation/interpolationPoint.h | 1 - interpolation/interpolationSettings.cpp | 77 +++++++++++ interpolation/interpolationSettings.h | 29 ++++ interpolation/spatialControl.cpp | 1 - meteo/meteoPoint.cpp | 1 - meteo/meteoPoint.h | 1 - pragaProject/pragaProject.cpp | 2 +- project/project.cpp | 171 +++++++++++++++++------- project/project.h | 5 +- proxyWidget/localProxyWidget.cpp | 25 ++-- 14 files changed, 308 insertions(+), 98 deletions(-) diff --git a/gis/gis.cpp b/gis/gis.cpp index e877c8bed..f447beaf2 100644 --- a/gis/gis.cpp +++ b/gis/gis.cpp @@ -593,6 +593,37 @@ namespace gis return sqrtf((dx * dx) + (dy * dy)); } + std::vector computeEuclideanDistanceStation2Area(std::vector>& cells,std::vector>& stations) + { + // è possibile sapere in quale cella (row,col) si trova la stazione? + std::vector distance(stations.size()); + for (int i=0;i computeMetropolisDistanceStation2Area(std::vector>& cells,std::vector>& stations) + { + // è possibile sapere in quale cella (row,col) si trova la stazione? + std::vector distance(stations.size()); + for (int i=0;i computeEuclideanDistanceStation2Area(std::vector>& cells,std::vector>& stations); + std::vector computeMetropolisDistanceStation2Area(std::vector>& cells,std::vector>& stations); bool updateMinMaxRasterGrid(Crit3DRasterGrid *rasterGrid); void convertFlagToNodata(Crit3DRasterGrid& myGrid); bool updateColorScale(Crit3DRasterGrid* rasterGrid, int row0, int col0, int row1, int col1); diff --git a/interpolation/interpolation.cpp b/interpolation/interpolation.cpp index f5e7522af..a3139d781 100644 --- a/interpolation/interpolation.cpp +++ b/interpolation/interpolation.cpp @@ -1122,15 +1122,9 @@ void localSelection(vector &inputPoints, vector < r1 += stepRadius; } - double temp = 0; - if (maxDistance != 0) for (i=0; i< selectedPoints.size(); i++) { - /*temp = (MAXVALUE(1 - selectedPoints[i].distance / maxDistance, EPSILON))*1000; - temp = round(temp); - selectedPoints[i].regressionWeight = temp/1000;*/ - selectedPoints[i].regressionWeight = MAXVALUE(1 - selectedPoints[i].distance / maxDistance, EPSILON); //selectedPoints[i].regressionWeight = MAXVALUE(std::exp(-selectedPoints[i].distance*selectedPoints[i].distance/((0.5*maxDistance)*(0.5*maxDistance))),EPSILON); @@ -1518,7 +1512,7 @@ void calculateFirstGuessCombinations(Crit3DProxy* myProxy) std::vector tempParam = myProxy->getFittingParametersRange(); std::vector firstGuessPosition = myProxy->getFittingFirstGuess(); std::vector tempFirstGuess; - int numSteps = 13; + int numSteps = 15; std::vector stepSize; int nrParam = int(tempParam.size()/2); @@ -1792,7 +1786,7 @@ bool multipleDetrendingElevationFitting(int elevationPos, std::vector parametersDelta; std::vector parameters; std::vector stepSize; - int numSteps = 10; + int numSteps = 15; std::function&)> myFunc; @@ -1801,7 +1795,7 @@ bool multipleDetrendingElevationFitting(int elevationPos, std::vector > firstGuessCombinations = mySettings->getProxy(elevationPos)->getFirstGuessCombinations(); // multiple non linear fitting @@ -2068,43 +2062,38 @@ bool glocalDetrendingFitting(std::vector &myPoint //matrice avrà, per ogni area, i parametri di fitting + la combination + il valore delle stazioni detrendate secondo il fitting corrispondente std::vector subsetPoints; int i = 0; - int j = 0; - int areasNr = int(myPoints.front().macroAreaCode); - int start = areasNr; + std::vector macroAreas = mySettings->getMacroAreas(); + int areasNr = macroAreas.size(); - //look for highest and lowest "macroAreaCode", that's the areas we'll have to iterate the fitting for - for (i = 0; i < myPoints.size(); i++) + int elevationPos = NODATA; + for (unsigned int pos=0; pos < mySettings->getCurrentCombination().getProxySize(); pos++) { - if (myPoints[i].macroAreaCode > areasNr) - areasNr = int(myPoints[i].macroAreaCode); - if (myPoints[i].macroAreaCode < start && !(myPoints[i].macroAreaCode == NODATA)) - start = int(myPoints[i].macroAreaCode); + if (getProxyPragaName(mySettings->getProxy(pos)->getName()) == proxyHeight) + elevationPos = pos; } std::vector>> allAreaParameters; std::vector allAreaCombinations; //now create the subset of points, including only the right points, and then call multipleDetrending - for (i = start-1; i < areasNr+1; i++) //ATTENZIONE INDICI + for (i = 0; i < areasNr; i++) //ATTENZIONE INDICI { + Crit3DMacroArea myArea = macroAreas[i]; - for (j = 0; j < myPoints.size(); j++) + std::vector temp = myArea.getMeteoPoints(); + if (! temp.empty()) { - if (myPoints[j].macroAreaCode == i) - subsetPoints.push_back(myPoints[j]); - } + std::vector subsetPoints; + for (int l = 0; l < temp.size(); l++) + { + + for (int k = 0; k < myPoints.size(); k++) + if (myPoints[k].index == temp[l]) + subsetPoints.push_back(myPoints[k]); + } - if (! subsetPoints.empty()) - { mySettings->setCurrentCombination(mySettings->getSelectedCombination()); mySettings->clearFitting(); - int elevationPos = NODATA; - for (unsigned int pos=0; pos < mySettings->getCurrentCombination().getProxySize(); pos++) - { - if (getProxyPragaName(mySettings->getProxy(pos)->getName()) == proxyHeight) - elevationPos = pos; //SPOSTA - } - if (elevationPos != NODATA && mySettings->getCurrentCombination().isProxyActive(elevationPos)) { if (!multipleDetrendingElevationFitting(elevationPos, subsetPoints, mySettings, myVar, errorStr)) return false; @@ -2117,7 +2106,7 @@ bool glocalDetrendingFitting(std::vector &myPoint detrendingOtherProxies(elevationPos, subsetPoints, mySettings); - //salvataggio parametri e valori staz(?) parlane con gabri pk dovresti salvare tutti i crit3dintpoint(?) + // allAreaParameters.push_back(mySettings->getFittingParameters()); allAreaCombinations.push_back(mySettings->getCurrentCombination()); } @@ -2245,8 +2234,8 @@ bool preInterpolation(std::vector &myPoints, Crit { if (!mySettings->getUseLocalDetrending()) setMultipleDetrendingHeightTemperatureRange(mySettings); - - if (mySettings->getUseGlocalDetrending()) + + if (mySettings->getUseGlocalDetrending()) { glocalDetrendingFitting(myPoints, mySettings, myVar, errorStr); } diff --git a/interpolation/interpolationPoint.cpp b/interpolation/interpolationPoint.cpp index 5e9aaff39..c456f72fd 100644 --- a/interpolation/interpolationPoint.cpp +++ b/interpolation/interpolationPoint.cpp @@ -39,7 +39,6 @@ Crit3DInterpolationDataPoint::Crit3DInterpolationDataPoint() regressionWeight = NODATA; lapseRateCode = primary; - macroAreaCode = NODATA; topographicDistance = nullptr; diff --git a/interpolation/interpolationPoint.h b/interpolation/interpolationPoint.h index cd641a803..82caaef55 100644 --- a/interpolation/interpolationPoint.h +++ b/interpolation/interpolationPoint.h @@ -25,7 +25,6 @@ lapseRateCodeType lapseRateCode; gis::Crit3DRasterGrid* topographicDistance; std::vector proxyValues; - int macroAreaCode; Crit3DInterpolationDataPoint(); diff --git a/interpolation/interpolationSettings.cpp b/interpolation/interpolationSettings.cpp index 51ffb12c5..e6c5860e5 100644 --- a/interpolation/interpolationSettings.cpp +++ b/interpolation/interpolationSettings.cpp @@ -461,6 +461,9 @@ void Crit3DInterpolationSettings::initialize() fittingFunction.clear(); fittingParameters.clear(); + macroAreas.clear(); + areaCombination.clear(); + areaParameters.clear(); isKrigingReady = false; precipitationAllZero = false; @@ -538,6 +541,20 @@ void Crit3DInterpolationSettings::setAreaCombination(std::vector Crit3DInterpolationSettings::getMacroAreas() +{ + return macroAreas; +} +void Crit3DInterpolationSettings::addMacroArea(Crit3DMacroArea myArea) +{ + macroAreas.push_back(myArea); +} + +void Crit3DInterpolationSettings::setMacroAreas(std::vector myAreas) +{ + macroAreas = myAreas; +} + TInterpolationMethod Crit3DInterpolationSettings::getInterpolationMethod() { return interpolationMethod;} @@ -1003,3 +1020,63 @@ bool Crit3DInterpolationSettings::getCombination(int combinationInteger, Crit3DP return true; } +Crit3DMacroArea::Crit3DMacroArea() +{ + meteoPoints.clear(); + areaParameters.clear(); + areaCombination.clear(); + areaCells.clear(); +} + +void Crit3DMacroArea::setMeteoPoints (std::vector myMeteoPoints) +{ + meteoPoints = myMeteoPoints; + return; +} + +std::vector Crit3DMacroArea::getMeteoPoints() +{ + return meteoPoints; +} + +void Crit3DMacroArea::setAreaCells (std::vector myCells) +{ + areaCells = myCells; + return; +} + +std::vector Crit3DMacroArea::getAreaCells() +{ + return areaCells; +} + +void Crit3DMacroArea::setParameters (std::vector> myParameters) +{ + areaParameters = myParameters; + return; +} + +std::vector> Crit3DMacroArea::getParameters() +{ + return areaParameters; +} + +void Crit3DMacroArea::setCombination (Crit3DProxyCombination myCombination) +{ + areaCombination = myCombination; + return; +} + +void Crit3DMacroArea::clear() +{ + areaCells.clear(); + areaParameters.clear(); + areaCombination.clear(); + meteoPoints.clear(); +} + +Crit3DProxyCombination Crit3DMacroArea::getCombination() +{ + return areaCombination; +} + diff --git a/interpolation/interpolationSettings.h b/interpolation/interpolationSettings.h index 42357ab54..cc665edd8 100644 --- a/interpolation/interpolationSettings.h +++ b/interpolation/interpolationSettings.h @@ -132,6 +132,29 @@ void setUseThermalInversion(bool value) { _useThermalInversion = value; } }; + class Crit3DMacroArea + { + private: + Crit3DProxyCombination areaCombination; + std::vector> areaParameters; + std::vector meteoPoints; + std::vector areaCells; + + public: + Crit3DMacroArea(); + + void setMeteoPoints (std::vector myMeteoPoints); + std::vector getMeteoPoints(); + void setParameters (std::vector> myParameters); + std::vector> getParameters(); + void setCombination (Crit3DProxyCombination myCombination); + Crit3DProxyCombination getCombination(); + void setAreaCells (std::vector myCells); + std::vector getAreaCells(); + void clear(); + + }; + class Crit3DInterpolationSettings { @@ -180,6 +203,8 @@ std::vector>> areaParameters; std::vector areaCombination; + std::vector macroAreas; + public: Crit3DInterpolationSettings(); @@ -215,6 +240,10 @@ void setAreaParameters(std::vector>> myAreaParameters); std::vector getAreaCombination(); void setAreaCombination(std::vector myAreaCombination); + std::vector getMacroAreas(); + void addMacroArea(Crit3DMacroArea myArea); + void setMacroAreas(std::vector myAreas); + void setUseDoNotRetrend(bool myValue); bool getUseDoNotRetrend(); diff --git a/interpolation/spatialControl.cpp b/interpolation/spatialControl.cpp index 5f9a9d532..c450e51da 100644 --- a/interpolation/spatialControl.cpp +++ b/interpolation/spatialControl.cpp @@ -429,7 +429,6 @@ bool passDataToInterpolation(Crit3DMeteoPoint* meteoPoints, int nrMeteoPoints, myPoint.topographicDistance = meteoPoints[i].topographicDistance; myPoint.isActive = true; myPoint.isMarked = meteoPoints[i].marked; - myPoint.macroAreaCode = meteoPoints[i].macroAreaCode; if (isEqual(xMin, NODATA)) { diff --git a/meteo/meteoPoint.cpp b/meteo/meteoPoint.cpp index 97ae29fc4..e32f7f411 100644 --- a/meteo/meteoPoint.cpp +++ b/meteo/meteoPoint.cpp @@ -61,7 +61,6 @@ void Crit3DMeteoPoint::clear() this->latInt = NODATA; this->lonInt = NODATA; this->isInsideDem = false; - this->macroAreaCode = NODATA; this->nrObsDataDaysH = 0; this->nrObsDataDaysD = 0; diff --git a/meteo/meteoPoint.h b/meteo/meteoPoint.h index dfbde836f..2c86a6da4 100644 --- a/meteo/meteoPoint.h +++ b/meteo/meteoPoint.h @@ -91,7 +91,6 @@ int latInt; int lonInt; bool isInsideDem; - int macroAreaCode; bool isUTC; bool isForecast; diff --git a/pragaProject/pragaProject.cpp b/pragaProject/pragaProject.cpp index 78ed1b01a..0ba7aa533 100644 --- a/pragaProject/pragaProject.cpp +++ b/pragaProject/pragaProject.cpp @@ -2639,7 +2639,7 @@ bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QL if (interpolationSettings.getUseGlocalDetrending()) { logInfoGUI("Loading macro areas map for glocal detrending..."); - if (loadMacroAreaGlocalMap() < 1) + if (loadGlocalFiles() < 1) return false; } diff --git a/project/project.cpp b/project/project.cpp index 4517beaeb..9db339fef 100644 --- a/project/project.cpp +++ b/project/project.cpp @@ -2142,7 +2142,7 @@ bool Project::loadTopographicDistanceMaps(bool onlyWithData, bool showInfo) return true; } -int Project::loadMacroAreaGlocalMap() +int Project::loadGlocalFiles() { //TODO: add a check for the code values? @@ -2152,40 +2152,116 @@ int Project::loadMacroAreaGlocalMap() //errore? } + /*QString mapsFolder = projectPath + PATH_GLOCAL; + if (! QDir(mapsFolder).exists()) + { + QDir().mkdir(mapsFolder); + }*/ + std::string myError; gis::Crit3DRasterGrid* macroAreasGrid = new gis::Crit3DRasterGrid(); - std::string fileName = mapsFolder.toStdString() + glocalMapName.toStdString(); + std::string fileNameMap = mapsFolder.toStdString() + glocalMapName.toStdString() + "_map"; + std::string fileNameStations = mapsFolder.toStdString() + glocalMapName.toStdString() + "_stations"; - if (!QFile::exists(QString::fromStdString(fileName + ".flt"))) + if (!QFile::exists(QString::fromStdString(fileNameMap + ".flt"))) { - myError = "Couldn't find " + fileName + " file in " + mapsFolder.toStdString(); + myError = "Couldn't find " + fileNameMap + " file in " + mapsFolder.toStdString(); logError(QString::fromStdString(myError)); return false; } - if (gis::readEsriGrid(fileName, macroAreasGrid, myError)) + if (gis::readEsriGrid(fileNameMap, macroAreasGrid, myError)) interpolationSettings.setMacroAreasMap(macroAreasGrid); else { - myError = "Error loading:\n" + fileName; + myError = "Error loading:\n" + fileNameMap; return false; } + //leggi csv aree + std::vector> areaPoints; + loadGlocalStationsCsv(QString::fromStdString(fileNameStations), areaPoints); - //associa il codice a ogni stazione e trova il numero di aree tot - int numAreas = 1; - for (int i=0; i < nrMeteoPoints; i++) + //salva in vettore gli indici dei meteopoints che appartengono a area n, infine salva quel vettore + + std::vector temp; + std::vector myAreas; + Crit3DMacroArea myArea; + + //per ogni area, rappresentata da righe di areaPoints, si fa ciclo sui meteopoints + for (int j = 0; j < areaPoints.size(); j++) + { + for (int i=0; i < nrMeteoPoints; i++) + { + //controlla se id si trova nel vettore areaPoints[j] e salva l'indice del meteopoint + for (int k = 0; k < areaPoints[j].size(); k++) + if (areaPoints[j][k] == std::stoi(meteoPoints[i].id)) temp.push_back(i); + } + //salvataggio di temp dentro il corrispettivo crit3dmacroarea + myArea.setMeteoPoints(temp); + myAreas.push_back(myArea); + temp.clear(); + myArea.clear(); + } + + interpolationSettings.setMacroAreas(myAreas); + + return areaPoints.size(); +} + +bool Project::loadGlocalStationsCsv(QString fileName, std::vector> &areaPoints) +{ + std::string myError; + + if (fileName == "") + { + myError = "Missing csv filename"; + logError(QString::fromStdString(myError)); + return false; + } + areaPoints.clear(); + + if (! QFile(fileName + ".csv").exists() || ! QFileInfo(fileName + ".csv").isFile()) + { + myError = "Missing csv file: " + fileName.toStdString(); + logError(QString::fromStdString(myError)); + return false; + } + + QFile myFile(fileName + ".csv"); + if (!myFile.open(QIODevice::ReadOnly)) { - int myValue = gis::getValueFromXY(*macroAreasGrid, meteoPoints[i].point.utm.x, meteoPoints[i].point.utm.y); + myError = "Open CSV failed: " + fileName.toStdString(); + logError(QString::fromStdString(myError)); + return false; + } - if (myValue > numAreas) - numAreas = myValue; + QTextStream myStream (&myFile); + QList line; - if (myValue != (*macroAreasGrid).header->flag) - meteoPoints[i].macroAreaCode = myValue; + if (myStream.atEnd()) + { + myError = "\nFile is void."; + myFile.close(); + return false; } - return numAreas+1; + int nrLine = 0; + std::vector temp; + while(!myStream.atEnd()) + { + nrLine++; + line = myStream.readLine().split(','); + for (int i = 1; i < line.size(); i++) + temp.push_back(line[i].toInt()); + + areaPoints.resize(line[0].toInt() + 1); + areaPoints[line[0].toInt()] = temp; + temp.clear(); + } + myFile.close(); + + return true; } bool Project::groupCellsInArea(std::vector &areaPoints, unsigned int index, bool isGrid) @@ -2581,7 +2657,7 @@ bool Project::interpolationDemLocalDetrending(meteoVariable myVar, const Crit3DT return true; } -bool Project::interpolationDemGlocalDetrending(meteoVariable myVar, const Crit3DTime& myTime, gis::Crit3DRasterGrid *myRaster) +bool Project::interpolationDemGlocalDetrending(int numAreas, meteoVariable myVar, const Crit3DTime& myTime, gis::Crit3DRasterGrid *myRaster) { if (!getUseDetrendingVar(myVar) || !interpolationSettings.getUseGlocalDetrending()) return false; @@ -2606,13 +2682,6 @@ bool Project::interpolationDemGlocalDetrending(meteoVariable myVar, const Crit3D Crit3DProxyCombination myCombination = interpolationSettings.getSelectedCombination(); interpolationSettings.setCurrentCombination(myCombination); - int elevationPos = NODATA; - for (unsigned int pos=0; pos < interpolationSettings.getCurrentCombination().getProxySize(); pos++) - { - if (getProxyPragaName(interpolationSettings.getProxy(pos)->getName()) == proxyHeight) - elevationPos = pos; - } - //obtain fitting parameters for every macro area if (! preInterpolation(interpolationPoints, &interpolationSettings, meteoSettings, &climateParameters, meteoPoints, nrMeteoPoints, myVar, myTime, errorStdStr)) @@ -2640,9 +2709,14 @@ bool Project::interpolationDemGlocalDetrending(meteoVariable myVar, const Crit3D unsigned row, col; std::vector>> allAreaFittingParameters = interpolationSettings.getAreaParameters(); std::vector allAreaCombinations = interpolationSettings.getAreaCombination(); - std::vector subsetInterpolationPoints; std::vector areaCells; - unsigned numAreas = allAreaFittingParameters.size(); + + int elevationPos = NODATA; + for (unsigned int pos=0; pos < interpolationSettings.getCurrentCombination().getProxySize(); pos++) + { + if (getProxyPragaName(interpolationSettings.getProxy(pos)->getName()) == proxyHeight) + elevationPos = pos; + } if (allAreaFittingParameters.size() < numAreas || allAreaCombinations.size() < numAreas) { @@ -2653,8 +2727,9 @@ bool Project::interpolationDemGlocalDetrending(meteoVariable myVar, const Crit3D for (unsigned areaIndex = 0; areaIndex < numAreas; areaIndex++) { //save all the cells in that specific area in a vector (areaCells) + Crit3DMacroArea myArea = interpolationSettings.getMacroAreas()[areaIndex]; areaCells.clear(); - groupCellsInArea(areaCells, areaIndex, false); //mettere qui calcolo buffer, riutilizzo funzione distanza da area + groupCellsInArea(areaCells, areaIndex, false); //take the parameters+combination for that area std::vector &)>> fittingFunction; @@ -2676,11 +2751,13 @@ bool Project::interpolationDemGlocalDetrending(meteoVariable myVar, const Crit3D //create group of macro area interpolation points interpolationSettings.setFittingFunction(fittingFunction); - subsetInterpolationPoints.clear(); - for (int l = 0; l < interpolationPoints.size(); l++) + std::vector temp = myArea.getMeteoPoints(); + std::vector subsetInterpolationPoints; + for (int l = 0; l < temp.size(); l++) { - if (interpolationPoints[l].macroAreaCode == areaIndex) - subsetInterpolationPoints.push_back(interpolationPoints[l]); + for (int k = 0; k < interpolationPoints.size(); k++) + if (interpolationPoints[k].index == l) + subsetInterpolationPoints.push_back(interpolationPoints[k]); } //detrending (done once for every macro area) @@ -2858,9 +2935,11 @@ bool Project::interpolationDemMain(meteoVariable myVar, const Crit3DTime& myTime if (interpolationSettings.getUseMultipleDetrending()) interpolationSettings.clearFitting(); + int numAreas = 1; if (interpolationSettings.getUseGlocalDetrending()) { - if (loadMacroAreaGlocalMap() < 1) return false; + numAreas = loadGlocalFiles(); + if (numAreas < 1) return false; } // dynamic lapserate if (getUseDetrendingVar(myVar) && interpolationSettings.getUseLocalDetrending()) @@ -2869,7 +2948,7 @@ bool Project::interpolationDemMain(meteoVariable myVar, const Crit3DTime& myTime } else if (getUseDetrendingVar(myVar) && interpolationSettings.getUseGlocalDetrending()) { - return interpolationDemGlocalDetrending(myVar, myTime, myRaster); + return interpolationDemGlocalDetrending(numAreas, myVar, myTime, myRaster); } else { @@ -2919,7 +2998,7 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) unsigned numAreas = 1; if (interpolationSettings.getUseGlocalDetrending()) { - numAreas = loadMacroAreaGlocalMap(); + numAreas = loadGlocalFiles(); if (numAreas < 1) return false; } @@ -3057,20 +3136,22 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) return false; } + int elevationPos = NODATA; + for (unsigned int pos=0; pos < interpolationSettings.getCurrentCombination().getProxySize(); pos++) + { + if (getProxyPragaName(interpolationSettings.getProxy(pos)->getName()) == proxyHeight) + elevationPos = pos; + } + for (unsigned areaIndex = 0; areaIndex < numAreas; areaIndex++) { + Crit3DMacroArea myArea = interpolationSettings.getMacroAreas()[areaIndex]; areaCells.clear(); groupCellsInArea(areaCells, areaIndex, true); std::vector &)> > fittingFunction; interpolationSettings.setFittingParameters(allAreaFittingParameters[areaIndex]); interpolationSettings.setCurrentCombination(allAreaCombinations[areaIndex]); - int elevationPos = NODATA; - for (unsigned int pos=0; pos < interpolationSettings.getCurrentCombination().getProxySize(); pos++) - { - if (getProxyPragaName(interpolationSettings.getProxy(pos)->getName()) == proxyHeight) - elevationPos = pos; - } for (int l = 0; l < allAreaFittingParameters[areaIndex].size(); l++) { @@ -3085,13 +3166,13 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) } interpolationSettings.setFittingFunction(fittingFunction); - + std::vector temp = myArea.getMeteoPoints(); std::vector subsetInterpolationPoints; - - for (int l = 0; l < interpolationPoints.size(); l++) + for (int l = 0; l < temp.size(); l++) { - if (interpolationPoints[l].macroAreaCode == areaIndex) - subsetInterpolationPoints.push_back(interpolationPoints[l]); + for (int k = 0; k < interpolationPoints.size(); k++) + if (interpolationPoints[k].index == l) + subsetInterpolationPoints.push_back(interpolationPoints[k]); } if (elevationPos != NODATA && myCombination.isProxyActive(elevationPos)) @@ -4097,7 +4178,7 @@ void Project::showLocalProxyGraph(gis::Crit3DGeoPoint myPoint) } if (interpolationSettings.getUseGlocalDetrending()) - loadMacroAreaGlocalMap(); + loadGlocalFiles(); localProxyWidget = new Crit3DLocalProxyWidget(myUtm.x, myUtm.y, myZDEM, myZGrid, this->gisSettings, &interpolationSettings, meteoPoints, nrMeteoPoints, currentVariable, currentFrequency, currentDate, currentHour, quality, &qualityInterpolationSettings, meteoSettings, &climateParameters, checkSpatialQuality); return; diff --git a/project/project.h b/project/project.h index 2a9790021..2b1ce4d15 100644 --- a/project/project.h +++ b/project/project.h @@ -260,7 +260,8 @@ bool writeTopographicDistanceMap(int pointIndex, const gis::Crit3DRasterGrid& demMap, QString pathTd); bool loadTopographicDistanceMaps(bool onlyWithData, bool showInfo); void passInterpolatedTemperatureToHumidityPoints(Crit3DTime myTime, Crit3DMeteoSettings *meteoSettings); - int loadMacroAreaGlocalMap(); + int loadGlocalFiles(); + bool loadGlocalStationsCsv(QString fileName, std::vector> &areaPoints); bool groupCellsInArea(std::vector &areaPoints, unsigned index, bool isGrid); bool checkInterpolation(meteoVariable myVar); @@ -269,7 +270,7 @@ bool interpolationDemMain(meteoVariable myVar, const Crit3DTime& myTime, gis::Crit3DRasterGrid *myRaster); bool interpolationDem(meteoVariable myVar, const Crit3DTime& myTime, gis::Crit3DRasterGrid *myRaster); bool interpolationDemLocalDetrending(meteoVariable myVar, const Crit3DTime& myTime, gis::Crit3DRasterGrid *myRaster); - bool interpolationDemGlocalDetrending(meteoVariable myVar, const Crit3DTime& myTime, gis::Crit3DRasterGrid *myRaster); + bool interpolationDemGlocalDetrending(int numAreas, meteoVariable myVar, const Crit3DTime& myTime, gis::Crit3DRasterGrid *myRaster); bool interpolateDemRadiation(const Crit3DTime& myTime, gis::Crit3DRasterGrid *myRaster); bool interpolationOutputPoints(std::vector &interpolationPoints, gis::Crit3DRasterGrid *outputGrid, meteoVariable myVar); diff --git a/proxyWidget/localProxyWidget.cpp b/proxyWidget/localProxyWidget.cpp index 9870cf168..cb9600bfd 100644 --- a/proxyWidget/localProxyWidget.cpp +++ b/proxyWidget/localProxyWidget.cpp @@ -354,6 +354,8 @@ void Crit3DLocalProxyWidget::plot() else if (interpolationSettings->getUseGlocalDetrending()) { areaCode = gis::getValueFromXY(*interpolationSettings->getMacroAreasMap(), x, y); + Crit3DMacroArea myArea = interpolationSettings->getMacroAreas()[areaCode]; + std::vector stations = myArea.getMeteoPoints(); if (detrended.isChecked()) { outInterpolationPoints.clear(); @@ -362,28 +364,31 @@ void Crit3DLocalProxyWidget::plot() interpolationSettings, meteoSettings, climateParam, outInterpolationPoints, checkSpatialQuality, errorStdStr); - for (int k = 0; k < outInterpolationPoints.size(); k++) + for (int k = 0; k < stations.size(); k++) { - if (outInterpolationPoints[k].macroAreaCode == areaCode) - { - subsetInterpolationPoints.push_back(outInterpolationPoints[k]); - } + for (int j = 0; j < outInterpolationPoints.size(); j++) + if (outInterpolationPoints[j].index == stations[k]) + { + subsetInterpolationPoints.push_back(outInterpolationPoints[j]); + } } detrending(subsetInterpolationPoints, interpolationSettings->getSelectedCombination(), interpolationSettings, climateParam, myVar, getCurrentTime()); } else { + outInterpolationPoints.clear(); checkAndPassDataToInterpolation(quality, myVar, meteoPoints, nrMeteoPoints, getCurrentTime(), SQinterpolationSettings, interpolationSettings, meteoSettings, climateParam, outInterpolationPoints, checkSpatialQuality, errorStdStr); - for (int k = 0; k < outInterpolationPoints.size(); k++) + for (int k = 0; k < stations.size(); k++) { - if (outInterpolationPoints[k].macroAreaCode == areaCode) - { - subsetInterpolationPoints.push_back(outInterpolationPoints[k]); - } + for (int j = 0; j < outInterpolationPoints.size(); j++) + if (outInterpolationPoints[j].index == stations[k]) + { + subsetInterpolationPoints.push_back(outInterpolationPoints[j]); + } } } } From 72b5ee9ebc3b94b3812490312bd63798c2973a6c Mon Sep 17 00:00:00 2001 From: Caterina Toscano Date: Mon, 28 Oct 2024 13:12:13 +0100 Subject: [PATCH 07/34] glocal detrending with buffer areas --- interpolation/interpolation.cpp | 16 ++- interpolation/interpolation.h | 2 +- interpolation/interpolationSettings.cpp | 4 +- interpolation/interpolationSettings.h | 6 +- pragaProject/pragaProject.cpp | 7 -- project/project.cpp | 156 ++++++++++++++++++++---- project/project.h | 4 +- proxyWidget/localProxyWidget.cpp | 11 +- proxyWidget/proxyWidget.cpp | 2 +- 9 files changed, 162 insertions(+), 46 deletions(-) diff --git a/interpolation/interpolation.cpp b/interpolation/interpolation.cpp index a3139d781..0cd8d5e86 100644 --- a/interpolation/interpolation.cpp +++ b/interpolation/interpolation.cpp @@ -1723,7 +1723,7 @@ bool multipleDetrendingMain(std::vector &myPoints if (elevationPos != NODATA && mySettings->getCurrentCombination().isProxyActive(elevationPos)) { - if (!multipleDetrendingElevationFitting(elevationPos, myPoints, mySettings, myVar, errorStr)) + if (!multipleDetrendingElevationFitting(elevationPos, myPoints, mySettings, myVar, errorStr, true)) return false; detrendingElevation(elevationPos, myPoints, mySettings); @@ -1739,7 +1739,7 @@ bool multipleDetrendingMain(std::vector &myPoints } bool multipleDetrendingElevationFitting(int elevationPos, std::vector &myPoints, - Crit3DInterpolationSettings* mySettings, meteoVariable myVar, std::string &errorStr) + Crit3DInterpolationSettings* mySettings, meteoVariable myVar, std::string &errorStr, bool isWeighted) { mySettings->getProxy(elevationPos)->setRegressionR2(NODATA); if (! getUseDetrendingVar(myVar)) return true; @@ -1799,8 +1799,13 @@ bool multipleDetrendingElevationFitting(int elevationPos, std::vector > firstGuessCombinations = mySettings->getProxy(elevationPos)->getFirstGuessCombinations(); // multiple non linear fitting - double R2 = interpolation::bestFittingMarquardt_nDimension_singleFunction(*(myFunc.target&)>()), 400, 4, parametersMin, parametersMax, parameters, parametersDelta, + double R2 = NODATA; + if (isWeighted) + R2 = interpolation::bestFittingMarquardt_nDimension_singleFunction(*(myFunc.target&)>()), 400, 4, parametersMin, parametersMax, parameters, parametersDelta, stepSize, numSteps, 1000, 0.002, 0.005, predictors, predictands, weights,firstGuessCombinations); + else + R2 = interpolation::bestFittingMarquardt_nDimension_singleFunction(*(myFunc.target&)>()), 400, 4, parametersMin, parametersMax, parameters, parametersDelta, + stepSize, numSteps, 1000, 0.002, 0.005, predictors, predictands,firstGuessCombinations); mySettings->getProxy(elevationPos)->setRegressionR2(R2); @@ -2088,7 +2093,10 @@ bool glocalDetrendingFitting(std::vector &myPoint for (int k = 0; k < myPoints.size(); k++) if (myPoints[k].index == temp[l]) + { subsetPoints.push_back(myPoints[k]); + break; + } } mySettings->setCurrentCombination(mySettings->getSelectedCombination()); @@ -2096,7 +2104,7 @@ bool glocalDetrendingFitting(std::vector &myPoint if (elevationPos != NODATA && mySettings->getCurrentCombination().isProxyActive(elevationPos)) { - if (!multipleDetrendingElevationFitting(elevationPos, subsetPoints, mySettings, myVar, errorStr)) return false; + if (!multipleDetrendingElevationFitting(elevationPos, subsetPoints, mySettings, myVar, errorStr, false)) return false; detrendingElevation(elevationPos, subsetPoints, mySettings); diff --git a/interpolation/interpolation.h b/interpolation/interpolation.h index 63b6c6a22..99fdfcf57 100644 --- a/interpolation/interpolation.h +++ b/interpolation/interpolation.h @@ -74,7 +74,7 @@ meteoVariable myVar, std::string &errorStr); bool multipleDetrendingElevationFitting(int elevationPos, std::vector &myPoints, - Crit3DInterpolationSettings* mySettings, meteoVariable myVar, std::string &errorStr); + Crit3DInterpolationSettings* mySettings, meteoVariable myVar, std::string &errorStr, bool isWeighted); void detrendingElevation(int elevationPos, std::vector &myPoints, Crit3DInterpolationSettings* mySettings); diff --git a/interpolation/interpolationSettings.cpp b/interpolation/interpolationSettings.cpp index e6c5860e5..f286df8d0 100644 --- a/interpolation/interpolationSettings.cpp +++ b/interpolation/interpolationSettings.cpp @@ -1039,13 +1039,13 @@ std::vector Crit3DMacroArea::getMeteoPoints() return meteoPoints; } -void Crit3DMacroArea::setAreaCells (std::vector myCells) +void Crit3DMacroArea::setAreaCells (std::vector myCells) { areaCells = myCells; return; } -std::vector Crit3DMacroArea::getAreaCells() +std::vector Crit3DMacroArea::getAreaCells() { return areaCells; } diff --git a/interpolation/interpolationSettings.h b/interpolation/interpolationSettings.h index cc665edd8..4371d82bb 100644 --- a/interpolation/interpolationSettings.h +++ b/interpolation/interpolationSettings.h @@ -138,7 +138,7 @@ Crit3DProxyCombination areaCombination; std::vector> areaParameters; std::vector meteoPoints; - std::vector areaCells; + std::vector areaCells; public: Crit3DMacroArea(); @@ -149,8 +149,8 @@ std::vector> getParameters(); void setCombination (Crit3DProxyCombination myCombination); Crit3DProxyCombination getCombination(); - void setAreaCells (std::vector myCells); - std::vector getAreaCells(); + void setAreaCells (std::vector myCells); + std::vector getAreaCells(); void clear(); }; diff --git a/pragaProject/pragaProject.cpp b/pragaProject/pragaProject.cpp index 0ba7aa533..27b4aeecc 100644 --- a/pragaProject/pragaProject.cpp +++ b/pragaProject/pragaProject.cpp @@ -2636,13 +2636,6 @@ bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QL return false; } - if (interpolationSettings.getUseGlocalDetrending()) - { - logInfoGUI("Loading macro areas map for glocal detrending..."); - if (loadGlocalFiles() < 1) - return false; - } - // save also time aggregated variables foreach (myVar, aggrVariables) varToSave.push_back(myVar); diff --git a/project/project.cpp b/project/project.cpp index c13a0bd65..2e69f48dc 100644 --- a/project/project.cpp +++ b/project/project.cpp @@ -2234,7 +2234,7 @@ bool Project::glocalWeightsMaps(float windowWidth) return true; } -int Project::loadGlocalFiles() +bool Project::loadGlocalAreasMap() { //TODO: add a check for the code values? @@ -2253,7 +2253,6 @@ int Project::loadGlocalFiles() std::string myError; gis::Crit3DRasterGrid* macroAreasGrid = new gis::Crit3DRasterGrid(); std::string fileNameMap = mapsFolder.toStdString() + glocalMapName.toStdString() + "_map"; - std::string fileNameStations = mapsFolder.toStdString() + glocalMapName.toStdString() + "_stations"; if (!QFile::exists(QString::fromStdString(fileNameMap + ".flt"))) { @@ -2270,12 +2269,20 @@ int Project::loadGlocalFiles() return false; } + return true; +} + +int Project::loadGlocalStationsAndCells(bool isGrid) +{ //leggi csv aree + QString mapsFolder = defaultPath + PATH_GEO; + std::string fileNameStations = mapsFolder.toStdString() + glocalMapName.toStdString() + "_stations"; std::vector> areaPoints; + std::string myError; + loadGlocalStationsCsv(QString::fromStdString(fileNameStations), areaPoints); //salva in vettore gli indici dei meteopoints che appartengono a area n, infine salva quel vettore - std::vector temp; std::vector myAreas; Crit3DMacroArea myArea; @@ -2283,11 +2290,15 @@ int Project::loadGlocalFiles() //per ogni area, rappresentata da righe di areaPoints, si fa ciclo sui meteopoints for (int j = 0; j < areaPoints.size(); j++) { - for (int i=0; i < nrMeteoPoints; i++) + for (int k = 0; k < areaPoints[j].size(); k++) { //controlla se id si trova nel vettore areaPoints[j] e salva l'indice del meteopoint - for (int k = 0; k < areaPoints[j].size(); k++) - if (areaPoints[j][k] == std::stoi(meteoPoints[i].id)) temp.push_back(i); + for (int i=0; i < nrMeteoPoints; i++) + if (areaPoints[j][k] == std::stoi(meteoPoints[i].id)) + { + temp.push_back(i); + break; + } } //salvataggio di temp dentro il corrispettivo crit3dmacroarea myArea.setMeteoPoints(temp); @@ -2296,11 +2307,86 @@ int Project::loadGlocalFiles() myArea.clear(); } + //assegnazione pesi a ogni cella di ogni area + if (!loadGlocalWeightMaps(myAreas, isGrid)) + { + myError = "error while loading the glocal weight maps."; + return false; + } + interpolationSettings.setMacroAreas(myAreas); return areaPoints.size(); } +bool Project::loadGlocalWeightMaps(std::vector &myAreas, bool isGrid) +{ + QString mapsFolder = projectPath + PATH_GLOCAL; + if (! QDir(mapsFolder).exists()) + glocalWeightsMaps(10000); + + std::string myError; + gis::Crit3DRasterGrid* macroAreasGrid = new gis::Crit3DRasterGrid(); + std::string fileName = mapsFolder.toStdString() + "glocalWeight_"; + + std::vector singleCell(2); + std::vector areaCells; + int nrCols, nrRows; + double myLat, myLon, myX, myY; + float myValue = NODATA; + + if (!isGrid) + { + nrCols = DEM.header->nrCols; + nrRows = DEM.header->nrCols; + } + else + { + nrCols = meteoGridDbHandler->gridStructure().header().nrCols; + nrRows = meteoGridDbHandler->gridStructure().header().nrCols; + } + + for (int i = 0; i < myAreas.size(); i++) + { + if (QFile::exists(QString::fromStdString(fileName) + QString::number(i) + ".flt")) + { + if (gis::readEsriGrid(fileName + std::to_string(i), macroAreasGrid, myError)) + { + for (int row = 0; row < nrRows; row++) + { + for (int col = 0; col < nrCols; col++) + { + if (isGrid) + { + gis::getLatLonFromRowCol(meteoGridDbHandler->gridStructure().header(), row, col, &myLat, &myLon); + gis::getUtmFromLatLon(gisSettings, myLat, myLon, &myX, &myY); + } + else + DEM.getXY(row, col, myX, myY); + + myValue = macroAreasGrid->getValueFromXY(myX, myY); + + if (!isEqual(myValue, NODATA) && !isEqual(myValue, 0)) + { + areaCells.push_back(row*nrCols+col); + areaCells.push_back(myValue); + } + } + } + myAreas[i].setAreaCells(areaCells); + areaCells.clear(); + } + else { + macroAreasGrid->clear(); + return false; + } + } + } + macroAreasGrid->clear(); + + return true; +} + bool Project::loadGlocalStationsCsv(QString fileName, std::vector> &areaPoints) { std::string myError; @@ -2801,7 +2887,7 @@ bool Project::interpolationDemGlocalDetrending(int numAreas, meteoVariable myVar unsigned row, col; std::vector>> allAreaFittingParameters = interpolationSettings.getAreaParameters(); std::vector allAreaCombinations = interpolationSettings.getAreaCombination(); - std::vector areaCells; + std::vector areaCells; int elevationPos = NODATA; for (unsigned int pos=0; pos < interpolationSettings.getCurrentCombination().getProxySize(); pos++) @@ -2820,8 +2906,8 @@ bool Project::interpolationDemGlocalDetrending(int numAreas, meteoVariable myVar { //save all the cells in that specific area in a vector (areaCells) Crit3DMacroArea myArea = interpolationSettings.getMacroAreas()[areaIndex]; - areaCells.clear(); - groupCellsInArea(areaCells, areaIndex, false); + areaCells = myArea.getAreaCells(); + //groupCellsInArea(areaCells, areaIndex, false); //take the parameters+combination for that area std::vector &)>> fittingFunction; @@ -2848,7 +2934,7 @@ bool Project::interpolationDemGlocalDetrending(int numAreas, meteoVariable myVar for (int l = 0; l < temp.size(); l++) { for (int k = 0; k < interpolationPoints.size(); k++) - if (interpolationPoints[k].index == l) + if (interpolationPoints[k].index == temp[l]) subsetInterpolationPoints.push_back(interpolationPoints[k]); } @@ -2858,11 +2944,10 @@ bool Project::interpolationDemGlocalDetrending(int numAreas, meteoVariable myVar detrendingOtherProxies(elevationPos, subsetInterpolationPoints, &interpolationSettings); - for (unsigned cellIndex = 0; cellIndex < areaCells.size(); cellIndex++) + for (unsigned cellIndex = 0; cellIndex < areaCells.size(); cellIndex = cellIndex + 2) { - row = unsigned(areaCells[cellIndex]/DEM.header->nrCols); - col = areaCells[cellIndex]%DEM.header->nrCols; + col = int(areaCells[cellIndex])%DEM.header->nrCols; float z = DEM.value[row][col]; if (! isEqual(z, myHeader.flag)) @@ -2871,9 +2956,14 @@ bool Project::interpolationDemGlocalDetrending(int numAreas, meteoVariable myVar getProxyValuesXY(x, y, &interpolationSettings, proxyValues); - myRaster->value[row][col] = interpolate(subsetInterpolationPoints, &interpolationSettings, meteoSettings, - myVar, x, y, z, proxyValues, true); + if (isEqual(myRaster->value[row][col], NODATA)) + myRaster->value[row][col] = interpolate(subsetInterpolationPoints, &interpolationSettings, meteoSettings, + myVar, x, y, z, proxyValues, true)*areaCells[cellIndex+1]; + else + myRaster->value[row][col] += interpolate(subsetInterpolationPoints, &interpolationSettings, meteoSettings, + myVar, x, y, z, proxyValues, true)*areaCells[cellIndex+1]; } + } } } @@ -3030,8 +3120,13 @@ bool Project::interpolationDemMain(meteoVariable myVar, const Crit3DTime& myTime int numAreas = 1; if (interpolationSettings.getUseGlocalDetrending()) { - numAreas = loadGlocalFiles(); - if (numAreas < 1) return false; + if (loadGlocalAreasMap()) + { + numAreas = loadGlocalStationsAndCells(false); + if (numAreas < 1) return false; + } + else + return false; } // dynamic lapserate if (getUseDetrendingVar(myVar) && interpolationSettings.getUseLocalDetrending()) @@ -3090,8 +3185,13 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) unsigned numAreas = 1; if (interpolationSettings.getUseGlocalDetrending()) { - numAreas = loadGlocalFiles(); - if (numAreas < 1) return false; + if (loadGlocalAreasMap()) + { + numAreas = loadGlocalStationsAndCells(true); + if (numAreas < 1) return false; + } + else + return false; } // check quality and pass data to interpolation @@ -3220,7 +3320,7 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) unsigned row, col; std::vector>> allAreaFittingParameters = interpolationSettings.getAreaParameters(); std::vector allAreaCombinations = interpolationSettings.getAreaCombination(); - std::vector areaCells; + std::vector areaCells; if (allAreaCombinations.size() != numAreas) { @@ -3238,8 +3338,8 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) for (unsigned areaIndex = 0; areaIndex < numAreas; areaIndex++) { Crit3DMacroArea myArea = interpolationSettings.getMacroAreas()[areaIndex]; - areaCells.clear(); - groupCellsInArea(areaCells, areaIndex, true); + areaCells = myArea.getAreaCells(); + //groupCellsInArea(areaCells, areaIndex, true); std::vector &)> > fittingFunction; interpolationSettings.setFittingParameters(allAreaFittingParameters[areaIndex]); interpolationSettings.setCurrentCombination(allAreaCombinations[areaIndex]); @@ -3274,10 +3374,10 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) unsigned nrCols = meteoGridDbHandler->meteoGrid()->gridStructure().header().nrCols; - for (unsigned cellIndex = 0; cellIndex < areaCells.size(); cellIndex++) + for (unsigned cellIndex = 0; cellIndex < areaCells.size(); cellIndex = cellIndex + 2) { row = unsigned(areaCells[cellIndex]/nrCols); - col = areaCells[cellIndex]%nrCols; + col = int(areaCells[cellIndex])%nrCols; if (meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->active) @@ -3306,7 +3406,8 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) proxyIndex++; } } - interpolatedValue = interpolate(subsetInterpolationPoints, &interpolationSettings, meteoSettings, myVar, myX, myY, myZ, proxyValues, true); + interpolatedValue += interpolate(subsetInterpolationPoints, &interpolationSettings, meteoSettings, myVar, myX, myY, myZ, proxyValues, true) + * areaCells[cellIndex + 1]; } if (freq == hourly) @@ -4270,7 +4371,10 @@ void Project::showLocalProxyGraph(gis::Crit3DGeoPoint myPoint) } if (interpolationSettings.getUseGlocalDetrending()) - loadGlocalFiles(); + { + if (loadGlocalAreasMap()) + loadGlocalStationsAndCells(false); + } localProxyWidget = new Crit3DLocalProxyWidget(myUtm.x, myUtm.y, myZDEM, myZGrid, this->gisSettings, &interpolationSettings, meteoPoints, nrMeteoPoints, currentVariable, currentFrequency, currentDate, currentHour, quality, &qualityInterpolationSettings, meteoSettings, &climateParameters, checkSpatialQuality); return; diff --git a/project/project.h b/project/project.h index c39ef7a4a..46302b3d1 100644 --- a/project/project.h +++ b/project/project.h @@ -260,7 +260,9 @@ bool writeTopographicDistanceMap(int pointIndex, const gis::Crit3DRasterGrid& demMap, QString pathTd); bool loadTopographicDistanceMaps(bool onlyWithData, bool showInfo); void passInterpolatedTemperatureToHumidityPoints(Crit3DTime myTime, Crit3DMeteoSettings *meteoSettings); - int loadGlocalFiles(); + bool loadGlocalAreasMap(); + int loadGlocalStationsAndCells(bool isGrid); + bool loadGlocalWeightMaps(std::vector &myAreas, bool isGrid); bool loadGlocalStationsCsv(QString fileName, std::vector> &areaPoints); bool groupCellsInArea(std::vector &areaPoints, unsigned index, bool isGrid); bool glocalWeightsMaps(float windowWidth); diff --git a/proxyWidget/localProxyWidget.cpp b/proxyWidget/localProxyWidget.cpp index cb9600bfd..93ef34072 100644 --- a/proxyWidget/localProxyWidget.cpp +++ b/proxyWidget/localProxyWidget.cpp @@ -353,7 +353,10 @@ void Crit3DLocalProxyWidget::plot() } else if (interpolationSettings->getUseGlocalDetrending()) { + areaCode = gis::getValueFromXY(*interpolationSettings->getMacroAreasMap(), x, y); + if (areaCode < interpolationSettings->getMacroAreas().size()) + { Crit3DMacroArea myArea = interpolationSettings->getMacroAreas()[areaCode]; std::vector stations = myArea.getMeteoPoints(); if (detrended.isChecked()) @@ -391,6 +394,7 @@ void Crit3DLocalProxyWidget::plot() } } } + } } QList pointListPrimary, pointListSecondary, pointListSupplemental, pointListMarked; QMap< QString, QPointF > idPointMap1; @@ -565,7 +569,12 @@ void Crit3DLocalProxyWidget::modelLRClicked(int toggled) setMultipleDetrendingHeightTemperatureRange(interpolationSettings); interpolationSettings->setCurrentCombination(interpolationSettings->getSelectedCombination()); interpolationSettings->clearFitting(); - if (! multipleDetrendingElevationFitting(proxyPos, subsetInterpolationPoints, interpolationSettings, myVar, errorStr)) return; + if (interpolationSettings->getUseLocalDetrending()) + { + if (! multipleDetrendingElevationFitting(proxyPos, subsetInterpolationPoints, interpolationSettings, myVar, errorStr, true)) return; + } + else if (interpolationSettings->getUseGlocalDetrending()) + if (! multipleDetrendingElevationFitting(proxyPos, subsetInterpolationPoints, interpolationSettings, myVar, errorStr, false)) return; std::vector> parameters = interpolationSettings->getFittingParameters(); diff --git a/proxyWidget/proxyWidget.cpp b/proxyWidget/proxyWidget.cpp index 3409034b4..4140ec1c3 100644 --- a/proxyWidget/proxyWidget.cpp +++ b/proxyWidget/proxyWidget.cpp @@ -481,7 +481,7 @@ void Crit3DProxyWidget::modelLRClicked(int toggled) setMultipleDetrendingHeightTemperatureRange(interpolationSettings); interpolationSettings->setCurrentCombination(interpolationSettings->getSelectedCombination()); interpolationSettings->clearFitting(); - if (! multipleDetrendingElevationFitting(proxyPos, outInterpolationPoints, interpolationSettings, myVar, errorStr)) return; + if (! multipleDetrendingElevationFitting(proxyPos, outInterpolationPoints, interpolationSettings, myVar, errorStr, true)) return; std::vector> parameters = interpolationSettings->getFittingParameters(); From 7cd4b0d001119633599bf0fec543811462376e0c Mon Sep 17 00:00:00 2001 From: cate-to Date: Mon, 28 Oct 2024 15:49:29 +0100 Subject: [PATCH 08/34] changes in glocal detrending (grid) --- interpolation/interpolation.cpp | 4 ++-- project/project.cpp | 26 +++++++++++++++++++------- 2 files changed, 21 insertions(+), 9 deletions(-) diff --git a/interpolation/interpolation.cpp b/interpolation/interpolation.cpp index 0cd8d5e86..c3671b467 100644 --- a/interpolation/interpolation.cpp +++ b/interpolation/interpolation.cpp @@ -2071,7 +2071,7 @@ bool glocalDetrendingFitting(std::vector &myPoint int areasNr = macroAreas.size(); int elevationPos = NODATA; - for (unsigned int pos=0; pos < mySettings->getCurrentCombination().getProxySize(); pos++) + for (unsigned int pos=0; pos < mySettings->getSelectedCombination().getProxySize(); pos++) { if (getProxyPragaName(mySettings->getProxy(pos)->getName()) == proxyHeight) elevationPos = pos; @@ -2102,7 +2102,7 @@ bool glocalDetrendingFitting(std::vector &myPoint mySettings->setCurrentCombination(mySettings->getSelectedCombination()); mySettings->clearFitting(); - if (elevationPos != NODATA && mySettings->getCurrentCombination().isProxyActive(elevationPos)) + if (elevationPos != NODATA && mySettings->getSelectedCombination().isProxyActive(elevationPos)) { if (!multipleDetrendingElevationFitting(elevationPos, subsetPoints, mySettings, myVar, errorStr, false)) return false; diff --git a/project/project.cpp b/project/project.cpp index 2e69f48dc..b045d3013 100644 --- a/project/project.cpp +++ b/project/project.cpp @@ -2329,7 +2329,7 @@ bool Project::loadGlocalWeightMaps(std::vector &myAreas, bool i gis::Crit3DRasterGrid* macroAreasGrid = new gis::Crit3DRasterGrid(); std::string fileName = mapsFolder.toStdString() + "glocalWeight_"; - std::vector singleCell(2); + std::vector areaCells; int nrCols, nrRows; double myLat, myLon, myX, myY; @@ -2343,7 +2343,7 @@ bool Project::loadGlocalWeightMaps(std::vector &myAreas, bool i else { nrCols = meteoGridDbHandler->gridStructure().header().nrCols; - nrRows = meteoGridDbHandler->gridStructure().header().nrCols; + nrRows = meteoGridDbHandler->gridStructure().header().nrRows; } for (int i = 0; i < myAreas.size(); i++) @@ -3321,6 +3321,7 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) std::vector>> allAreaFittingParameters = interpolationSettings.getAreaParameters(); std::vector allAreaCombinations = interpolationSettings.getAreaCombination(); std::vector areaCells; + float myValue = NODATA; if (allAreaCombinations.size() != numAreas) { @@ -3406,7 +3407,7 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) proxyIndex++; } } - interpolatedValue += interpolate(subsetInterpolationPoints, &interpolationSettings, meteoSettings, myVar, myX, myY, myZ, proxyValues, true) + interpolatedValue = interpolate(subsetInterpolationPoints, &interpolationSettings, meteoSettings, myVar, myX, myY, myZ, proxyValues, true) * areaCells[cellIndex + 1]; } @@ -3415,16 +3416,27 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) if (meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->nrObsDataDaysH == 0) meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->initializeObsDataH(1, 1, myTime.date); - meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->setMeteoPointValueH(myTime.date, myTime.getHour(), myTime.getMinutes(), myVar, float(interpolatedValue)); - meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->currentValue = float(interpolatedValue); + myValue = meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->getMeteoPointValueH(myTime.date, myTime.getHour(), myTime.getMinutes(), myVar); + + if (isEqual(myValue, NODATA)) + myValue = 0; + + meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->setMeteoPointValueH(myTime.date, myTime.getHour(), myTime.getMinutes(), myVar, float(interpolatedValue+myValue)); + meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->currentValue = float(interpolatedValue+myValue); + } else if (freq == daily) { if (meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->nrObsDataDaysD == 0) meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->initializeObsDataD(1, myTime.date); - meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->setMeteoPointValueD(myTime.date, myVar, float(interpolatedValue)); - meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->currentValue = float(interpolatedValue); + myValue = meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->getMeteoPointValueD(myTime.date, myVar); + + if (isEqual(myValue, NODATA)) + myValue = 0; + + meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->setMeteoPointValueD(myTime.date, myVar, float(interpolatedValue+myValue)); + meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->currentValue = float(interpolatedValue+myValue); } } From 166f46dfd663ac85e7635fefb7930ba6d2134206 Mon Sep 17 00:00:00 2001 From: cate-to Date: Tue, 29 Oct 2024 09:43:40 +0100 Subject: [PATCH 09/34] fixed grid glocal detrending --- project/project.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/project/project.cpp b/project/project.cpp index b045d3013..cd9333a9e 100644 --- a/project/project.cpp +++ b/project/project.cpp @@ -2332,7 +2332,7 @@ bool Project::loadGlocalWeightMaps(std::vector &myAreas, bool i std::vector areaCells; int nrCols, nrRows; - double myLat, myLon, myX, myY; + double myX, myY; float myValue = NODATA; if (!isGrid) @@ -2358,8 +2358,8 @@ bool Project::loadGlocalWeightMaps(std::vector &myAreas, bool i { if (isGrid) { - gis::getLatLonFromRowCol(meteoGridDbHandler->gridStructure().header(), row, col, &myLat, &myLon); - gis::getUtmFromLatLon(gisSettings, myLat, myLon, &myX, &myY); + myX = meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->point.utm.x; + myY = meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->point.utm.y; } else DEM.getXY(row, col, myX, myY); From 8644929050bd0b630332eac137b063b59efbfd41 Mon Sep 17 00:00:00 2001 From: cate-to Date: Tue, 29 Oct 2024 11:41:37 +0100 Subject: [PATCH 10/34] fix glocal detrending for multiple proxies --- interpolation/interpolation.cpp | 12 +++++++++++- project/project.cpp | 1 - 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/interpolation/interpolation.cpp b/interpolation/interpolation.cpp index c3671b467..dd8fd6f41 100644 --- a/interpolation/interpolation.cpp +++ b/interpolation/interpolation.cpp @@ -2039,6 +2039,9 @@ void detrendingOtherProxies(int elevationPos, std::vector> fullParameters = parameters; + std::vector&)>> fullFunc = myFunc; + // detrending float detrendValue; for (int i = 0; i < myPoints.size(); i++) @@ -2050,13 +2053,20 @@ void detrendingOtherProxies(int elevationPos, std::vectorgetCurrentCombination().isProxyActive(pos) && mySettings->getCurrentCombination().isProxySignificant(pos)) { proxyValue = myPoints[i].getProxyValue(pos); - proxyValues.push_back(double(proxyValue)); + if (!isEqual(proxyValue, NODATA)) + proxyValues.push_back(double(proxyValue)); + else { + parameters.erase(parameters.begin()+proxyValues.size()); + myFunc.erase(myFunc.begin()+proxyValues.size()); + } } } detrendValue = float(functionSum(myFunc, proxyValues, parameters)); myPoints[i].value -= detrendValue; + myFunc = fullFunc; + parameters = fullParameters; } } diff --git a/project/project.cpp b/project/project.cpp index cd9333a9e..adc04782e 100644 --- a/project/project.cpp +++ b/project/project.cpp @@ -2963,7 +2963,6 @@ bool Project::interpolationDemGlocalDetrending(int numAreas, meteoVariable myVar myRaster->value[row][col] += interpolate(subsetInterpolationPoints, &interpolationSettings, meteoSettings, myVar, x, y, z, proxyValues, true)*areaCells[cellIndex+1]; } - } } } From 4f17a0a90a6c2ccb77bcc2441972aead7baf1fad Mon Sep 17 00:00:00 2001 From: cate-to Date: Tue, 29 Oct 2024 14:54:53 +0100 Subject: [PATCH 11/34] fix grid glocal detrending --- interpolation/interpolation.cpp | 6 +- project/project.cpp | 6 +- proxyWidget/localProxyWidget.cpp | 133 +++++++++++++++++++------------ 3 files changed, 92 insertions(+), 53 deletions(-) diff --git a/interpolation/interpolation.cpp b/interpolation/interpolation.cpp index dd8fd6f41..21fc5dd4a 100644 --- a/interpolation/interpolation.cpp +++ b/interpolation/interpolation.cpp @@ -1984,7 +1984,10 @@ bool multipleDetrendingOtherProxiesFitting(int elevationPos, std::vector getUseLocalDetrending() && othersPoints.size() < mySettings->getMinPointsLocalDetrending()) @@ -2037,6 +2040,7 @@ void detrendingOtherProxies(int elevationPos, std::vector> fullParameters = parameters; diff --git a/project/project.cpp b/project/project.cpp index adc04782e..894665dfe 100644 --- a/project/project.cpp +++ b/project/project.cpp @@ -2918,7 +2918,7 @@ bool Project::interpolationDemGlocalDetrending(int numAreas, meteoVariable myVar for (int l = 0; l < allAreaFittingParameters[areaIndex].size(); l++) { if (allAreaFittingParameters[areaIndex][l].size() == 2) - fittingFunction.push_back(functionLinear); + fittingFunction.push_back(functionLinear_intercept); else if (allAreaFittingParameters[areaIndex][l].size() == 4) fittingFunction.push_back(lapseRatePiecewise_two); else if (allAreaFittingParameters[areaIndex][l].size() == 5) @@ -3348,7 +3348,7 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) for (int l = 0; l < allAreaFittingParameters[areaIndex].size(); l++) { if (allAreaFittingParameters[areaIndex][l].size() == 2) - fittingFunction.push_back(functionLinear); + fittingFunction.push_back(functionLinear_intercept); else if (allAreaFittingParameters[areaIndex][l].size() == 4) fittingFunction.push_back(lapseRatePiecewise_two); else if (allAreaFittingParameters[areaIndex][l].size() == 5) @@ -3363,7 +3363,7 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) for (int l = 0; l < temp.size(); l++) { for (int k = 0; k < interpolationPoints.size(); k++) - if (interpolationPoints[k].index == l) + if (interpolationPoints[k].index == temp[l]) subsetInterpolationPoints.push_back(interpolationPoints[k]); } diff --git a/proxyWidget/localProxyWidget.cpp b/proxyWidget/localProxyWidget.cpp index 93ef34072..2e3b76a66 100644 --- a/proxyWidget/localProxyWidget.cpp +++ b/proxyWidget/localProxyWidget.cpp @@ -563,78 +563,113 @@ void Crit3DLocalProxyWidget::modelLRClicked(int toggled) if (toggled && subsetInterpolationPoints.size() != 0) { - if (comboAxisX.currentText() == "elevation") + int elevationPos = NODATA; + for (unsigned int pos=0; pos < interpolationSettings->getSelectedCombination().getProxySize(); pos++) { + if (getProxyPragaName(interpolationSettings->getProxy(pos)->getName()) == proxyHeight) + elevationPos = pos; + } std::string errorStr; setMultipleDetrendingHeightTemperatureRange(interpolationSettings); interpolationSettings->setCurrentCombination(interpolationSettings->getSelectedCombination()); interpolationSettings->clearFitting(); if (interpolationSettings->getUseLocalDetrending()) { - if (! multipleDetrendingElevationFitting(proxyPos, subsetInterpolationPoints, interpolationSettings, myVar, errorStr, true)) return; + if (! multipleDetrendingElevationFitting(elevationPos, subsetInterpolationPoints, interpolationSettings, myVar, errorStr, true)) return; } else if (interpolationSettings->getUseGlocalDetrending()) - if (! multipleDetrendingElevationFitting(proxyPos, subsetInterpolationPoints, interpolationSettings, myVar, errorStr, false)) return; + if (! multipleDetrendingElevationFitting(elevationPos, subsetInterpolationPoints, interpolationSettings, myVar, errorStr, false)) return; std::vector> parameters = interpolationSettings->getFittingParameters(); - if (!parameters.empty()) + if (comboAxisX.currentText() == "elevation") { - if (parameters.front().size() > 3 && parameters.front().size() < 7) + if (!parameters.empty()) { - xMin = getZmin(subsetInterpolationPoints); - xMax = getZmax(subsetInterpolationPoints); - - if (interpolationSettings->getUseMultipleDetrending()) + if (parameters.front().size() > 3 && parameters.front().size() < 7) { - std::vector xVector; - for (int m = xMin; m < xMax; m += 5) - xVector.push_back(m); + xMin = getZmin(subsetInterpolationPoints); + xMax = getZmax(subsetInterpolationPoints); - for (int p = 0; p < int(xVector.size()); p++) + if (interpolationSettings->getUseMultipleDetrending()) { - point.setX(xVector[p]); - if (parameters.front().size() == 4) - point.setY(lapseRatePiecewise_two(xVector[p], parameters.front())); - else if (parameters.front().size() == 5) - point.setY(lapseRatePiecewise_three(xVector[p], parameters.front())); - else if (parameters.front().size() == 6) - point.setY(lapseRatePiecewise_three_free(xVector[p], parameters.front())); - point_vector.append(point); + std::vector xVector; + for (int m = xMin; m < xMax; m += 5) + xVector.push_back(m); + + for (int p = 0; p < int(xVector.size()); p++) + { + point.setX(xVector[p]); + if (parameters.front().size() == 4) + point.setY(lapseRatePiecewise_two(xVector[p], parameters.front())); + else if (parameters.front().size() == 5) + point.setY(lapseRatePiecewise_three(xVector[p], parameters.front())); + else if (parameters.front().size() == 6) + point.setY(lapseRatePiecewise_three_free(xVector[p], parameters.front())); + point_vector.append(point); + } + + if (interpolationSettings->getProxy(elevationPos)->getRegressionR2() != NODATA) + r2.setText(QString("%1").arg(interpolationSettings->getProxy(elevationPos)->getRegressionR2(), 0, 'f', 3)); + + if (parameters.front().size() > 3) + { + par0.setText(QString("%1").arg(parameters.front()[0], 0, 'f', 4)); + par1.setText(QString("%1").arg(parameters.front()[1], 0, 'f', 4)); + par2.setText(QString("%1").arg(parameters.front()[2], 0, 'f', 4)); + par3.setText(QString("%1").arg(parameters.front()[3], 0, 'f', 4)); + } + if (parameters.front().size() > 4) + par4.setText(QString("%1").arg(parameters.front()[4], 0, 'f', 4)); + if (parameters.front().size() > 5) + par5.setText(QString("%1").arg(parameters.front()[5], 0, 'f', 4)); } - - if (interpolationSettings->getProxy(proxyPos)->getRegressionR2() != NODATA) - r2.setText(QString("%1").arg(interpolationSettings->getProxy(proxyPos)->getRegressionR2(), 0, 'f', 3)); - - if (parameters.front().size() > 3) - { - par0.setText(QString("%1").arg(parameters.front()[0], 0, 'f', 4)); - par1.setText(QString("%1").arg(parameters.front()[1], 0, 'f', 4)); - par2.setText(QString("%1").arg(parameters.front()[2], 0, 'f', 4)); - par3.setText(QString("%1").arg(parameters.front()[3], 0, 'f', 4)); - } - if (parameters.front().size() > 4) - par4.setText(QString("%1").arg(parameters.front()[4], 0, 'f', 4)); - if (parameters.front().size() > 5) - par5.setText(QString("%1").arg(parameters.front()[5], 0, 'f', 4)); } } } - } - else - { - //TODO: local proxy graph for other proxies - /*interpolationSettings->setCurrentCombination(interpolationSettings->getSelectedCombination()); - interpolationSettings->clearFitting(); - std::string errorStr; - if (! multipleDetrendingOtherProxiesFitting(proxyPos, subsetInterpolationPoints, interpolationSettings, myVar, errorStr)) return; + else + { + //TODO: OTHER PROXIES + /*std::string errorStr; + //detrendingElevation(elevationPos, outInterpolationPoints, interpolationSettings); + if (! multipleDetrendingOtherProxiesFitting(elevationPos, subsetInterpolationPoints, interpolationSettings, myVar, errorStr)) return; - std::vector> parameters = interpolationSettings->getFittingParameters(); - if (!interpolationSettings->getFittingParameters().empty() && !interpolationSettings->getFittingFunction().empty()) - detrendingOtherProxies(3, outInterpolationPoints, interpolationSettings); - */ + parameters = interpolationSettings->getFittingParameters(); + Crit3DProxyCombination myCombination = interpolationSettings->getCurrentCombination(); + + int pos = 0; + for (int i = 0; i < proxyPos + 1; i++) + if (myCombination.isProxyActive(i) && myCombination.isProxySignificant(i)) pos++; + pos -= 1; + + xMin = getZmin(subsetInterpolationPoints); + xMax = getZmax(subsetInterpolationPoints); + + if (interpolationSettings->getUseMultipleDetrending() && pos < parameters.size()) + { + std::vector xVector; + for (int m = xMin; m < xMax; m += 5) + xVector.push_back(m); + + for (int p = 0; p < int(xVector.size()); p++) + { + point.setX(xVector[p]); + point.setY(functionLinear_intercept(xVector[p], parameters[pos])); + point_vector.append(point); + } + + if (interpolationSettings->getProxy(proxyPos)->getRegressionR2() != NODATA) + r2.setText(QString("%1").arg(interpolationSettings->getProxy(proxyPos)->getRegressionR2(), 0, 'f', 3)); + + if (parameters[pos].size() == 2) + { + par0.setText(QString("%1").arg(parameters[pos][0], 0, 'f', 4)); + par1.setText(QString("%1").arg(parameters[pos][1], 0, 'f', 4)); + } + }*/ + + } - } chartView->drawModelLapseRate(point_vector); } } From 086354cc33e8e43be724a050cf952a2c51ff947a Mon Sep 17 00:00:00 2001 From: cate-to Date: Wed, 30 Oct 2024 11:42:30 +0100 Subject: [PATCH 12/34] cleanup glocal detremnding --- interpolation/interpolation.cpp | 39 +--- interpolation/interpolationSettings.cpp | 32 +-- interpolation/interpolationSettings.h | 13 +- project/project.cpp | 259 ++++++++++++------------ 4 files changed, 148 insertions(+), 195 deletions(-) diff --git a/interpolation/interpolation.cpp b/interpolation/interpolation.cpp index 21fc5dd4a..5fa3dd33c 100644 --- a/interpolation/interpolation.cpp +++ b/interpolation/interpolation.cpp @@ -2076,13 +2076,10 @@ void detrendingOtherProxies(int elevationPos, std::vector &myPoints, Crit3DInterpolationSettings* mySettings, meteoVariable myVar, std::string& errorStr) { - //per n aree, si creano gruppi di stazioni in base al codice (sorta di selection) - //poi si chiama fitting per ogni area, il fitting va salvato in una matrice - //matrice avrà, per ogni area, i parametri di fitting + la combination + il valore delle stazioni detrendate secondo il fitting corrispondente std::vector subsetPoints; - int i = 0; std::vector macroAreas = mySettings->getMacroAreas(); int areasNr = macroAreas.size(); + int i = 0; int elevationPos = NODATA; for (unsigned int pos=0; pos < mySettings->getSelectedCombination().getProxySize(); pos++) @@ -2091,14 +2088,10 @@ bool glocalDetrendingFitting(std::vector &myPoint elevationPos = pos; } - std::vector>> allAreaParameters; - std::vector allAreaCombinations; - //now create the subset of points, including only the right points, and then call multipleDetrending - for (i = 0; i < areasNr; i++) //ATTENZIONE INDICI + //create the subset of points starting from the meteopoints vector saved in every macro area + for (i = 0; i < areasNr; i++) //ATTENZIONE INDICI ? { - Crit3DMacroArea myArea = macroAreas[i]; - - std::vector temp = myArea.getMeteoPoints(); + std::vector temp = macroAreas[i].getMeteoPoints(); if (! temp.empty()) { std::vector subsetPoints; @@ -2116,35 +2109,23 @@ bool glocalDetrendingFitting(std::vector &myPoint mySettings->setCurrentCombination(mySettings->getSelectedCombination()); mySettings->clearFitting(); + //fitting and elevation detrending (necessary for the fitting of other proxies) if (elevationPos != NODATA && mySettings->getSelectedCombination().isProxyActive(elevationPos)) { if (!multipleDetrendingElevationFitting(elevationPos, subsetPoints, mySettings, myVar, errorStr, false)) return false; - detrendingElevation(elevationPos, subsetPoints, mySettings); - } if (!multipleDetrendingOtherProxiesFitting(elevationPos, subsetPoints, mySettings, myVar, errorStr)) return false; - detrendingOtherProxies(elevationPos, subsetPoints, mySettings); - - // - allAreaParameters.push_back(mySettings->getFittingParameters()); - allAreaCombinations.push_back(mySettings->getCurrentCombination()); - } - else - { - std::vector> emptyParameters; - Crit3DProxyCombination emptyCombination; - allAreaParameters.push_back(emptyParameters); - allAreaCombinations.push_back(emptyCombination); + //save parameters and combination in the macro area + macroAreas[i].setParameters(mySettings->getFittingParameters()); + macroAreas[i].setCombination(mySettings->getCurrentCombination()); } subsetPoints.clear(); } - - mySettings->setAreaParameters(allAreaParameters); - mySettings->setAreaCombination(allAreaCombinations); - + //update macro areas in settings + mySettings->setMacroAreas(macroAreas); return true; } diff --git a/interpolation/interpolationSettings.cpp b/interpolation/interpolationSettings.cpp index f286df8d0..8618b983b 100644 --- a/interpolation/interpolationSettings.cpp +++ b/interpolation/interpolationSettings.cpp @@ -462,8 +462,6 @@ void Crit3DInterpolationSettings::initialize() fittingFunction.clear(); fittingParameters.clear(); macroAreas.clear(); - areaCombination.clear(); - areaParameters.clear(); isKrigingReady = false; precipitationAllZero = false; @@ -519,36 +517,10 @@ gis::Crit3DRasterGrid* Crit3DInterpolationSettings::getMacroAreasMap() return macroAreasMap; } -std::vector>> Crit3DInterpolationSettings::getAreaParameters() -{ - return areaParameters; -} - -void Crit3DInterpolationSettings::setAreaParameters(std::vector > > myAreaParameters) -{ - areaParameters.clear(); - areaParameters = myAreaParameters; -} - -std::vector Crit3DInterpolationSettings::getAreaCombination() -{ - return areaCombination; -} - -void Crit3DInterpolationSettings::setAreaCombination(std::vector myAreaCombination) -{ - areaCombination.clear(); - areaCombination = myAreaCombination; -} - std::vector Crit3DInterpolationSettings::getMacroAreas() { return macroAreas; } -void Crit3DInterpolationSettings::addMacroArea(Crit3DMacroArea myArea) -{ - macroAreas.push_back(myArea); -} void Crit3DInterpolationSettings::setMacroAreas(std::vector myAreas) { @@ -1050,13 +1022,13 @@ std::vector Crit3DMacroArea::getAreaCells() return areaCells; } -void Crit3DMacroArea::setParameters (std::vector> myParameters) +void Crit3DMacroArea::setParameters (std::vector> myParameters) { areaParameters = myParameters; return; } -std::vector> Crit3DMacroArea::getParameters() +std::vector> Crit3DMacroArea::getParameters() { return areaParameters; } diff --git a/interpolation/interpolationSettings.h b/interpolation/interpolationSettings.h index 4371d82bb..6f8072af0 100644 --- a/interpolation/interpolationSettings.h +++ b/interpolation/interpolationSettings.h @@ -136,7 +136,7 @@ { private: Crit3DProxyCombination areaCombination; - std::vector> areaParameters; + std::vector> areaParameters; std::vector meteoPoints; std::vector areaCells; @@ -145,8 +145,8 @@ void setMeteoPoints (std::vector myMeteoPoints); std::vector getMeteoPoints(); - void setParameters (std::vector> myParameters); - std::vector> getParameters(); + void setParameters (std::vector> myParameters); + std::vector> getParameters(); void setCombination (Crit3DProxyCombination myCombination); Crit3DProxyCombination getCombination(); void setAreaCells (std::vector myCells); @@ -200,8 +200,6 @@ std::vector > fittingParameters; std::vector&)>> fittingFunction; std::vector pointsRange; - std::vector>> areaParameters; - std::vector areaCombination; std::vector macroAreas; @@ -236,12 +234,7 @@ bool getUseGlocalDetrending(); void setMacroAreasMap(gis::Crit3DRasterGrid *value); gis::Crit3DRasterGrid *getMacroAreasMap(); - std::vector > > getAreaParameters(); - void setAreaParameters(std::vector>> myAreaParameters); - std::vector getAreaCombination(); - void setAreaCombination(std::vector myAreaCombination); std::vector getMacroAreas(); - void addMacroArea(Crit3DMacroArea myArea); void setMacroAreas(std::vector myAreas); void setUseDoNotRetrend(bool myValue); diff --git a/project/project.cpp b/project/project.cpp index 894665dfe..830371cbc 100644 --- a/project/project.cpp +++ b/project/project.cpp @@ -2860,7 +2860,7 @@ bool Project::interpolationDemGlocalDetrending(int numAreas, meteoVariable myVar Crit3DProxyCombination myCombination = interpolationSettings.getSelectedCombination(); interpolationSettings.setCurrentCombination(myCombination); - //obtain fitting parameters for every macro area + //obtain fitting parameters and combination for every macro area if (! preInterpolation(interpolationPoints, &interpolationSettings, meteoSettings, &climateParameters, meteoPoints, nrMeteoPoints, myVar, myTime, errorStdStr)) { @@ -2885,8 +2885,7 @@ bool Project::interpolationDemGlocalDetrending(int numAreas, meteoVariable myVar //preparing to start a "global" interpolation for every single macro area unsigned row, col; - std::vector>> allAreaFittingParameters = interpolationSettings.getAreaParameters(); - std::vector allAreaCombinations = interpolationSettings.getAreaCombination(); + float z; std::vector areaCells; int elevationPos = NODATA; @@ -2896,72 +2895,76 @@ bool Project::interpolationDemGlocalDetrending(int numAreas, meteoVariable myVar elevationPos = pos; } - if (allAreaFittingParameters.size() < numAreas || allAreaCombinations.size() < numAreas) + if (interpolationSettings.getMacroAreas().size() != numAreas) { errorString = "Error in function interpolationGrid: there is something wrong with the macro areas raster file."; return false; } - for (unsigned areaIndex = 0; areaIndex < numAreas; areaIndex++) + for (unsigned areaIndex = 4; areaIndex < 5; areaIndex++) { - //save all the cells in that specific area in a vector (areaCells) + //load macro area and its cells Crit3DMacroArea myArea = interpolationSettings.getMacroAreas()[areaIndex]; areaCells = myArea.getAreaCells(); - //groupCellsInArea(areaCells, areaIndex, false); - - //take the parameters+combination for that area - std::vector &)>> fittingFunction; - interpolationSettings.setFittingParameters(allAreaFittingParameters[areaIndex]); - interpolationSettings.setCurrentCombination(allAreaCombinations[areaIndex]); - //find the fitting functions vector based on the length of the parameters vector for every proxy - for (int l = 0; l < allAreaFittingParameters[areaIndex].size(); l++) + if (!areaCells.empty()) { - if (allAreaFittingParameters[areaIndex][l].size() == 2) - fittingFunction.push_back(functionLinear_intercept); - else if (allAreaFittingParameters[areaIndex][l].size() == 4) - fittingFunction.push_back(lapseRatePiecewise_two); - else if (allAreaFittingParameters[areaIndex][l].size() == 5) - fittingFunction.push_back(lapseRatePiecewise_three); - else if (allAreaFittingParameters[areaIndex][l].size() == 6) - fittingFunction.push_back(lapseRatePiecewise_three_free); - } + //take the parameters+combination for that area + interpolationSettings.setFittingParameters(myArea.getParameters()); + interpolationSettings.setCurrentCombination(myArea.getCombination()); - //create group of macro area interpolation points - interpolationSettings.setFittingFunction(fittingFunction); - std::vector temp = myArea.getMeteoPoints(); - std::vector subsetInterpolationPoints; - for (int l = 0; l < temp.size(); l++) - { - for (int k = 0; k < interpolationPoints.size(); k++) - if (interpolationPoints[k].index == temp[l]) - subsetInterpolationPoints.push_back(interpolationPoints[k]); - } + //find the fitting functions vector based on the length of the parameters vector for every proxy + std::vector &)>> fittingFunction; - //detrending (done once for every macro area) - if (elevationPos != NODATA && myCombination.isProxyActive(elevationPos)) - detrendingElevation(elevationPos, subsetInterpolationPoints, &interpolationSettings); + for (int l = 0; l < myArea.getParameters().size(); l++) + { + if (myArea.getParameters()[l].size() == 2) + fittingFunction.push_back(functionLinear_intercept); + else if (myArea.getParameters()[l].size() == 4) + fittingFunction.push_back(lapseRatePiecewise_two); + else if (myArea.getParameters()[l].size() == 5) + fittingFunction.push_back(lapseRatePiecewise_three); + else if (myArea.getParameters()[l].size() == 6) + fittingFunction.push_back(lapseRatePiecewise_three_free); + } + interpolationSettings.setFittingFunction(fittingFunction); - detrendingOtherProxies(elevationPos, subsetInterpolationPoints, &interpolationSettings); + //create group of macro area interpolation points + std::vector temp = myArea.getMeteoPoints(); + std::vector subsetInterpolationPoints; + for (int l = 0; l < temp.size(); l++) + { + for (int k = 0; k < interpolationPoints.size(); k++) + if (interpolationPoints[k].index == temp[l]) + subsetInterpolationPoints.push_back(interpolationPoints[k]); + } - for (unsigned cellIndex = 0; cellIndex < areaCells.size(); cellIndex = cellIndex + 2) - { - row = unsigned(areaCells[cellIndex]/DEM.header->nrCols); - col = int(areaCells[cellIndex])%DEM.header->nrCols; + //detrending (done once for every macro area) + if (elevationPos != NODATA && myCombination.isProxyActive(elevationPos)) + detrendingElevation(elevationPos, subsetInterpolationPoints, &interpolationSettings); - float z = DEM.value[row][col]; - if (! isEqual(z, myHeader.flag)) + detrendingOtherProxies(elevationPos, subsetInterpolationPoints, &interpolationSettings); + + //calculate value for every cell + for (unsigned cellIndex = 0; cellIndex < areaCells.size(); cellIndex = cellIndex + 2) { - gis::getUtmXYFromRowCol(myHeader, row, col, &x, &y); + row = unsigned(areaCells[cellIndex]/DEM.header->nrCols); + col = int(areaCells[cellIndex])%DEM.header->nrCols; + z = DEM.value[row][col]; - getProxyValuesXY(x, y, &interpolationSettings, proxyValues); + if (! isEqual(z, myHeader.flag)) + { + gis::getUtmXYFromRowCol(myHeader, row, col, &x, &y); - if (isEqual(myRaster->value[row][col], NODATA)) - myRaster->value[row][col] = interpolate(subsetInterpolationPoints, &interpolationSettings, meteoSettings, - myVar, x, y, z, proxyValues, true)*areaCells[cellIndex+1]; - else - myRaster->value[row][col] += interpolate(subsetInterpolationPoints, &interpolationSettings, meteoSettings, - myVar, x, y, z, proxyValues, true)*areaCells[cellIndex+1]; + getProxyValuesXY(x, y, &interpolationSettings, proxyValues); + + if (isEqual(myRaster->value[row][col], NODATA)) + myRaster->value[row][col] = interpolate(subsetInterpolationPoints, &interpolationSettings, meteoSettings, + myVar, x, y, z, proxyValues, true)*areaCells[cellIndex+1]; + else + myRaster->value[row][col] += interpolate(subsetInterpolationPoints, &interpolationSettings, meteoSettings, + myVar, x, y, z, proxyValues, true)*areaCells[cellIndex+1]; + } } } } @@ -3317,12 +3320,10 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) else { unsigned row, col; - std::vector>> allAreaFittingParameters = interpolationSettings.getAreaParameters(); - std::vector allAreaCombinations = interpolationSettings.getAreaCombination(); std::vector areaCells; float myValue = NODATA; - if (allAreaCombinations.size() != numAreas) + if (interpolationSettings.getMacroAreas().size() != numAreas) { errorString = "Error in function interpolationGrid: there is something wrong with the macro areas raster file."; return false; @@ -3337,106 +3338,112 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) for (unsigned areaIndex = 0; areaIndex < numAreas; areaIndex++) { + //load macro area and its cells Crit3DMacroArea myArea = interpolationSettings.getMacroAreas()[areaIndex]; areaCells = myArea.getAreaCells(); - //groupCellsInArea(areaCells, areaIndex, true); - std::vector &)> > fittingFunction; - interpolationSettings.setFittingParameters(allAreaFittingParameters[areaIndex]); - interpolationSettings.setCurrentCombination(allAreaCombinations[areaIndex]); - - for (int l = 0; l < allAreaFittingParameters[areaIndex].size(); l++) + if (!areaCells.empty()) { - if (allAreaFittingParameters[areaIndex][l].size() == 2) - fittingFunction.push_back(functionLinear_intercept); - else if (allAreaFittingParameters[areaIndex][l].size() == 4) - fittingFunction.push_back(lapseRatePiecewise_two); - else if (allAreaFittingParameters[areaIndex][l].size() == 5) - fittingFunction.push_back(lapseRatePiecewise_three); - else if (allAreaFittingParameters[areaIndex][l].size() == 6) - fittingFunction.push_back(lapseRatePiecewise_three_free); - } + //take the parameters+combination for that area + interpolationSettings.setFittingParameters(myArea.getParameters()); + interpolationSettings.setCurrentCombination(myArea.getCombination()); - interpolationSettings.setFittingFunction(fittingFunction); - std::vector temp = myArea.getMeteoPoints(); - std::vector subsetInterpolationPoints; - for (int l = 0; l < temp.size(); l++) - { - for (int k = 0; k < interpolationPoints.size(); k++) - if (interpolationPoints[k].index == temp[l]) - subsetInterpolationPoints.push_back(interpolationPoints[k]); - } - - if (elevationPos != NODATA && myCombination.isProxyActive(elevationPos)) - detrendingElevation(elevationPos, subsetInterpolationPoints, &interpolationSettings); - - detrendingOtherProxies(elevationPos, subsetInterpolationPoints, &interpolationSettings); + //find the fitting functions vector based on the length of the parameters vector for every proxy + std::vector &)> > fittingFunction; + for (int l = 0; l < myArea.getParameters().size(); l++) + { + if (myArea.getParameters()[l].size() == 2) + fittingFunction.push_back(functionLinear_intercept); + else if (myArea.getParameters()[l].size() == 4) + fittingFunction.push_back(lapseRatePiecewise_two); + else if (myArea.getParameters()[l].size() == 5) + fittingFunction.push_back(lapseRatePiecewise_three); + else if (myArea.getParameters()[l].size() == 6) + fittingFunction.push_back(lapseRatePiecewise_three_free); + } + interpolationSettings.setFittingFunction(fittingFunction); - unsigned nrCols = meteoGridDbHandler->meteoGrid()->gridStructure().header().nrCols; + //create vector of macro area interpolation points + std::vector temp = myArea.getMeteoPoints(); + std::vector subsetInterpolationPoints; + for (int l = 0; l < temp.size(); l++) + { + for (int k = 0; k < interpolationPoints.size(); k++) + if (interpolationPoints[k].index == temp[l]) + subsetInterpolationPoints.push_back(interpolationPoints[k]); + } - for (unsigned cellIndex = 0; cellIndex < areaCells.size(); cellIndex = cellIndex + 2) - { - row = unsigned(areaCells[cellIndex]/nrCols); - col = int(areaCells[cellIndex])%nrCols; + //detrending + if (elevationPos != NODATA && myCombination.isProxyActive(elevationPos)) + detrendingElevation(elevationPos, subsetInterpolationPoints, &interpolationSettings); + detrendingOtherProxies(elevationPos, subsetInterpolationPoints, &interpolationSettings); - if (meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->active) + unsigned nrCols = meteoGridDbHandler->meteoGrid()->gridStructure().header().nrCols; + //calculate value for every cell + for (unsigned cellIndex = 0; cellIndex < areaCells.size(); cellIndex = cellIndex + 2) { - myX = meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->point.utm.x; - myY = meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->point.utm.y; - myZ = meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->point.z; + row = unsigned(areaCells[cellIndex]/nrCols); + col = int(areaCells[cellIndex])%nrCols; - if (getUseDetrendingVar(myVar)) + if (meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->active) { - proxyIndex = 0; + myX = meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->point.utm.x; + myY = meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->point.utm.y; + myZ = meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->point.z; - for (i=0; i < interpolationSettings.getProxyNr(); i++) + if (getUseDetrendingVar(myVar)) { - proxyValues[i] = NODATA; + proxyIndex = 0; - if (myCombination.isProxyActive(i)) + for (i=0; i < interpolationSettings.getProxyNr(); i++) { - if (proxyIndex < meteoGridProxies.size()) + proxyValues[i] = NODATA; + + if (interpolationSettings.getCurrentCombination().isProxyActive(i)) { - float proxyValue = gis::getValueFromXY(*meteoGridProxies[proxyIndex], myX, myY); - if (proxyValue != meteoGridProxies[proxyIndex]->header->flag) - proxyValues[i] = double(proxyValue); + if (proxyIndex < meteoGridProxies.size()) + { + float proxyValue = gis::getValueFromXY(*meteoGridProxies[proxyIndex], myX, myY); + if (proxyValue != meteoGridProxies[proxyIndex]->header->flag) + proxyValues[i] = double(proxyValue); + } + + proxyIndex++; } - - proxyIndex++; } - } - interpolatedValue = interpolate(subsetInterpolationPoints, &interpolationSettings, meteoSettings, myVar, myX, myY, myZ, proxyValues, true) - * areaCells[cellIndex + 1]; + interpolatedValue = interpolate(subsetInterpolationPoints, &interpolationSettings, meteoSettings, myVar, myX, myY, myZ, proxyValues, true) + * areaCells[cellIndex + 1]; - } - if (freq == hourly) - { - if (meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->nrObsDataDaysH == 0) - meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->initializeObsDataH(1, 1, myTime.date); + } + if (freq == hourly) + { + if (meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->nrObsDataDaysH == 0) + meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->initializeObsDataH(1, 1, myTime.date); - myValue = meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->getMeteoPointValueH(myTime.date, myTime.getHour(), myTime.getMinutes(), myVar); + myValue = meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->getMeteoPointValueH(myTime.date, myTime.getHour(), myTime.getMinutes(), myVar); - if (isEqual(myValue, NODATA)) - myValue = 0; + if (isEqual(myValue, NODATA)) + myValue = 0; - meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->setMeteoPointValueH(myTime.date, myTime.getHour(), myTime.getMinutes(), myVar, float(interpolatedValue+myValue)); - meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->currentValue = float(interpolatedValue+myValue); + meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->setMeteoPointValueH(myTime.date, myTime.getHour(), myTime.getMinutes(), myVar, float(interpolatedValue+myValue)); + meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->currentValue = float(interpolatedValue+myValue); - } - else if (freq == daily) - { - if (meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->nrObsDataDaysD == 0) - meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->initializeObsDataD(1, myTime.date); + } + else if (freq == daily) + { + if (meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->nrObsDataDaysD == 0) + meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->initializeObsDataD(1, myTime.date); - myValue = meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->getMeteoPointValueD(myTime.date, myVar); + myValue = meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->getMeteoPointValueD(myTime.date, myVar); - if (isEqual(myValue, NODATA)) - myValue = 0; + if (isEqual(myValue, NODATA)) + myValue = 0; - meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->setMeteoPointValueD(myTime.date, myVar, float(interpolatedValue+myValue)); - meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->currentValue = float(interpolatedValue+myValue); + meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->setMeteoPointValueD(myTime.date, myVar, float(interpolatedValue+myValue)); + meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->currentValue = float(interpolatedValue+myValue); + } } } } From ca405274a19972b5edf4a05a5b9778e7ffe0d8a5 Mon Sep 17 00:00:00 2001 From: cate-to Date: Wed, 30 Oct 2024 12:02:40 +0100 Subject: [PATCH 13/34] debug --- project/project.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/project/project.cpp b/project/project.cpp index 830371cbc..1b1d38fd1 100644 --- a/project/project.cpp +++ b/project/project.cpp @@ -2901,7 +2901,7 @@ bool Project::interpolationDemGlocalDetrending(int numAreas, meteoVariable myVar return false; } - for (unsigned areaIndex = 4; areaIndex < 5; areaIndex++) + for (unsigned areaIndex = 0; areaIndex < numAreas; areaIndex++) { //load macro area and its cells Crit3DMacroArea myArea = interpolationSettings.getMacroAreas()[areaIndex]; From 50672543d149533f11dab2a3faf0b0c0c85ed17f Mon Sep 17 00:00:00 2001 From: cate-to Date: Wed, 30 Oct 2024 12:56:17 +0100 Subject: [PATCH 14/34] fix grid glocal detrending --- project/project.cpp | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/project/project.cpp b/project/project.cpp index 1b1d38fd1..7f49c30f6 100644 --- a/project/project.cpp +++ b/project/project.cpp @@ -3319,6 +3319,20 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) } else { + //set all cells to NODATA + for (unsigned col = 0; col < unsigned(meteoGridDbHandler->meteoGrid()->gridStructure().header().nrCols); col++) + { + for (unsigned row = 0; row < unsigned(meteoGridDbHandler->meteoGrid()->gridStructure().header().nrRows); row++) + { + if (meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->active) + { + if (freq == hourly) + meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->currentValue = NODATA; + else if (freq == daily) + meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->currentValue = NODATA; + } + } + } unsigned row, col; std::vector areaCells; float myValue = NODATA; @@ -3414,14 +3428,13 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) } interpolatedValue = interpolate(subsetInterpolationPoints, &interpolationSettings, meteoSettings, myVar, myX, myY, myZ, proxyValues, true) * areaCells[cellIndex + 1]; - } if (freq == hourly) { if (meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->nrObsDataDaysH == 0) meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->initializeObsDataH(1, 1, myTime.date); - myValue = meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->getMeteoPointValueH(myTime.date, myTime.getHour(), myTime.getMinutes(), myVar); + myValue = meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->currentValue; if (isEqual(myValue, NODATA)) myValue = 0; @@ -3435,7 +3448,7 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) if (meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->nrObsDataDaysD == 0) meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->initializeObsDataD(1, myTime.date); - myValue = meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->getMeteoPointValueD(myTime.date, myVar); + myValue = meteoGridDbHandler->meteoGrid()->meteoPoints()[row][col]->currentValue; if (isEqual(myValue, NODATA)) myValue = 0; From 17b366aeba37756f20534c94cf22af71e921bef1 Mon Sep 17 00:00:00 2001 From: cate-to Date: Thu, 31 Oct 2024 11:39:31 +0100 Subject: [PATCH 15/34] glocal detrending period --- pragaProject/pragaProject.cpp | 7 +++++++ project/project.cpp | 37 +++++------------------------------ project/project.h | 2 +- 3 files changed, 13 insertions(+), 33 deletions(-) diff --git a/pragaProject/pragaProject.cpp b/pragaProject/pragaProject.cpp index 27b4aeecc..76c426a39 100644 --- a/pragaProject/pragaProject.cpp +++ b/pragaProject/pragaProject.cpp @@ -2635,6 +2635,13 @@ bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QL if (! loadTopographicDistanceMaps(true, false)) return false; } + if (interpolationSettings.getUseGlocalDetrending()) + { + if (loadGlocalAreasMap()) + if (loadGlocalStationsAndCells(!interpolationSettings.getMeteoGridUpscaleFromDem()) < 1) return false; + else + return false; + } // save also time aggregated variables foreach (myVar, aggrVariables) diff --git a/project/project.cpp b/project/project.cpp index 7f49c30f6..cf6ff2329 100644 --- a/project/project.cpp +++ b/project/project.cpp @@ -2835,7 +2835,7 @@ bool Project::interpolationDemLocalDetrending(meteoVariable myVar, const Crit3DT return true; } -bool Project::interpolationDemGlocalDetrending(int numAreas, meteoVariable myVar, const Crit3DTime& myTime, gis::Crit3DRasterGrid *myRaster) +bool Project::interpolationDemGlocalDetrending(meteoVariable myVar, const Crit3DTime& myTime, gis::Crit3DRasterGrid *myRaster) { if (!getUseDetrendingVar(myVar) || !interpolationSettings.getUseGlocalDetrending()) return false; @@ -2895,13 +2895,7 @@ bool Project::interpolationDemGlocalDetrending(int numAreas, meteoVariable myVar elevationPos = pos; } - if (interpolationSettings.getMacroAreas().size() != numAreas) - { - errorString = "Error in function interpolationGrid: there is something wrong with the macro areas raster file."; - return false; - } - - for (unsigned areaIndex = 0; areaIndex < numAreas; areaIndex++) + for (unsigned areaIndex = 0; areaIndex < interpolationSettings.getMacroAreas().size(); areaIndex++) { //load macro area and its cells Crit3DMacroArea myArea = interpolationSettings.getMacroAreas()[areaIndex]; @@ -3119,17 +3113,6 @@ bool Project::interpolationDemMain(meteoVariable myVar, const Crit3DTime& myTime if (interpolationSettings.getUseMultipleDetrending()) interpolationSettings.clearFitting(); - int numAreas = 1; - if (interpolationSettings.getUseGlocalDetrending()) - { - if (loadGlocalAreasMap()) - { - numAreas = loadGlocalStationsAndCells(false); - if (numAreas < 1) return false; - } - else - return false; - } // dynamic lapserate if (getUseDetrendingVar(myVar) && interpolationSettings.getUseLocalDetrending()) { @@ -3137,7 +3120,7 @@ bool Project::interpolationDemMain(meteoVariable myVar, const Crit3DTime& myTime } else if (getUseDetrendingVar(myVar) && interpolationSettings.getUseGlocalDetrending()) { - return interpolationDemGlocalDetrending(numAreas, myVar, myTime, myRaster); + return interpolationDemGlocalDetrending(myVar, myTime, myRaster); } else { @@ -3184,14 +3167,10 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) if (interpolationSettings.getUseMultipleDetrending()) interpolationSettings.clearFitting(); - unsigned numAreas = 1; if (interpolationSettings.getUseGlocalDetrending()) { if (loadGlocalAreasMap()) - { - numAreas = loadGlocalStationsAndCells(true); - if (numAreas < 1) return false; - } + if (loadGlocalStationsAndCells(true) < 1) return false; else return false; } @@ -3337,12 +3316,6 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) std::vector areaCells; float myValue = NODATA; - if (interpolationSettings.getMacroAreas().size() != numAreas) - { - errorString = "Error in function interpolationGrid: there is something wrong with the macro areas raster file."; - return false; - } - int elevationPos = NODATA; for (unsigned int pos=0; pos < interpolationSettings.getCurrentCombination().getProxySize(); pos++) { @@ -3350,7 +3323,7 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) elevationPos = pos; } - for (unsigned areaIndex = 0; areaIndex < numAreas; areaIndex++) + for (unsigned areaIndex = 0; areaIndex < interpolationSettings.getMacroAreas().size(); areaIndex++) { //load macro area and its cells Crit3DMacroArea myArea = interpolationSettings.getMacroAreas()[areaIndex]; diff --git a/project/project.h b/project/project.h index 46302b3d1..43bc937b2 100644 --- a/project/project.h +++ b/project/project.h @@ -273,7 +273,7 @@ bool interpolationDemMain(meteoVariable myVar, const Crit3DTime& myTime, gis::Crit3DRasterGrid *myRaster); bool interpolationDem(meteoVariable myVar, const Crit3DTime& myTime, gis::Crit3DRasterGrid *myRaster); bool interpolationDemLocalDetrending(meteoVariable myVar, const Crit3DTime& myTime, gis::Crit3DRasterGrid *myRaster); - bool interpolationDemGlocalDetrending(int numAreas, meteoVariable myVar, const Crit3DTime& myTime, gis::Crit3DRasterGrid *myRaster); + bool interpolationDemGlocalDetrending(meteoVariable myVar, const Crit3DTime& myTime, gis::Crit3DRasterGrid *myRaster); bool interpolateDemRadiation(const Crit3DTime& myTime, gis::Crit3DRasterGrid *myRaster); bool interpolationOutputPoints(std::vector &interpolationPoints, gis::Crit3DRasterGrid *outputGrid, meteoVariable myVar); From b8648b0dd61d495962432727dc1d5284887f6033 Mon Sep 17 00:00:00 2001 From: Gabriele Antolini Date: Thu, 31 Oct 2024 11:45:51 +0100 Subject: [PATCH 16/34] td for glocal --- interpolation/interpolation.cpp | 2 +- project/dialogInterpolation.cpp | 3 --- 2 files changed, 1 insertion(+), 4 deletions(-) diff --git a/interpolation/interpolation.cpp b/interpolation/interpolation.cpp index 5fa3dd33c..60cf172a5 100644 --- a/interpolation/interpolation.cpp +++ b/interpolation/interpolation.cpp @@ -2256,7 +2256,7 @@ bool preInterpolation(std::vector &myPoints, Crit } } - if (mySettings->getUseTD() && getUseTdVar(myVar)) + if (mySettings->getUseTD() && getUseTdVar(myVar) && ! mySettings->getUseGlocalDetrending()) { topographicDistanceOptimize(myVar, myMeteoPoints, nrMeteoPoints, myPoints, mySettings, meteoSettings); } diff --git a/project/dialogInterpolation.cpp b/project/dialogInterpolation.cpp index 8b4151a93..8e799db08 100644 --- a/project/dialogInterpolation.cpp +++ b/project/dialogInterpolation.cpp @@ -239,9 +239,6 @@ void DialogInterpolation::multipleDetrendingChanged(int active) void DialogInterpolation::localDetrendingChanged(int active) { - if (active == Qt::Checked) topographicDistanceEdit->setChecked(Qt::Unchecked); - topographicDistanceEdit->setEnabled(active == Qt::Unchecked); - maxTdMultiplierEdit.setEnabled(active == Qt::Unchecked); minPointsLocalDetrendingEdit.setEnabled(active == Qt::Checked); if (active == Qt::Checked) optimalDetrendingEdit->setChecked(Qt::Unchecked); optimalDetrendingEdit->setEnabled(active == Qt::Unchecked); From fe7e973a3c7bea7f8e6a60cbe827e619195113de Mon Sep 17 00:00:00 2001 From: Gabriele Antolini Date: Thu, 31 Oct 2024 18:05:19 +0100 Subject: [PATCH 17/34] interpolation grid period without upscaling 1st part --- pragaProject/pragaProject.cpp | 46 +++++++++++++++++++++++++---------- project/project.cpp | 44 ++++++++++++++++++++++++++++----- project/project.h | 3 ++- 3 files changed, 73 insertions(+), 20 deletions(-) diff --git a/pragaProject/pragaProject.cpp b/pragaProject/pragaProject.cpp index 76c426a39..9c379fd63 100644 --- a/pragaProject/pragaProject.cpp +++ b/pragaProject/pragaProject.cpp @@ -2578,8 +2578,11 @@ bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QL QList varToSave; int countDaysSaving = 0; - if (pragaDailyMaps == nullptr) pragaDailyMaps = new Crit3DDailyMeteoMaps(DEM); - if (pragaHourlyMaps == nullptr) pragaHourlyMaps = new PragaHourlyMeteoMaps(DEM); + if (interpolationSettings.getMeteoGridUpscaleFromDem()) + { + if (pragaDailyMaps == nullptr) pragaDailyMaps = new Crit3DDailyMeteoMaps(DEM); + if (pragaHourlyMaps == nullptr) pragaHourlyMaps = new PragaHourlyMeteoMaps(DEM); + } // find out needed frequency foreach (myVar, variables) @@ -2635,12 +2638,11 @@ bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QL if (! loadTopographicDistanceMaps(true, false)) return false; } + if (interpolationSettings.getUseGlocalDetrending()) { - if (loadGlocalAreasMap()) - if (loadGlocalStationsAndCells(!interpolationSettings.getMeteoGridUpscaleFromDem()) < 1) return false; - else - return false; + if (! loadGlocalAreasMap()) return false; + if (! loadGlocalStationsAndCells(!interpolationSettings.getMeteoGridUpscaleFromDem())) return false; } // save also time aggregated variables @@ -2712,11 +2714,21 @@ bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QL if (myVar == airRelHumidity && interpolationSettings.getUseDewPoint()) { if (interpolationSettings.getUseInterpolatedTForRH()) - passInterpolatedTemperatureToHumidityPoints(getCrit3DTime(myDate, myHour), meteoSettings); - if (! interpolationDemMain(airDewTemperature, getCrit3DTime(myDate, myHour), hourlyMeteoMaps->mapHourlyTdew)) return false; - hourlyMeteoMaps->computeRelativeHumidityMap(hourlyMeteoMaps->mapHourlyRelHum); + { + if (interpolationSettings.getMeteoGridUpscaleFromDem()) + { + passInterpolatedTemperatureToHumidityPoints(getCrit3DTime(myDate, myHour), meteoSettings); + if (! interpolationDemMain(airDewTemperature, getCrit3DTime(myDate, myHour), hourlyMeteoMaps->mapHourlyTdew)) return false; + hourlyMeteoMaps->computeRelativeHumidityMap(hourlyMeteoMaps->mapHourlyRelHum); + } + else + { + passGridTemperatureToHumidityPoints(getCrit3DTime(myDate, myHour), meteoSettings); + if (! interpolationMeteoGrid(airDewTemperature, hourly, getCrit3DTime(myDate, myHour))) return false; + } + } - if (saveRasters) + if (interpolationSettings.getMeteoGridUpscaleFromDem() && saveRasters) { myGrid = getPragaMapFromVar(airDewTemperature); rasterName = getMapFileOutName(airDewTemperature, myDate, myHour); @@ -2724,9 +2736,17 @@ bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QL } } else if (myVar == windVectorDirection || myVar == windVectorIntensity) { - if (! interpolationDemMain(windVectorX, getCrit3DTime(myDate, myHour), getPragaMapFromVar(windVectorX))) return false; - if (! interpolationDemMain(windVectorY, getCrit3DTime(myDate, myHour), getPragaMapFromVar(windVectorY))) return false; - if (! pragaHourlyMaps->computeWindVector()) return false; + if (interpolationSettings.getMeteoGridUpscaleFromDem()) + { + if (! interpolationDemMain(windVectorX, getCrit3DTime(myDate, myHour), getPragaMapFromVar(windVectorX))) return false; + if (! interpolationDemMain(windVectorY, getCrit3DTime(myDate, myHour), getPragaMapFromVar(windVectorY))) return false; + if (! pragaHourlyMaps->computeWindVector()) return false; + } + else + { + if (! interpolationMeteoGrid(windVectorX, hourly, getCrit3DTime(myDate, myHour))) return false; + if (! interpolationMeteoGrid(windVectorY, hourly, getCrit3DTime(myDate, myHour))) return false; + } } else if (myVar == leafWetness) { hourlyMeteoMaps->computeLeafWetnessMap() ; diff --git a/project/project.cpp b/project/project.cpp index cf6ff2329..6b62c9065 100644 --- a/project/project.cpp +++ b/project/project.cpp @@ -2272,7 +2272,7 @@ bool Project::loadGlocalAreasMap() return true; } -int Project::loadGlocalStationsAndCells(bool isGrid) +bool Project::loadGlocalStationsAndCells(bool isGrid) { //leggi csv aree QString mapsFolder = defaultPath + PATH_GEO; @@ -2316,7 +2316,7 @@ int Project::loadGlocalStationsAndCells(bool isGrid) interpolationSettings.setMacroAreas(myAreas); - return areaPoints.size(); + return (areaPoints.size() > 0); } bool Project::loadGlocalWeightMaps(std::vector &myAreas, bool isGrid) @@ -2519,6 +2519,40 @@ void Project::passInterpolatedTemperatureToHumidityPoints(Crit3DTime myTime, Cri } } +void Project::passGridTemperatureToHumidityPoints(Crit3DTime myTime, Crit3DMeteoSettings* meteoSettings) +{ + float airRelHum, airT; + int row, col; + double lat, lon; + + for (int i = 0; i < nrMeteoPoints; i++) + { + airRelHum = meteoPoints[i].getMeteoPointValue(myTime, airRelHumidity, meteoSettings); + airT = meteoPoints[i].getMeteoPointValue(myTime, airTemperature, meteoSettings); + + if (! isEqual(airRelHum, NODATA) && isEqual(airT, NODATA)) + { + if (meteoGridDbHandler->meteoGrid()->gridStructure().isUTM()) + { + gis::getRowColFromLonLat(meteoGridDbHandler->meteoGrid()->gridStructure().header(), + meteoPoints[i].point.utm.x, meteoPoints[i].point.utm.y, &row, &col); + } + else + { + gis::getLatLonFromUtm(gisSettings, meteoPoints[i].point.utm.x, meteoPoints[i].point.utm.y, &lat, &lon); + gis::getRowColFromLonLat(meteoGridDbHandler->meteoGrid()->gridStructure().header(), + lon, lat, &row, &col); + } + + if (! gis::isOutOfGridRowCol(row, col, meteoGridDbHandler->meteoGrid()->gridStructure().header())) + { + meteoPoints[i].setMeteoPointValueH(myTime.date, myTime.getHour(), myTime.getMinutes(), airTemperature, + meteoGridDbHandler->meteoGrid()->meteoPoint(row,col).getMeteoPointValueH(myTime.date, myTime.getHour(), myTime.getMinutes(), airTemperature)); + + } + } + } +} bool Project::interpolationOutputPoints(std::vector &interpolationPoints, gis::Crit3DRasterGrid *outputGrid, meteoVariable myVar) @@ -3169,10 +3203,8 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) if (interpolationSettings.getUseGlocalDetrending()) { - if (loadGlocalAreasMap()) - if (loadGlocalStationsAndCells(true) < 1) return false; - else - return false; + if (! loadGlocalAreasMap()) return false; + if (! loadGlocalStationsAndCells(!interpolationSettings.getMeteoGridUpscaleFromDem())) return false; } // check quality and pass data to interpolation diff --git a/project/project.h b/project/project.h index 43bc937b2..42d2386bc 100644 --- a/project/project.h +++ b/project/project.h @@ -260,8 +260,9 @@ bool writeTopographicDistanceMap(int pointIndex, const gis::Crit3DRasterGrid& demMap, QString pathTd); bool loadTopographicDistanceMaps(bool onlyWithData, bool showInfo); void passInterpolatedTemperatureToHumidityPoints(Crit3DTime myTime, Crit3DMeteoSettings *meteoSettings); + void passGridTemperatureToHumidityPoints(Crit3DTime myTime, Crit3DMeteoSettings* meteoSettings); bool loadGlocalAreasMap(); - int loadGlocalStationsAndCells(bool isGrid); + bool loadGlocalStationsAndCells(bool isGrid); bool loadGlocalWeightMaps(std::vector &myAreas, bool isGrid); bool loadGlocalStationsCsv(QString fileName, std::vector> &areaPoints); bool groupCellsInArea(std::vector &areaPoints, unsigned index, bool isGrid); From 81d78255d38d965b4098ad6b8776f545c9d76bbf Mon Sep 17 00:00:00 2001 From: Caterina Toscano Date: Mon, 4 Nov 2024 12:50:34 +0100 Subject: [PATCH 18/34] td for glocal detrending --- interpolation/interpolation.h | 2 + project/project.cpp | 125 +++++++++++++++------------------- project/project.h | 2 + 3 files changed, 58 insertions(+), 71 deletions(-) diff --git a/interpolation/interpolation.h b/interpolation/interpolation.h index 99fdfcf57..f145ccb33 100644 --- a/interpolation/interpolation.h +++ b/interpolation/interpolation.h @@ -49,6 +49,8 @@ bool preInterpolation(std::vector &myPoints, Crit3DInterpolationSettings *mySettings, Crit3DMeteoSettings *meteoSettings, Crit3DClimateParameters* myClimate, Crit3DMeteoPoint *myMeteoPoints, int nrMeteoPoints, meteoVariable myVar, Crit3DTime myTime, std::string &errorStr); + void topographicDistanceOptimize(meteoVariable myVar, Crit3DMeteoPoint* &myMeteoPoints, int nrMeteoPoints, std::vector &interpolationPoints, Crit3DInterpolationSettings* mySettings, Crit3DMeteoSettings* meteoSettings); + bool krigingEstimateVariogram(float *myDist, float *mySemiVar,int sizeMyVar, int nrMyPoints,float myMaxDistance, double *mySill, double *myNugget, double *myRange, double *mySlope, TkrigingMode *myMode, int nrPointData); bool krigLinearPrep(double *mySlope, double *myNugget, int nrPointData); diff --git a/project/project.cpp b/project/project.cpp index 6b62c9065..8343eedba 100644 --- a/project/project.cpp +++ b/project/project.cpp @@ -2934,45 +2934,11 @@ bool Project::interpolationDemGlocalDetrending(meteoVariable myVar, const Crit3D //load macro area and its cells Crit3DMacroArea myArea = interpolationSettings.getMacroAreas()[areaIndex]; areaCells = myArea.getAreaCells(); + std::vector subsetInterpolationPoints; if (!areaCells.empty()) { - //take the parameters+combination for that area - interpolationSettings.setFittingParameters(myArea.getParameters()); - interpolationSettings.setCurrentCombination(myArea.getCombination()); - - //find the fitting functions vector based on the length of the parameters vector for every proxy - std::vector &)>> fittingFunction; - - for (int l = 0; l < myArea.getParameters().size(); l++) - { - if (myArea.getParameters()[l].size() == 2) - fittingFunction.push_back(functionLinear_intercept); - else if (myArea.getParameters()[l].size() == 4) - fittingFunction.push_back(lapseRatePiecewise_two); - else if (myArea.getParameters()[l].size() == 5) - fittingFunction.push_back(lapseRatePiecewise_three); - else if (myArea.getParameters()[l].size() == 6) - fittingFunction.push_back(lapseRatePiecewise_three_free); - } - interpolationSettings.setFittingFunction(fittingFunction); - - //create group of macro area interpolation points - std::vector temp = myArea.getMeteoPoints(); - std::vector subsetInterpolationPoints; - for (int l = 0; l < temp.size(); l++) - { - for (int k = 0; k < interpolationPoints.size(); k++) - if (interpolationPoints[k].index == temp[l]) - subsetInterpolationPoints.push_back(interpolationPoints[k]); - } - - //detrending (done once for every macro area) - if (elevationPos != NODATA && myCombination.isProxyActive(elevationPos)) - detrendingElevation(elevationPos, subsetInterpolationPoints, &interpolationSettings); - - detrendingOtherProxies(elevationPos, subsetInterpolationPoints, &interpolationSettings); - + macroAreaDetrending(myArea, myVar, interpolationPoints, subsetInterpolationPoints, elevationPos); //calculate value for every cell for (unsigned cellIndex = 0; cellIndex < areaCells.size(); cellIndex = cellIndex + 2) { @@ -3360,44 +3326,11 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) //load macro area and its cells Crit3DMacroArea myArea = interpolationSettings.getMacroAreas()[areaIndex]; areaCells = myArea.getAreaCells(); + std::vector subsetInterpolationPoints; if (!areaCells.empty()) { - //take the parameters+combination for that area - interpolationSettings.setFittingParameters(myArea.getParameters()); - interpolationSettings.setCurrentCombination(myArea.getCombination()); - - //find the fitting functions vector based on the length of the parameters vector for every proxy - std::vector &)> > fittingFunction; - for (int l = 0; l < myArea.getParameters().size(); l++) - { - if (myArea.getParameters()[l].size() == 2) - fittingFunction.push_back(functionLinear_intercept); - else if (myArea.getParameters()[l].size() == 4) - fittingFunction.push_back(lapseRatePiecewise_two); - else if (myArea.getParameters()[l].size() == 5) - fittingFunction.push_back(lapseRatePiecewise_three); - else if (myArea.getParameters()[l].size() == 6) - fittingFunction.push_back(lapseRatePiecewise_three_free); - } - interpolationSettings.setFittingFunction(fittingFunction); - - //create vector of macro area interpolation points - std::vector temp = myArea.getMeteoPoints(); - std::vector subsetInterpolationPoints; - for (int l = 0; l < temp.size(); l++) - { - for (int k = 0; k < interpolationPoints.size(); k++) - if (interpolationPoints[k].index == temp[l]) - subsetInterpolationPoints.push_back(interpolationPoints[k]); - } - - //detrending - if (elevationPos != NODATA && myCombination.isProxyActive(elevationPos)) - detrendingElevation(elevationPos, subsetInterpolationPoints, &interpolationSettings); - - detrendingOtherProxies(elevationPos, subsetInterpolationPoints, &interpolationSettings); - + macroAreaDetrending(myArea, myVar, interpolationPoints, subsetInterpolationPoints, elevationPos); unsigned nrCols = meteoGridDbHandler->meteoGrid()->gridStructure().header().nrCols; //calculate value for every cell for (unsigned cellIndex = 0; cellIndex < areaCells.size(); cellIndex = cellIndex + 2) @@ -3471,6 +3404,56 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) return true; } +void Project::macroAreaDetrending(Crit3DMacroArea myArea, meteoVariable myVar, std::vector interpolationPoints, std::vector &subsetInterpolationPoints, int elevationPos) +{ + //take the parameters+combination for that area + interpolationSettings.setFittingParameters(myArea.getParameters()); + interpolationSettings.setCurrentCombination(myArea.getCombination()); + + //find the fitting functions vector based on the length of the parameters vector for every proxy + std::vector &)> > fittingFunction; + for (int l = 0; l < myArea.getParameters().size(); l++) + { + if (myArea.getParameters()[l].size() == 2) + fittingFunction.push_back(functionLinear_intercept); + else if (myArea.getParameters()[l].size() == 4) + fittingFunction.push_back(lapseRatePiecewise_two); + else if (myArea.getParameters()[l].size() == 5) + fittingFunction.push_back(lapseRatePiecewise_three); + else if (myArea.getParameters()[l].size() == 6) + fittingFunction.push_back(lapseRatePiecewise_three_free); + } + interpolationSettings.setFittingFunction(fittingFunction); + + //create vector of macro area interpolation points + std::vector temp = myArea.getMeteoPoints(); + for (int l = 0; l < temp.size(); l++) + { + for (int k = 0; k < interpolationPoints.size(); k++) + if (interpolationPoints[k].index == temp[l]) + subsetInterpolationPoints.push_back(interpolationPoints[k]); + } + + //detrending + if (elevationPos != NODATA && myArea.getCombination().isProxyActive(elevationPos)) + detrendingElevation(elevationPos, subsetInterpolationPoints, &interpolationSettings); + + detrendingOtherProxies(elevationPos, subsetInterpolationPoints, &interpolationSettings); + + Crit3DMeteoPoint* myMeteoPoints = new Crit3DMeteoPoint[unsigned(myArea.getMeteoPoints().size())]; + std::vector meteoPointsList = myArea.getMeteoPoints(); + + for (int i = 0; i < meteoPointsList.size(); i++) + myMeteoPoints[i] = meteoPoints[meteoPointsList[i]]; + + if (interpolationSettings.getUseTD() && getUseTdVar(myVar)) + topographicDistanceOptimize(myVar, myMeteoPoints, myArea.getMeteoPoints().size(), subsetInterpolationPoints, &interpolationSettings, meteoSettings); + + myMeteoPoints->clear(); + + return; +} + float Project::meteoDataConsistency(meteoVariable myVar, const Crit3DTime& timeIni, const Crit3DTime& timeFin) { diff --git a/project/project.h b/project/project.h index 42d2386bc..75cf4c809 100644 --- a/project/project.h +++ b/project/project.h @@ -267,6 +267,8 @@ bool loadGlocalStationsCsv(QString fileName, std::vector> &areaPoints); bool groupCellsInArea(std::vector &areaPoints, unsigned index, bool isGrid); bool glocalWeightsMaps(float windowWidth); + void macroAreaDetrending(Crit3DMacroArea myArea, meteoVariable myVar, std::vector interpolationPoints, std::vector &subsetInterpolationPoints, int elevationPos); + bool checkInterpolation(meteoVariable myVar); bool checkInterpolationGrid(meteoVariable myVar); From 2a540f66f0cfa157a365aea07b7a2b6352368657 Mon Sep 17 00:00:00 2001 From: Gabriele Antolini Date: Mon, 4 Nov 2024 18:41:01 +0100 Subject: [PATCH 19/34] revised interpolation grid period --- meteo/meteo.cpp | 3 + meteo/meteoGrid.cpp | 37 ++++- meteo/meteoGrid.h | 4 +- meteo/meteoPoint.cpp | 71 +++++--- meteo/meteoPoint.h | 3 +- pragaProject/pragaProject.cpp | 302 +++++++++++++++++----------------- pragaProject/pragaProject.h | 4 +- pragaProject/pragaShell.cpp | 27 ++- 8 files changed, 268 insertions(+), 183 deletions(-) diff --git a/meteo/meteo.cpp b/meteo/meteo.cpp index 2c930b88b..514bf3749 100644 --- a/meteo/meteo.cpp +++ b/meteo/meteo.cpp @@ -631,6 +631,9 @@ double ET0_Penman_hourly_net_rad(double heigth, double netIrradiance, double air double gamma; /*!< psychrometric constant (kPa C-1) */ double firstTerm, secondTerm, denominator; + if (heigth == NODATA || netIrradiance == NODATA || airTemp == NODATA || airHum == NODATA || windSpeed10 == NODATA) + return NODATA; + netRadiation = 3600 * netIrradiance; es = saturationVaporPressure(airTemp) / 1000.; diff --git a/meteo/meteoGrid.cpp b/meteo/meteoGrid.cpp index acf77e5e1..f3ca36391 100644 --- a/meteo/meteoGrid.cpp +++ b/meteo/meteoGrid.cpp @@ -794,6 +794,24 @@ void Crit3DMeteoGrid::emptyGridData(Crit3DDate dateIni, Crit3DDate dateFin) } } +void Crit3DMeteoGrid::computeRelativeHumidityFromTd(const Crit3DDate myDate, const int myHour) +{ + float t,td,rh; + + for (unsigned row = 0; row < unsigned(gridStructure().header().nrRows); row++) + for (unsigned col = 0; col < unsigned(gridStructure().header().nrCols); col++) + { + t = _meteoPoints[row][col]->getMeteoPointValueH(myDate, myHour, 0, airTemperature); + td = _meteoPoints[row][col]->getMeteoPointValueH(myDate, myHour, 0, airDewTemperature); + + if (! isEqual(t, NODATA) && ! isEqual(td, NODATA)) + { + rh = relHumFromTdew(td, t); + _meteoPoints[row][col]->setMeteoPointValueH(myDate, myHour, 0, airRelHumidity, rh); + } + } +} + void Crit3DMeteoGrid::computeWindVectorHourly(const Crit3DDate myDate, const int myHour) { float intensity = NODATA, direction = NODATA; @@ -1012,7 +1030,22 @@ void Crit3DMeteoGrid::saveRowColfromZone(gis::Crit3DRasterGrid* zoneGrid, std::v } -void Crit3DMeteoGrid::computeHourlyDerivedVariables(Crit3DTime dateTime) +void Crit3DMeteoGrid::computeHourlyDerivedVariables(Crit3DTime dateTime, std::vector variables, bool useNetRad) +{ + + for (unsigned row = 0; row < unsigned(gridStructure().header().nrRows); row++) + { + for (unsigned col = 0; col < unsigned(gridStructure().header().nrCols); col++) + { + if (_meteoPoints[row][col]->active) + { + _meteoPoints[row][col]->computeHourlyDerivedVariables(dateTime, variables, useNetRad); + } + } + } +} + +void Crit3DMeteoGrid::computeDailyDerivedVariables(Crit3DDate date, std::vector variables, Crit3DMeteoSettings& meteoSettings) { for (unsigned row = 0; row < unsigned(gridStructure().header().nrRows); row++) @@ -1021,7 +1054,7 @@ void Crit3DMeteoGrid::computeHourlyDerivedVariables(Crit3DTime dateTime) { if (_meteoPoints[row][col]->active) { - _meteoPoints[row][col]->computeDerivedVariables(dateTime); + _meteoPoints[row][col]->computeDailyDerivedVariables(date, variables, meteoSettings); } } } diff --git a/meteo/meteoGrid.h b/meteo/meteoGrid.h index c16fe4ea7..b89c1a490 100644 --- a/meteo/meteoGrid.h +++ b/meteo/meteoGrid.h @@ -147,8 +147,10 @@ void saveRowColfromZone(gis::Crit3DRasterGrid* zoneGrid, std::vector > &meteoGridRow, std::vector > &meteoGridCol); + void computeRelativeHumidityFromTd(const Crit3DDate myDate, const int myHour); void computeWindVectorHourly(const Crit3DDate myDate, const int myHour); - void computeHourlyDerivedVariables(Crit3DTime dateTime); + void computeHourlyDerivedVariables(Crit3DTime dateTime, std::vector variables, bool useNetRad); + void computeDailyDerivedVariables(Crit3DDate date, std::vector variables, Crit3DMeteoSettings &meteoSettings); private: diff --git a/meteo/meteoPoint.cpp b/meteo/meteoPoint.cpp index e32f7f411..5fe71f248 100644 --- a/meteo/meteoPoint.cpp +++ b/meteo/meteoPoint.cpp @@ -1218,38 +1218,69 @@ std::vector Crit3DMeteoPoint::getProxyValues() return myValues; } -bool Crit3DMeteoPoint::computeDerivedVariables(Crit3DTime dateTime) +bool Crit3DMeteoPoint::computeHourlyDerivedVariables(Crit3DTime dateTime, std::vector variables, bool useNetRad) { - short leafW; - float temperature, windSpeed, height, netRadiation; - Crit3DDate myDate = dateTime.date; int myHour = dateTime.getHour(); - bool leafWres = false; - bool et0res = false; + for (auto var : variables) + { + if (var == leafWetness) + { + short leafW; + if (computeLeafWetness(getMeteoPointValueH(myDate, myHour, 0, precipitation), + getMeteoPointValueH(myDate, myHour, 0, airRelHumidity), &leafW)) + setMeteoPointValueH(myDate, myHour, 0, leafWetness, leafW); + } + else if (var == referenceEvapotranspiration) + { + float et0; - double relHumidity = double(getMeteoPointValueH(myDate, myHour, 0, airRelHumidity)); - double prec = double(getMeteoPointValueH(myDate, myHour, 0, precipitation)); + if (useNetRad) + { + et0 = float(ET0_Penman_hourly_net_rad(double(point.z), + double(getMeteoPointValueH(myDate, myHour, 0, netIrradiance)), + double(getMeteoPointValueH(myDate, myHour, 0, airTemperature)), + double(getMeteoPointValueH(myDate, myHour, 0, airRelHumidity)), + double(getMeteoPointValueH(myDate, myHour, 0, windScalarIntensity)))); + + } + else + { + //da sistemare + et0 = ET0_Penman_hourly(double(point.z), + double(getMeteoPointValueH(myDate, myHour, 0, atmTransmissivity) / 0.75), + double(getMeteoPointValueH(myDate, myHour, 0, globalIrradiance)), + double(getMeteoPointValueH(myDate, myHour, 0, airTemperature)), + double(getMeteoPointValueH(myDate, myHour, 0, airRelHumidity)), + double(getMeteoPointValueH(myDate, myHour, 0, windScalarIntensity))); + } - if (computeLeafWetness(prec, relHumidity, &leafW)) - leafWres = setMeteoPointValueH(myDate, myHour, 0, leafWetness, leafW); + setMeteoPointValueH(myDate, myHour, 0, referenceEvapotranspiration, et0); + } + } - temperature = getMeteoPointValueH(myDate, myHour, 0, airTemperature); - windSpeed = getMeteoPointValueH(myDate, myHour, 0, windScalarIntensity); - netRadiation = getMeteoPointValueH(myDate, myHour, 0, netIrradiance); - height = float(this->point.z); - float et0; + return true; +} - if (! isEqual(temperature, NODATA) && ! isEqual(relHumidity, NODATA) && ! isEqual(windSpeed, NODATA)) +bool Crit3DMeteoPoint::computeDailyDerivedVariables(Crit3DDate date, std::vector variables, Crit3DMeteoSettings& meteoSettings) +{ + for (auto var : variables) { - et0 = float(ET0_Penman_hourly_net_rad(double(height), double(netRadiation), - double(temperature), double(relHumidity), double(windSpeed))); - et0res = setMeteoPointValueH(myDate, myHour, 0, referenceEvapotranspiration, et0); + if (var == dailyReferenceEvapotranspirationHS) + { + float et0 = dailyEtpHargreaves(getMeteoPointValueD(date, dailyAirTemperatureMin), + getMeteoPointValueD(date, dailyAirTemperatureMax), + date, latitude, &meteoSettings); + + setMeteoPointValueD(date, dailyReferenceEvapotranspirationHS, et0); + } } - return (leafWres && et0res); + + return true; } + bool Crit3DMeteoPoint::computeMonthlyAggregate(Crit3DDate firstDate, Crit3DDate lastDate, meteoVariable dailyMeteoVar, Crit3DMeteoSettings* meteoSettings, Crit3DQuality* qualityCheck, Crit3DClimateParameters* climateParam) diff --git a/meteo/meteoPoint.h b/meteo/meteoPoint.h index 2c86a6da4..f8480069e 100644 --- a/meteo/meteoPoint.h +++ b/meteo/meteoPoint.h @@ -171,7 +171,8 @@ void setAltitude(double altitude); void setLapseRateCode(std::string lapseRateCode); - bool computeDerivedVariables(Crit3DTime dateTime); + bool computeHourlyDerivedVariables(Crit3DTime dateTime, std::vector variables, bool useNetRad); + bool computeDailyDerivedVariables(Crit3DDate date, std::vector variables, Crit3DMeteoSettings &meteoSettings); bool computeMonthlyAggregate(Crit3DDate firstDate, Crit3DDate lastDate, meteoVariable dailyMeteoVar, Crit3DMeteoSettings *meteoSettings, Crit3DQuality *qualityCheck, Crit3DClimateParameters *climateParam); TObsDataH *getObsDataH() const; void initializeObsDataDFromMp(unsigned int numberOfDays, const Crit3DDate& firstDate, Crit3DMeteoPoint mp); diff --git a/pragaProject/pragaProject.cpp b/pragaProject/pragaProject.cpp index 9c379fd63..16d416903 100644 --- a/pragaProject/pragaProject.cpp +++ b/pragaProject/pragaProject.cpp @@ -2241,7 +2241,7 @@ bool PragaProject::timeAggregateGrid(QDate dateIni, QDate dateFin, QList variables, bool useNetRad) { // check meteo grid @@ -2261,32 +2261,66 @@ bool PragaProject::hourlyDerivedVariablesGrid(QDate first, QDate last, bool load QDateTime firstDateTime = QDateTime(first, QTime(1,0), Qt::UTC); QDateTime lastDateTime = QDateTime(last.addDays(1), QTime(0,0), Qt::UTC); - // now only hourly-->daily - if (loadData) + std::vector hourlyVars; + std::vector dailyVars; + + bool isHourly = false, isDaily = false; + meteoVariable myVar; + frequencyType freq; + + // find out needed frequency + foreach (myVar, variables) { - logInfoGUI("Loading grid data: " + first.toString("yyyy-MM-dd") + "-" + last.toString("yyyy-MM-dd")); - loadMeteoGridHourlyData(firstDateTime, lastDateTime, false); + freq = getVarFrequency(myVar); + + if (freq == noFrequency) + { + errorString = "Unknown variable: " + QString::fromStdString(getMeteoVarName(myVar)); + return false; + } + else if (freq == hourly) + { + isHourly = true; + hourlyVars.push_back(myVar); + } + else if (freq == daily) + { + isDaily = true; + dailyVars.push_back(myVar); + } } - while(firstDateTime <= lastDateTime) + logInfoGUI("Loading grid data: " + first.toString("yyyy-MM-dd") + "-" + last.toString("yyyy-MM-dd")); + + if (isHourly) loadMeteoGridHourlyData(firstDateTime, lastDateTime, false); + if (isDaily) loadMeteoGridDailyData(firstDateTime.date(), lastDateTime.date(), false); + + if (isHourly) { - meteoGridDbHandler->meteoGrid()->computeHourlyDerivedVariables(getCrit3DTime(firstDateTime)); - firstDateTime = firstDateTime.addSecs(3600); + while(firstDateTime <= lastDateTime) + { + meteoGridDbHandler->meteoGrid()->computeHourlyDerivedVariables(getCrit3DTime(firstDateTime), hourlyVars, useNetRad); + firstDateTime = firstDateTime.addSecs(3600); + } } - firstDateTime = QDateTime(first, QTime(1,0), Qt::UTC); - // saving hourly meteo grid data to DB - if (saveData) + if (isDaily) { - - // save derived variables - QList variables; - variables << leafWetness << referenceEvapotranspiration; - QString myError; - logInfoGUI("Saving meteo grid data"); - if (! meteoGridDbHandler->saveGridData(&myError, firstDateTime, lastDateTime, variables, meteoSettings)) return false; + QDate myDate = firstDateTime.date(); + while (myDate <= lastDateTime.date()) + { + meteoGridDbHandler->meteoGrid()->computeDailyDerivedVariables(getCrit3DDate(myDate), dailyVars, *meteoSettings); + myDate = myDate.addDays(1); + } } + firstDateTime = QDateTime(first, QTime(1,0), Qt::UTC); + // save derived variables + variables << leafWetness << referenceEvapotranspiration; + QString myError; + logInfoGUI("Saving meteo grid data"); + if (! meteoGridDbHandler->saveGridData(&myError, firstDateTime, lastDateTime, variables, meteoSettings)) return false; + return true; } @@ -2533,9 +2567,101 @@ bool PragaProject::interpolationOutputPointsPeriod(QDate firstDate, QDate lastDa return true; } +bool PragaProject::interpolationMeteoGrid(meteoVariable myVar, frequencyType myFrequency, const Crit3DTime& myTime) +{ + if (meteoGridDbHandler == nullptr) + { + errorString = "Open a Meteo Grid before."; + return false; + } + + if (interpolationSettings.getMeteoGridUpscaleFromDem()) + { + if (myFrequency == hourly) + { + if (myVar == airRelHumidity && interpolationSettings.getUseDewPoint()) + { + passInterpolatedTemperatureToHumidityPoints(myTime, meteoSettings); + if (! interpolationDemMain(airDewTemperature, myTime, hourlyMeteoMaps->mapHourlyTdew)) return false; + hourlyMeteoMaps->computeRelativeHumidityMap(hourlyMeteoMaps->mapHourlyRelHum); + meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, myFrequency, myTime.date, myTime.getHour(), + myTime.getMinutes(), &DEM, hourlyMeteoMaps->mapHourlyRelHum, interpolationSettings.getMeteoGridAggrMethod()); + + } + else if (myVar == windVectorDirection || myVar == windVectorIntensity) + { + if (! interpolationDemMain(windVectorX, myTime, getPragaMapFromVar(windVectorX))) return false; + if (! interpolationDemMain(windVectorY, myTime, getPragaMapFromVar(windVectorY))) return false; + meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(windVectorX, hourly, myTime.date, myTime.getHour(), 0, &DEM, getPragaMapFromVar(windVectorX), interpolationSettings.getMeteoGridAggrMethod()); + meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(windVectorY, hourly, myTime.date, myTime.getHour(), 0, &DEM, getPragaMapFromVar(windVectorY), interpolationSettings.getMeteoGridAggrMethod()); + meteoGridDbHandler->meteoGrid()->computeWindVectorHourly(myTime.date, myTime.getHour()); + } + else if (myVar == leafWetness) + { + hourlyMeteoMaps->computeLeafWetnessMap() ; + meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, myFrequency, myTime.date, myTime.getHour(), + myTime.getMinutes(), &DEM, hourlyMeteoMaps->mapHourlyLeafW, interpolationSettings.getMeteoGridAggrMethod()); + + } + else if (myVar == referenceEvapotranspiration) + { + hourlyMeteoMaps->computeET0PMMap(DEM, radiationMaps); + meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, myFrequency, myTime.date, myTime.getHour(), + myTime.getMinutes(), &DEM, hourlyMeteoMaps->mapHourlyET0, interpolationSettings.getMeteoGridAggrMethod()); + } + else + { + if (!interpolationDemMain(myVar, myTime, getPragaMapFromVar(myVar))) return false; + } + + meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, hourly, myTime.date, 0, 0, &DEM, getPragaMapFromVar(myVar), interpolationSettings.getMeteoGridAggrMethod()); + + } + else if (myFrequency == daily) + { + if (myVar == dailyReferenceEvapotranspirationHS) { + pragaDailyMaps->computeHSET0Map(&gisSettings, myTime.date);} + else if (myVar == dailyAirTemperatureMax || myVar == dailyAirTemperatureMin) { + if (! pragaDailyMaps->fixDailyThermalConsistency()) return false;} + else + { + if (! interpolationDemMain(myVar, myTime, getPragaMapFromVar(myVar))) + return false; + } + + meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, daily, myTime.date, 0, 0, &DEM, getPragaMapFromVar(myVar), interpolationSettings.getMeteoGridAggrMethod()); + } + } + else + { + if (myFrequency == hourly) + { + if (myVar == airRelHumidity && interpolationSettings.getUseDewPoint()) + { + passGridTemperatureToHumidityPoints(myTime, meteoSettings); + if (! interpolationGrid(airDewTemperature, myTime)) return false; + meteoGridDbHandler->meteoGrid()->computeRelativeHumidityFromTd(myTime.date, myTime.getHour()); + } + else if (myVar == windVectorDirection || myVar == windVectorIntensity) + { + if (! interpolationMeteoGrid(windVectorX, hourly, myTime)) return false; + if (! interpolationMeteoGrid(windVectorY, hourly, myTime)) return false; + meteoGridDbHandler->meteoGrid()->computeWindVectorHourly(myTime.date, myTime.getHour()); + } + else + { + if (! interpolationGrid(myVar, myTime)) return false; + } + } + } + + meteoGridDbHandler->meteoGrid()->fillMeteoRaster(); + + return true; +} bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QList variables, - QList aggrVariables, bool saveRasters, + QList aggrVariables, int nrDaysLoading, int nrDaysSaving) { // check variables @@ -2571,7 +2697,6 @@ bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QL QString myError, rasterName, varName; int myHour; QDate myDate = dateIni; - gis::Crit3DRasterGrid* myGrid; meteoVariable myVar; frequencyType freq; bool isDaily = false, isHourly = false; @@ -2688,13 +2813,8 @@ bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QL if (useProxies && currentYear != myDate.year()) { logInfoGUI("Interpolating proxy grid series..."); - - if (! checkProxyGridSeries(&interpolationSettings, DEM, proxyGridSeries, myDate, errorString)) - return false; - - if (! readProxyValues()) - return false; - + if (! checkProxyGridSeries(&interpolationSettings, DEM, proxyGridSeries, myDate, errorString)) return false; + if (! readProxyValues()) return false; currentYear = myDate.year(); } @@ -2708,74 +2828,8 @@ bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QL { if (getVarFrequency(myVar) == hourly) { - varName = QString::fromStdString(getMeteoVarName(myVar)); - logInfo(varName); - - if (myVar == airRelHumidity && interpolationSettings.getUseDewPoint()) - { - if (interpolationSettings.getUseInterpolatedTForRH()) - { - if (interpolationSettings.getMeteoGridUpscaleFromDem()) - { - passInterpolatedTemperatureToHumidityPoints(getCrit3DTime(myDate, myHour), meteoSettings); - if (! interpolationDemMain(airDewTemperature, getCrit3DTime(myDate, myHour), hourlyMeteoMaps->mapHourlyTdew)) return false; - hourlyMeteoMaps->computeRelativeHumidityMap(hourlyMeteoMaps->mapHourlyRelHum); - } - else - { - passGridTemperatureToHumidityPoints(getCrit3DTime(myDate, myHour), meteoSettings); - if (! interpolationMeteoGrid(airDewTemperature, hourly, getCrit3DTime(myDate, myHour))) return false; - } - } - - if (interpolationSettings.getMeteoGridUpscaleFromDem() && saveRasters) - { - myGrid = getPragaMapFromVar(airDewTemperature); - rasterName = getMapFileOutName(airDewTemperature, myDate, myHour); - if (rasterName != "") gis::writeEsriGrid(getProjectPath().toStdString() + rasterName.toStdString(), myGrid, errString); - } - } - else if (myVar == windVectorDirection || myVar == windVectorIntensity) { - if (interpolationSettings.getMeteoGridUpscaleFromDem()) - { - if (! interpolationDemMain(windVectorX, getCrit3DTime(myDate, myHour), getPragaMapFromVar(windVectorX))) return false; - if (! interpolationDemMain(windVectorY, getCrit3DTime(myDate, myHour), getPragaMapFromVar(windVectorY))) return false; - if (! pragaHourlyMaps->computeWindVector()) return false; - } - else - { - if (! interpolationMeteoGrid(windVectorX, hourly, getCrit3DTime(myDate, myHour))) return false; - if (! interpolationMeteoGrid(windVectorY, hourly, getCrit3DTime(myDate, myHour))) return false; - } - } - else if (myVar == leafWetness) { - hourlyMeteoMaps->computeLeafWetnessMap() ; - } - else if (myVar == referenceEvapotranspiration) { - hourlyMeteoMaps->computeET0PMMap(DEM, radiationMaps); - } - else { - if (! interpolationDemMain(myVar, getCrit3DTime(myDate, myHour), getPragaMapFromVar(myVar))) return false; - } - - myGrid = getPragaMapFromVar(myVar); - if (myGrid == nullptr) return false; - - //save raster - if (saveRasters) - { - rasterName = getMapFileOutName(myVar, myDate, myHour); - if (rasterName != "") gis::writeEsriGrid(getProjectPath().toStdString() + rasterName.toStdString(), myGrid, errString); - } - - if (myVar == windVectorDirection || myVar == windVectorIntensity) - { - meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(windVectorX, hourly, getCrit3DDate(myDate), myHour, 0, &DEM, getPragaMapFromVar(windVectorX), interpolationSettings.getMeteoGridAggrMethod()); - meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(windVectorY, hourly, getCrit3DDate(myDate), myHour, 0, &DEM, getPragaMapFromVar(windVectorY), interpolationSettings.getMeteoGridAggrMethod()); - meteoGridDbHandler->meteoGrid()->computeWindVectorHourly(getCrit3DDate(myDate), myHour); - } - else - meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, hourly, getCrit3DDate(myDate), myHour, 0, &DEM, myGrid, interpolationSettings.getMeteoGridAggrMethod()); + logInfo(QString::fromStdString(getMeteoVarName(myVar))); + if (! interpolationMeteoGrid(myVar, hourly, getCrit3DTime(myDate, myHour))) return false; } } } @@ -2789,32 +2843,8 @@ bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QL { if (getVarFrequency(myVar) == daily) { - varName = QString::fromStdString(getMeteoVarName(myVar)); - logInfo(varName); - - if (myVar == dailyReferenceEvapotranspirationHS) { - pragaDailyMaps->computeHSET0Map(&gisSettings, getCrit3DDate(myDate)); - } - else { - if (! interpolationDemMain(myVar, getCrit3DTime(myDate, 0), getPragaMapFromVar(myVar))) return false; - } - - // fix daily temperatures consistency - if (myVar == dailyAirTemperatureMax || myVar == dailyAirTemperatureMin) { - if (! pragaDailyMaps->fixDailyThermalConsistency()) return false; - } - - myGrid = getPragaMapFromVar(myVar); - if (myGrid == nullptr) return false; - - //save raster - if (saveRasters) - { - rasterName = getMapFileOutName(myVar, myDate, 0); - gis::writeEsriGrid(getProjectPath().toStdString() + rasterName.toStdString(), myGrid, errString); - } - - meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, daily, getCrit3DDate(myDate), 0, 0, &DEM, myGrid, interpolationSettings.getMeteoGridAggrMethod()); + logInfo(QString::fromStdString(getMeteoVarName(myVar))); + if (! interpolationMeteoGrid(myVar, daily, getCrit3DTime(myDate, myHour))) return false; } } } @@ -2947,36 +2977,6 @@ bool PragaProject::interpolationCrossValidationPeriod(QDate dateIni, QDate dateF return true; } -bool PragaProject::interpolationMeteoGrid(meteoVariable myVar, frequencyType myFrequency, const Crit3DTime& myTime) -{ - if (meteoGridDbHandler != nullptr) - { - if (interpolationSettings.getMeteoGridUpscaleFromDem()) - { - gis::Crit3DRasterGrid *myRaster = new gis::Crit3DRasterGrid; - if (!interpolationDemMain(myVar, myTime, myRaster)) return false; - - meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, myFrequency, myTime.date, myTime.getHour(), - myTime.getMinutes(), &DEM, myRaster, interpolationSettings.getMeteoGridAggrMethod()); - } - else - { - if (! interpolationGrid(myVar, myTime)) - return false; - } - - meteoGridDbHandler->meteoGrid()->fillMeteoRaster(); - } - else - { - errorString = "Open a Meteo Grid before."; - return false; - } - - return true; -} - - bool PragaProject::dbMeteoPointDataCount(QDate myFirstDate, QDate myLastDate, meteoVariable myVar, QString dataset, std::vector &myCounter) { frequencyType myFreq = getVarFrequency(myVar); diff --git a/pragaProject/pragaProject.h b/pragaProject/pragaProject.h index 3c7eb83c6..cea938f6d 100644 --- a/pragaProject/pragaProject.h +++ b/pragaProject/pragaProject.h @@ -111,13 +111,13 @@ bool interpolationMeteoGrid(meteoVariable myVar, frequencyType myFrequency, const Crit3DTime& myTime); bool interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QList variables, - QList aggrVariables, bool saveRasters, int nrDaysLoading, int nrDaysSaving); + QList aggrVariables, int nrDaysLoading, int nrDaysSaving); bool interpolationCrossValidationPeriod(QDate dateIni, QDate dateFin, meteoVariable myVar, QString filename); bool saveGrid(meteoVariable myVar, frequencyType myFrequency, const Crit3DTime& myTime, bool showInfo); bool timeAggregateGridVarHourlyInDaily(meteoVariable dailyVar, Crit3DDate dateIni, Crit3DDate dateFin); bool timeAggregateGrid(QDate dateIni, QDate dateFin, QList variables, bool loadData, bool saveData); bool computeDailyVariablesPoint(Crit3DMeteoPoint *meteoPoint, QDate first, QDate last, QList variables); - bool hourlyDerivedVariablesGrid(QDate first, QDate last, bool loadData, bool saveData); + bool derivedVariablesGrid(QDate first, QDate last, QList variables, bool useNetRad); bool elaborationPointsCycle(bool isAnomaly, bool showInfo); bool elaborationPointsCycleGrid(bool isAnomaly, bool showInfo); bool elaborationCheck(bool isMeteoGrid, bool isAnomaly); diff --git a/pragaProject/pragaShell.cpp b/pragaProject/pragaShell.cpp index f4a28c316..36245d408 100644 --- a/pragaProject/pragaShell.cpp +++ b/pragaProject/pragaShell.cpp @@ -286,7 +286,6 @@ int cmdInterpolationGridPeriod(PragaProject* myProject, QList argumentL } QDate dateIni, dateFin; - bool saveRasters = false; QList varString; QList variables, aggrVariables; QString var; @@ -341,8 +340,6 @@ int cmdInterpolationGridPeriod(PragaProject* myProject, QList argumentL dateFin = QDate::currentDate().addDays(-1); dateIni = dateFin.addDays(-6); } - else if (argumentList.at(i).left(2) == "-r") - saveRasters = true; else if (argumentList.at(i).left(3) == "-s:") saveInterval = argumentList[i].right(argumentList[i].length()-3).toInt(&parseSaveInterval); else if (argumentList.at(i).left(3) == "-l:") @@ -374,7 +371,7 @@ int cmdInterpolationGridPeriod(PragaProject* myProject, QList argumentL return PRAGA_INVALID_COMMAND; } - if (! myProject->interpolationMeteoGridPeriod(dateIni, dateFin, variables, aggrVariables, saveRasters, loadInterval, saveInterval)) + if (! myProject->interpolationMeteoGridPeriod(dateIni, dateFin, variables, aggrVariables, loadInterval, saveInterval)) return PRAGA_ERROR; return PRAGA_OK; @@ -445,10 +442,24 @@ int cmdHourlyDerivedVariablesGrid(PragaProject* myProject, QList argume // default date QDate first = QDate::currentDate(); QDate last = first.addDays(9); + QList variables; + QList varString; + QString var; + meteoVariable meteoVar; + bool useNetRad = false; for (int i = 1; i < argumentList.size(); i++) { - if (argumentList.at(i).left(4) == "-d1:") + if (argumentList[i].left(3) == "-v:") + { + varString = argumentList[i].right(argumentList[i].length()-3).split(","); + foreach (var,varString) + { + meteoVar = getMeteoVar(var.toStdString()); + if (meteoVar != noMeteoVar) variables << meteoVar; + } + } + else if (argumentList.at(i).left(4) == "-d1:") { QString dateIniStr = argumentList[i].right(argumentList[i].length()-4); first = QDate::fromString(dateIniStr, "dd/MM/yyyy"); @@ -458,6 +469,10 @@ int cmdHourlyDerivedVariablesGrid(PragaProject* myProject, QList argume QString dateFinStr = argumentList[i].right(argumentList[i].length()-4); last = QDate::fromString(dateFinStr, "dd/MM/yyyy"); } + else if (argumentList.at(i).left(2) == "-n") + { + useNetRad = true; + } } @@ -473,7 +488,7 @@ int cmdHourlyDerivedVariablesGrid(PragaProject* myProject, QList argume return PRAGA_INVALID_COMMAND; } - if (! myProject->hourlyDerivedVariablesGrid(first, last, true, true)) + if (! myProject->derivedVariablesGrid(first, last, variables, useNetRad)) return PRAGA_ERROR; return PRAGA_OK; From d10e9f5b3415ae5d5a2db7254ec3cc3e0088bde8 Mon Sep 17 00:00:00 2001 From: Gabriele Antolini Date: Tue, 5 Nov 2024 11:58:40 +0100 Subject: [PATCH 20/34] fixed interpolation grid period --- pragaProject/pragaProject.cpp | 33 ++++++++++++++++++--------------- 1 file changed, 18 insertions(+), 15 deletions(-) diff --git a/pragaProject/pragaProject.cpp b/pragaProject/pragaProject.cpp index 16d416903..a8739e687 100644 --- a/pragaProject/pragaProject.cpp +++ b/pragaProject/pragaProject.cpp @@ -2584,52 +2584,55 @@ bool PragaProject::interpolationMeteoGrid(meteoVariable myVar, frequencyType myF passInterpolatedTemperatureToHumidityPoints(myTime, meteoSettings); if (! interpolationDemMain(airDewTemperature, myTime, hourlyMeteoMaps->mapHourlyTdew)) return false; hourlyMeteoMaps->computeRelativeHumidityMap(hourlyMeteoMaps->mapHourlyRelHum); - meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, myFrequency, myTime.date, myTime.getHour(), - myTime.getMinutes(), &DEM, hourlyMeteoMaps->mapHourlyRelHum, interpolationSettings.getMeteoGridAggrMethod()); + meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, myFrequency, myTime.date, myTime.getHour(), myTime.getMinutes(), + &DEM, hourlyMeteoMaps->mapHourlyRelHum, interpolationSettings.getMeteoGridAggrMethod()); } else if (myVar == windVectorDirection || myVar == windVectorIntensity) { if (! interpolationDemMain(windVectorX, myTime, getPragaMapFromVar(windVectorX))) return false; if (! interpolationDemMain(windVectorY, myTime, getPragaMapFromVar(windVectorY))) return false; - meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(windVectorX, hourly, myTime.date, myTime.getHour(), 0, &DEM, getPragaMapFromVar(windVectorX), interpolationSettings.getMeteoGridAggrMethod()); - meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(windVectorY, hourly, myTime.date, myTime.getHour(), 0, &DEM, getPragaMapFromVar(windVectorY), interpolationSettings.getMeteoGridAggrMethod()); + meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(windVectorX, hourly, myTime.date, myTime.getHour(), myTime.getMinutes(), + &DEM, getPragaMapFromVar(windVectorX), interpolationSettings.getMeteoGridAggrMethod()); + meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(windVectorY, hourly, myTime.date, myTime.getHour(), myTime.getMinutes(), + &DEM, getPragaMapFromVar(windVectorY), interpolationSettings.getMeteoGridAggrMethod()); meteoGridDbHandler->meteoGrid()->computeWindVectorHourly(myTime.date, myTime.getHour()); } else if (myVar == leafWetness) { hourlyMeteoMaps->computeLeafWetnessMap() ; - meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, myFrequency, myTime.date, myTime.getHour(), - myTime.getMinutes(), &DEM, hourlyMeteoMaps->mapHourlyLeafW, interpolationSettings.getMeteoGridAggrMethod()); + meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, myFrequency, myTime.date, myTime.getHour(), myTime.getMinutes(), + &DEM, hourlyMeteoMaps->mapHourlyLeafW, interpolationSettings.getMeteoGridAggrMethod()); } else if (myVar == referenceEvapotranspiration) { hourlyMeteoMaps->computeET0PMMap(DEM, radiationMaps); - meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, myFrequency, myTime.date, myTime.getHour(), - myTime.getMinutes(), &DEM, hourlyMeteoMaps->mapHourlyET0, interpolationSettings.getMeteoGridAggrMethod()); + meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, myFrequency, myTime.date, myTime.getHour(), myTime.getMinutes(), + &DEM, hourlyMeteoMaps->mapHourlyET0, interpolationSettings.getMeteoGridAggrMethod()); } else { if (!interpolationDemMain(myVar, myTime, getPragaMapFromVar(myVar))) return false; + meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, hourly, myTime.date, myTime.getHour(), myTime.getMinutes(), + &DEM, getPragaMapFromVar(myVar), interpolationSettings.getMeteoGridAggrMethod()); } - - meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, hourly, myTime.date, 0, 0, &DEM, getPragaMapFromVar(myVar), interpolationSettings.getMeteoGridAggrMethod()); - } else if (myFrequency == daily) { if (myVar == dailyReferenceEvapotranspirationHS) { pragaDailyMaps->computeHSET0Map(&gisSettings, myTime.date);} - else if (myVar == dailyAirTemperatureMax || myVar == dailyAirTemperatureMin) { - if (! pragaDailyMaps->fixDailyThermalConsistency()) return false;} else { if (! interpolationDemMain(myVar, myTime, getPragaMapFromVar(myVar))) return false; } - meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, daily, myTime.date, 0, 0, &DEM, getPragaMapFromVar(myVar), interpolationSettings.getMeteoGridAggrMethod()); + if (myVar == dailyAirTemperatureMax || myVar == dailyAirTemperatureMin) { + if (! pragaDailyMaps->fixDailyThermalConsistency()) return false;} + + meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, daily, myTime.date, myTime.getHour(), myTime.getMinutes(), + &DEM, getPragaMapFromVar(myVar), interpolationSettings.getMeteoGridAggrMethod()); } } else @@ -2694,7 +2697,7 @@ bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QL // order variables for derived computation std::string errString; - QString myError, rasterName, varName; + QString myError; int myHour; QDate myDate = dateIni; meteoVariable myVar; From 3b4a4daa63ce9e00726005db12616313c5a5a763 Mon Sep 17 00:00:00 2001 From: Gabriele Antolini Date: Tue, 5 Nov 2024 13:01:07 +0100 Subject: [PATCH 21/34] grid daily thermal consistency --- meteo/meteoGrid.cpp | 20 ++++++++++++++++++++ meteo/meteoGrid.h | 1 + pragaProject/pragaProject.cpp | 7 +++++++ project/meteoMaps.cpp | 2 ++ 4 files changed, 30 insertions(+) diff --git a/meteo/meteoGrid.cpp b/meteo/meteoGrid.cpp index f3ca36391..2b4cfdbe4 100644 --- a/meteo/meteoGrid.cpp +++ b/meteo/meteoGrid.cpp @@ -834,6 +834,26 @@ void Crit3DMeteoGrid::computeWindVectorHourly(const Crit3DDate myDate, const int } } +void Crit3DMeteoGrid::fixDailyThermalConsistency(const Crit3DDate myDate) +{ + float tmin = NODATA, tmax = NODATA; + + for (unsigned row = 0; row < unsigned(gridStructure().header().nrRows); row++) + for (unsigned col = 0; col < unsigned(gridStructure().header().nrCols); col++) + { + tmin = _meteoPoints[row][col]->getMeteoPointValueD(myDate, dailyAirTemperatureMin); + tmax = _meteoPoints[row][col]->getMeteoPointValueD(myDate, dailyAirTemperatureMax); + + if (! isEqual(tmin, NODATA) && ! isEqual(tmax, NODATA)) + { + if (tmin > tmax) + { + _meteoPoints[row][col]->setMeteoPointValueD(myDate, dailyAirTemperatureMin, tmax - 0.1); + } + } + } +} + void Crit3DMeteoGrid::spatialAggregateMeteoGrid(meteoVariable myVar, frequencyType freq, Crit3DDate date, int hour, int minute, gis::Crit3DRasterGrid* myDEM, gis::Crit3DRasterGrid *myRaster, aggregationMethod elab) { diff --git a/meteo/meteoGrid.h b/meteo/meteoGrid.h index b89c1a490..da4179b50 100644 --- a/meteo/meteoGrid.h +++ b/meteo/meteoGrid.h @@ -149,6 +149,7 @@ void computeRelativeHumidityFromTd(const Crit3DDate myDate, const int myHour); void computeWindVectorHourly(const Crit3DDate myDate, const int myHour); + void fixDailyThermalConsistency(const Crit3DDate myDate); void computeHourlyDerivedVariables(Crit3DTime dateTime, std::vector variables, bool useNetRad); void computeDailyDerivedVariables(Crit3DDate date, std::vector variables, Crit3DMeteoSettings &meteoSettings); diff --git a/pragaProject/pragaProject.cpp b/pragaProject/pragaProject.cpp index a8739e687..d1efa866c 100644 --- a/pragaProject/pragaProject.cpp +++ b/pragaProject/pragaProject.cpp @@ -2656,6 +2656,13 @@ bool PragaProject::interpolationMeteoGrid(meteoVariable myVar, frequencyType myF if (! interpolationGrid(myVar, myTime)) return false; } } + else if (myFrequency == daily) + { + if (! interpolationGrid(myVar, myTime)) return false; + if (myVar == dailyAirTemperatureMax || myVar == dailyAirTemperatureMin) { + meteoGridDbHandler->meteoGrid()->fixDailyThermalConsistency(myTime.date);} + + } } meteoGridDbHandler->meteoGrid()->fillMeteoRaster(); diff --git a/project/meteoMaps.cpp b/project/meteoMaps.cpp index e20d98631..0554dec20 100644 --- a/project/meteoMaps.cpp +++ b/project/meteoMaps.cpp @@ -102,6 +102,8 @@ bool Crit3DDailyMeteoMaps::computeHSET0Map(gis::Crit3DGisSettings* gisSettings, bool Crit3DDailyMeteoMaps::fixDailyThermalConsistency() { + // da migliorare! + if (! mapDailyTMax->isLoaded || ! mapDailyTMin->isLoaded) return true; if ( mapDailyTMax->getMapTime() != mapDailyTMin->getMapTime()) return true; From c75c805ba5c970d5586c369a7b84e256a7d38127 Mon Sep 17 00:00:00 2001 From: Gabriele Antolini Date: Tue, 5 Nov 2024 15:55:46 +0100 Subject: [PATCH 22/34] grid active cell for humidity --- meteo/meteoGrid.cpp | 47 +++++++++++++++++++++++++++------------------ project/project.cpp | 6 ++++-- 2 files changed, 32 insertions(+), 21 deletions(-) diff --git a/meteo/meteoGrid.cpp b/meteo/meteoGrid.cpp index 2b4cfdbe4..4a28ac979 100644 --- a/meteo/meteoGrid.cpp +++ b/meteo/meteoGrid.cpp @@ -801,13 +801,16 @@ void Crit3DMeteoGrid::computeRelativeHumidityFromTd(const Crit3DDate myDate, con for (unsigned row = 0; row < unsigned(gridStructure().header().nrRows); row++) for (unsigned col = 0; col < unsigned(gridStructure().header().nrCols); col++) { - t = _meteoPoints[row][col]->getMeteoPointValueH(myDate, myHour, 0, airTemperature); - td = _meteoPoints[row][col]->getMeteoPointValueH(myDate, myHour, 0, airDewTemperature); - - if (! isEqual(t, NODATA) && ! isEqual(td, NODATA)) + if (_meteoPoints[row][col]->active) { - rh = relHumFromTdew(td, t); - _meteoPoints[row][col]->setMeteoPointValueH(myDate, myHour, 0, airRelHumidity, rh); + t = _meteoPoints[row][col]->getMeteoPointValueH(myDate, myHour, 0, airTemperature); + td = _meteoPoints[row][col]->getMeteoPointValueH(myDate, myHour, 0, airDewTemperature); + + if (! isEqual(t, NODATA) && ! isEqual(td, NODATA)) + { + rh = relHumFromTdew(td, t); + _meteoPoints[row][col]->setMeteoPointValueH(myDate, myHour, 0, airRelHumidity, rh); + } } } } @@ -820,15 +823,18 @@ void Crit3DMeteoGrid::computeWindVectorHourly(const Crit3DDate myDate, const int for (unsigned row = 0; row < unsigned(gridStructure().header().nrRows); row++) for (unsigned col = 0; col < unsigned(gridStructure().header().nrCols); col++) { - u = _meteoPoints[row][col]->getMeteoPointValueH(myDate, myHour, 0, windVectorX); - v = _meteoPoints[row][col]->getMeteoPointValueH(myDate, myHour, 0, windVectorY); - - if (! isEqual(u, NODATA) && ! isEqual(v, NODATA)) + if (_meteoPoints[row][col]->active) { - if (computeWindPolar(u, v, &intensity, &direction)) + u = _meteoPoints[row][col]->getMeteoPointValueH(myDate, myHour, 0, windVectorX); + v = _meteoPoints[row][col]->getMeteoPointValueH(myDate, myHour, 0, windVectorY); + + if (! isEqual(u, NODATA) && ! isEqual(v, NODATA)) { - _meteoPoints[row][col]->setMeteoPointValueH(myDate, myHour, 0, windVectorIntensity, intensity); - _meteoPoints[row][col]->setMeteoPointValueH(myDate, myHour, 0, windVectorDirection, direction); + if (computeWindPolar(u, v, &intensity, &direction)) + { + _meteoPoints[row][col]->setMeteoPointValueH(myDate, myHour, 0, windVectorIntensity, intensity); + _meteoPoints[row][col]->setMeteoPointValueH(myDate, myHour, 0, windVectorDirection, direction); + } } } } @@ -841,14 +847,17 @@ void Crit3DMeteoGrid::fixDailyThermalConsistency(const Crit3DDate myDate) for (unsigned row = 0; row < unsigned(gridStructure().header().nrRows); row++) for (unsigned col = 0; col < unsigned(gridStructure().header().nrCols); col++) { - tmin = _meteoPoints[row][col]->getMeteoPointValueD(myDate, dailyAirTemperatureMin); - tmax = _meteoPoints[row][col]->getMeteoPointValueD(myDate, dailyAirTemperatureMax); - - if (! isEqual(tmin, NODATA) && ! isEqual(tmax, NODATA)) + if (_meteoPoints[row][col]->active) { - if (tmin > tmax) + tmin = _meteoPoints[row][col]->getMeteoPointValueD(myDate, dailyAirTemperatureMin); + tmax = _meteoPoints[row][col]->getMeteoPointValueD(myDate, dailyAirTemperatureMax); + + if (! isEqual(tmin, NODATA) && ! isEqual(tmax, NODATA)) { - _meteoPoints[row][col]->setMeteoPointValueD(myDate, dailyAirTemperatureMin, tmax - 0.1); + if (tmin > tmax) + { + _meteoPoints[row][col]->setMeteoPointValueD(myDate, dailyAirTemperatureMin, tmax - 0.1); + } } } } diff --git a/project/project.cpp b/project/project.cpp index 8343eedba..4c7bb2721 100644 --- a/project/project.cpp +++ b/project/project.cpp @@ -2546,9 +2546,11 @@ void Project::passGridTemperatureToHumidityPoints(Crit3DTime myTime, Crit3DMeteo if (! gis::isOutOfGridRowCol(row, col, meteoGridDbHandler->meteoGrid()->gridStructure().header())) { - meteoPoints[i].setMeteoPointValueH(myTime.date, myTime.getHour(), myTime.getMinutes(), airTemperature, + if (meteoGridDbHandler->meteoGrid()->meteoPoint(row,col).active) + { + meteoPoints[i].setMeteoPointValueH(myTime.date, myTime.getHour(), myTime.getMinutes(), airTemperature, meteoGridDbHandler->meteoGrid()->meteoPoint(row,col).getMeteoPointValueH(myTime.date, myTime.getHour(), myTime.getMinutes(), airTemperature)); - + } } } } From 36d81296b50d4643fa7f28c49c8421e903fe7611 Mon Sep 17 00:00:00 2001 From: Gabriele Antolini Date: Tue, 5 Nov 2024 16:34:08 +0100 Subject: [PATCH 23/34] interpolation grid radiation --- pragaProject/pragaProject.cpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/pragaProject/pragaProject.cpp b/pragaProject/pragaProject.cpp index d1efa866c..bf75dca6e 100644 --- a/pragaProject/pragaProject.cpp +++ b/pragaProject/pragaProject.cpp @@ -2651,6 +2651,13 @@ bool PragaProject::interpolationMeteoGrid(meteoVariable myVar, frequencyType myF if (! interpolationMeteoGrid(windVectorY, hourly, myTime)) return false; meteoGridDbHandler->meteoGrid()->computeWindVectorHourly(myTime.date, myTime.getHour()); } + else if (myVar == globalIrradiance) + { + // for radiation use dem and aggregate + if (! interpolateDemRadiation(myTime.addSeconds(-1800), radiationMaps->globalRadiationMap)) return false; + meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(globalIrradiance, daily, myTime.date, myTime.getHour(), myTime.getMinutes(), + &DEM, radiationMaps->globalRadiationMap, interpolationSettings.getMeteoGridAggrMethod()); + } else { if (! interpolationGrid(myVar, myTime)) return false; From d0c4d046d6cac0d1f45a63a4bda762c360061355 Mon Sep 17 00:00:00 2001 From: Gabriele Antolini Date: Tue, 5 Nov 2024 17:35:48 +0100 Subject: [PATCH 24/34] derived variables --- pragaProject/pragaProject.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pragaProject/pragaProject.cpp b/pragaProject/pragaProject.cpp index bf75dca6e..944847bf1 100644 --- a/pragaProject/pragaProject.cpp +++ b/pragaProject/pragaProject.cpp @@ -2306,8 +2306,8 @@ bool PragaProject::derivedVariablesGrid(QDate first, QDate last, QList meteoGrid()->computeDailyDerivedVariables(getCrit3DDate(myDate), dailyVars, *meteoSettings); myDate = myDate.addDays(1); From 13d3ba58018c89afc1cc4fbf2b3a9ea43d0ea322 Mon Sep 17 00:00:00 2001 From: Gabriele Antolini Date: Wed, 6 Nov 2024 15:51:21 +0100 Subject: [PATCH 25/34] derived variables meteogrid --- meteo/meteoGrid.cpp | 10 ++--- meteo/meteoGrid.h | 4 +- meteo/meteoPoint.cpp | 75 ++++++++++++++++------------------ meteo/meteoPoint.h | 4 +- pragaProject/pragaProject.cpp | 76 +++++++++++++++++++++++++++++++---- pragaProject/pragaProject.h | 7 ++-- pragaProject/pragaShell.cpp | 21 ++++++++-- 7 files changed, 134 insertions(+), 63 deletions(-) diff --git a/meteo/meteoGrid.cpp b/meteo/meteoGrid.cpp index 4a28ac979..a1f36e302 100644 --- a/meteo/meteoGrid.cpp +++ b/meteo/meteoGrid.cpp @@ -856,7 +856,7 @@ void Crit3DMeteoGrid::fixDailyThermalConsistency(const Crit3DDate myDate) { if (tmin > tmax) { - _meteoPoints[row][col]->setMeteoPointValueD(myDate, dailyAirTemperatureMin, tmax - 0.1); + _meteoPoints[row][col]->setMeteoPointValueD(myDate, dailyAirTemperatureMin, float(tmax - 0.1)); } } } @@ -1059,7 +1059,7 @@ void Crit3DMeteoGrid::saveRowColfromZone(gis::Crit3DRasterGrid* zoneGrid, std::v } -void Crit3DMeteoGrid::computeHourlyDerivedVariables(Crit3DTime dateTime, std::vector variables, bool useNetRad) +void Crit3DMeteoGrid::computeHourlyDerivedVar(Crit3DTime dateTime, meteoVariable myVar, bool useNetRad) { for (unsigned row = 0; row < unsigned(gridStructure().header().nrRows); row++) @@ -1068,13 +1068,13 @@ void Crit3DMeteoGrid::computeHourlyDerivedVariables(Crit3DTime dateTime, std::ve { if (_meteoPoints[row][col]->active) { - _meteoPoints[row][col]->computeHourlyDerivedVariables(dateTime, variables, useNetRad); + _meteoPoints[row][col]->computeHourlyDerivedVar(dateTime, myVar, useNetRad); } } } } -void Crit3DMeteoGrid::computeDailyDerivedVariables(Crit3DDate date, std::vector variables, Crit3DMeteoSettings& meteoSettings) +void Crit3DMeteoGrid::computeDailyDerivedVar(Crit3DDate date, meteoVariable myVar, Crit3DMeteoSettings& meteoSettings) { for (unsigned row = 0; row < unsigned(gridStructure().header().nrRows); row++) @@ -1083,7 +1083,7 @@ void Crit3DMeteoGrid::computeDailyDerivedVariables(Crit3DDate date, std::vector< { if (_meteoPoints[row][col]->active) { - _meteoPoints[row][col]->computeDailyDerivedVariables(date, variables, meteoSettings); + _meteoPoints[row][col]->computeDailyDerivedVar(date, myVar, meteoSettings); } } } diff --git a/meteo/meteoGrid.h b/meteo/meteoGrid.h index da4179b50..6b4a88de3 100644 --- a/meteo/meteoGrid.h +++ b/meteo/meteoGrid.h @@ -150,8 +150,8 @@ void computeRelativeHumidityFromTd(const Crit3DDate myDate, const int myHour); void computeWindVectorHourly(const Crit3DDate myDate, const int myHour); void fixDailyThermalConsistency(const Crit3DDate myDate); - void computeHourlyDerivedVariables(Crit3DTime dateTime, std::vector variables, bool useNetRad); - void computeDailyDerivedVariables(Crit3DDate date, std::vector variables, Crit3DMeteoSettings &meteoSettings); + void computeHourlyDerivedVar(Crit3DTime dateTime, meteoVariable myVar, bool useNetRad); + void computeDailyDerivedVar(Crit3DDate date, meteoVariable myVar, Crit3DMeteoSettings &meteoSettings); private: diff --git a/meteo/meteoPoint.cpp b/meteo/meteoPoint.cpp index 5fe71f248..62327dbac 100644 --- a/meteo/meteoPoint.cpp +++ b/meteo/meteoPoint.cpp @@ -1218,63 +1218,58 @@ std::vector Crit3DMeteoPoint::getProxyValues() return myValues; } -bool Crit3DMeteoPoint::computeHourlyDerivedVariables(Crit3DTime dateTime, std::vector variables, bool useNetRad) +bool Crit3DMeteoPoint::computeHourlyDerivedVar(Crit3DTime dateTime, meteoVariable myVar, bool useNetRad) { Crit3DDate myDate = dateTime.date; int myHour = dateTime.getHour(); + float value = NODATA; + short valueShort = NODATA; - for (auto var : variables) + if (myVar == leafWetness) { - if (var == leafWetness) + if (computeLeafWetness(getMeteoPointValueH(myDate, myHour, 0, precipitation), + getMeteoPointValueH(myDate, myHour, 0, airRelHumidity), &valueShort)) + setMeteoPointValueH(myDate, myHour, 0, leafWetness, float(valueShort)); + } + else if (myVar == referenceEvapotranspiration) + { + if (useNetRad) { - short leafW; - if (computeLeafWetness(getMeteoPointValueH(myDate, myHour, 0, precipitation), - getMeteoPointValueH(myDate, myHour, 0, airRelHumidity), &leafW)) - setMeteoPointValueH(myDate, myHour, 0, leafWetness, leafW); + value = float(ET0_Penman_hourly_net_rad(double(point.z), + double(getMeteoPointValueH(myDate, myHour, 0, netIrradiance)), + double(getMeteoPointValueH(myDate, myHour, 0, airTemperature)), + double(getMeteoPointValueH(myDate, myHour, 0, airRelHumidity)), + double(getMeteoPointValueH(myDate, myHour, 0, windScalarIntensity)))); + } - else if (var == referenceEvapotranspiration) + else { - float et0; - - if (useNetRad) - { - et0 = float(ET0_Penman_hourly_net_rad(double(point.z), - double(getMeteoPointValueH(myDate, myHour, 0, netIrradiance)), - double(getMeteoPointValueH(myDate, myHour, 0, airTemperature)), - double(getMeteoPointValueH(myDate, myHour, 0, airRelHumidity)), - double(getMeteoPointValueH(myDate, myHour, 0, windScalarIntensity)))); - - } - else - { - //da sistemare - et0 = ET0_Penman_hourly(double(point.z), - double(getMeteoPointValueH(myDate, myHour, 0, atmTransmissivity) / 0.75), - double(getMeteoPointValueH(myDate, myHour, 0, globalIrradiance)), - double(getMeteoPointValueH(myDate, myHour, 0, airTemperature)), - double(getMeteoPointValueH(myDate, myHour, 0, airRelHumidity)), - double(getMeteoPointValueH(myDate, myHour, 0, windScalarIntensity))); - } - - setMeteoPointValueH(myDate, myHour, 0, referenceEvapotranspiration, et0); + //da sistemare + value = ET0_Penman_hourly(double(point.z), + double(getMeteoPointValueH(myDate, myHour, 0, atmTransmissivity) / float(0.75)), + double(getMeteoPointValueH(myDate, myHour, 0, globalIrradiance)), + double(getMeteoPointValueH(myDate, myHour, 0, airTemperature)), + double(getMeteoPointValueH(myDate, myHour, 0, airRelHumidity)), + double(getMeteoPointValueH(myDate, myHour, 0, windScalarIntensity))); } } + setMeteoPointValueH(myDate, myHour, 0, myVar, value); + return true; } -bool Crit3DMeteoPoint::computeDailyDerivedVariables(Crit3DDate date, std::vector variables, Crit3DMeteoSettings& meteoSettings) +bool Crit3DMeteoPoint::computeDailyDerivedVar(Crit3DDate date, meteoVariable myVar, Crit3DMeteoSettings& meteoSettings) { - for (auto var : variables) + float value; + + if (myVar == dailyReferenceEvapotranspirationHS) { - if (var == dailyReferenceEvapotranspirationHS) - { - float et0 = dailyEtpHargreaves(getMeteoPointValueD(date, dailyAirTemperatureMin), - getMeteoPointValueD(date, dailyAirTemperatureMax), - date, latitude, &meteoSettings); + value = dailyEtpHargreaves(getMeteoPointValueD(date, dailyAirTemperatureMin), + getMeteoPointValueD(date, dailyAirTemperatureMax), + date, latitude, &meteoSettings); - setMeteoPointValueD(date, dailyReferenceEvapotranspirationHS, et0); - } + setMeteoPointValueD(date, dailyReferenceEvapotranspirationHS, value); } return true; diff --git a/meteo/meteoPoint.h b/meteo/meteoPoint.h index f8480069e..f8edf6adc 100644 --- a/meteo/meteoPoint.h +++ b/meteo/meteoPoint.h @@ -171,8 +171,8 @@ void setAltitude(double altitude); void setLapseRateCode(std::string lapseRateCode); - bool computeHourlyDerivedVariables(Crit3DTime dateTime, std::vector variables, bool useNetRad); - bool computeDailyDerivedVariables(Crit3DDate date, std::vector variables, Crit3DMeteoSettings &meteoSettings); + bool computeHourlyDerivedVar(Crit3DTime dateTime, meteoVariable myVar, bool useNetRad); + bool computeDailyDerivedVar(Crit3DDate date, meteoVariable myVar, Crit3DMeteoSettings &meteoSettings); bool computeMonthlyAggregate(Crit3DDate firstDate, Crit3DDate lastDate, meteoVariable dailyMeteoVar, Crit3DMeteoSettings *meteoSettings, Crit3DQuality *qualityCheck, Crit3DClimateParameters *climateParam); TObsDataH *getObsDataH() const; void initializeObsDataDFromMp(unsigned int numberOfDays, const Crit3DDate& firstDate, Crit3DMeteoPoint mp); diff --git a/pragaProject/pragaProject.cpp b/pragaProject/pragaProject.cpp index 944847bf1..ac860c5da 100644 --- a/pragaProject/pragaProject.cpp +++ b/pragaProject/pragaProject.cpp @@ -2241,7 +2241,7 @@ bool PragaProject::timeAggregateGrid(QDate dateIni, QDate dateFin, QList variables, bool useNetRad) +bool PragaProject::derivedVariablesMeteoGridPeriod(QDate first, QDate last, QList variables, bool useNetRad) { // check meteo grid @@ -2299,7 +2299,9 @@ bool PragaProject::derivedVariablesGrid(QDate first, QDate last, QList meteoGrid()->computeHourlyDerivedVariables(getCrit3DTime(firstDateTime), hourlyVars, useNetRad); + foreach (myVar, hourlyVars) + meteoGridDbHandler->meteoGrid()->computeHourlyDerivedVar(getCrit3DTime(firstDateTime), myVar, useNetRad); + firstDateTime = firstDateTime.addSecs(3600); } } @@ -2309,14 +2311,14 @@ bool PragaProject::derivedVariablesGrid(QDate first, QDate last, QList meteoGrid()->computeDailyDerivedVariables(getCrit3DDate(myDate), dailyVars, *meteoSettings); + foreach (myVar, dailyVars) + meteoGridDbHandler->meteoGrid()->computeDailyDerivedVar(getCrit3DDate(myDate), myVar, *meteoSettings); + myDate = myDate.addDays(1); } } firstDateTime = QDateTime(first, QTime(1,0), Qt::UTC); - // save derived variables - variables << leafWetness << referenceEvapotranspiration; QString myError; logInfoGUI("Saving meteo grid data"); if (! meteoGridDbHandler->saveGridData(&myError, firstDateTime, lastDateTime, variables, meteoSettings)) return false; @@ -2567,6 +2569,34 @@ bool PragaProject::interpolationOutputPointsPeriod(QDate firstDate, QDate lastDa return true; } +bool PragaProject::deriveVariableMeteoGrid(meteoVariable myVar, frequencyType myFrequency, const Crit3DTime& myTime) +{ + if (interpolationSettings.getMeteoGridUpscaleFromDem()) + { + if (myVar == leafWetness) + { + if (! hourlyMeteoMaps->computeLeafWetnessMap()) return false; + meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, myFrequency, myTime.date, myTime.getHour(), myTime.getMinutes(), + &DEM, hourlyMeteoMaps->mapHourlyLeafW, interpolationSettings.getMeteoGridAggrMethod()); + } + else if (myVar == referenceEvapotranspiration) + { + if (! hourlyMeteoMaps->computeET0PMMap(DEM, radiationMaps)) return false; + meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, myFrequency, myTime.date, myTime.getHour(), myTime.getMinutes(), + &DEM, hourlyMeteoMaps->mapHourlyET0, interpolationSettings.getMeteoGridAggrMethod()); + } + } + else + { + if (myFrequency == hourly) + meteoGridDbHandler->meteoGrid()->computeHourlyDerivedVar(myTime, myVar, false); + else if (myFrequency == daily) + meteoGridDbHandler->meteoGrid()->computeDailyDerivedVar(myTime.date, myVar, *meteoSettings); + } + + return true; +} + bool PragaProject::interpolationMeteoGrid(meteoVariable myVar, frequencyType myFrequency, const Crit3DTime& myTime) { if (meteoGridDbHandler == nullptr) @@ -2603,7 +2633,6 @@ bool PragaProject::interpolationMeteoGrid(meteoVariable myVar, frequencyType myF hourlyMeteoMaps->computeLeafWetnessMap() ; meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, myFrequency, myTime.date, myTime.getHour(), myTime.getMinutes(), &DEM, hourlyMeteoMaps->mapHourlyLeafW, interpolationSettings.getMeteoGridAggrMethod()); - } else if (myVar == referenceEvapotranspiration) { @@ -2655,7 +2684,9 @@ bool PragaProject::interpolationMeteoGrid(meteoVariable myVar, frequencyType myF { // for radiation use dem and aggregate if (! interpolateDemRadiation(myTime.addSeconds(-1800), radiationMaps->globalRadiationMap)) return false; - meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(globalIrradiance, daily, myTime.date, myTime.getHour(), myTime.getMinutes(), + meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(atmTransmissivity, hourly, myTime.date, myTime.getHour(), myTime.getMinutes(), + &DEM, radiationMaps->globalRadiationMap, interpolationSettings.getMeteoGridAggrMethod()); + meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(globalIrradiance, hourly, myTime.date, myTime.getHour(), myTime.getMinutes(), &DEM, radiationMaps->globalRadiationMap, interpolationSettings.getMeteoGridAggrMethod()); } else @@ -2678,7 +2709,7 @@ bool PragaProject::interpolationMeteoGrid(meteoVariable myVar, frequencyType myF } bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QList variables, - QList aggrVariables, + QList derivedVariables, QList aggrVariables, int nrDaysLoading, int nrDaysSaving) { // check variables @@ -2749,6 +2780,23 @@ bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QL varToSave.push_back(windVectorIntensity); } + // derivedVariables + std::vector hourlyDerivedVars, dailyDerivedVars; + foreach (myVar, derivedVariables) + { + freq = getVarFrequency(myVar); + + if (freq == noFrequency) + { + errorString = "Unknown variable: " + QString::fromStdString(getMeteoVarName(myVar)); + return false; + } + else if (freq == hourly) + hourlyDerivedVars.push_back(myVar); + else if (freq == daily) + dailyDerivedVars.push_back(myVar); + } + // find out if detrending needed bool useProxies = false; foreach (myVar, variables) @@ -2849,6 +2897,12 @@ bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QL if (! interpolationMeteoGrid(myVar, hourly, getCrit3DTime(myDate, myHour))) return false; } } + + foreach (myVar, hourlyDerivedVars) + { + logInfo(QString::fromStdString(getMeteoVarName(myVar))); + deriveVariableMeteoGrid(myVar, hourly, getCrit3DTime(myDate, myHour)); + } } } @@ -2864,6 +2918,12 @@ bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QL if (! interpolationMeteoGrid(myVar, daily, getCrit3DTime(myDate, myHour))) return false; } } + + foreach (myVar, dailyDerivedVars) + { + logInfo(QString::fromStdString(getMeteoVarName(myVar))); + deriveVariableMeteoGrid(myVar, daily, getCrit3DTime(myDate, myHour)); + } } if (countDaysSaving == nrDaysSaving || myDate == dateFin || myDate == loadDateFin) diff --git a/pragaProject/pragaProject.h b/pragaProject/pragaProject.h index cea938f6d..fbf318147 100644 --- a/pragaProject/pragaProject.h +++ b/pragaProject/pragaProject.h @@ -109,15 +109,16 @@ bool interpolationOutputPointsPeriod(QDate dateIni, QDate lastDate, QList variables); + bool deriveVariableMeteoGrid(meteoVariable myVar, frequencyType myFrequency, const Crit3DTime& myTime); bool interpolationMeteoGrid(meteoVariable myVar, frequencyType myFrequency, const Crit3DTime& myTime); - bool interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QList variables, - QList aggrVariables, int nrDaysLoading, int nrDaysSaving); + bool interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QList variables, QList aggrVariables, + QList derivedVariables, int nrDaysLoading, int nrDaysSaving); bool interpolationCrossValidationPeriod(QDate dateIni, QDate dateFin, meteoVariable myVar, QString filename); bool saveGrid(meteoVariable myVar, frequencyType myFrequency, const Crit3DTime& myTime, bool showInfo); bool timeAggregateGridVarHourlyInDaily(meteoVariable dailyVar, Crit3DDate dateIni, Crit3DDate dateFin); bool timeAggregateGrid(QDate dateIni, QDate dateFin, QList variables, bool loadData, bool saveData); bool computeDailyVariablesPoint(Crit3DMeteoPoint *meteoPoint, QDate first, QDate last, QList variables); - bool derivedVariablesGrid(QDate first, QDate last, QList variables, bool useNetRad); + bool derivedVariablesMeteoGridPeriod(QDate first, QDate last, QList variables, bool useNetRad); bool elaborationPointsCycle(bool isAnomaly, bool showInfo); bool elaborationPointsCycleGrid(bool isAnomaly, bool showInfo); bool elaborationCheck(bool isMeteoGrid, bool isAnomaly); diff --git a/pragaProject/pragaShell.cpp b/pragaProject/pragaShell.cpp index 36245d408..8f793543e 100644 --- a/pragaProject/pragaShell.cpp +++ b/pragaProject/pragaShell.cpp @@ -287,7 +287,7 @@ int cmdInterpolationGridPeriod(PragaProject* myProject, QList argumentL QDate dateIni, dateFin; QList varString; - QList variables, aggrVariables; + QList variables, aggrVariables, derivedVariables; QString var; meteoVariable meteoVar; int saveInterval = 1; @@ -311,6 +311,21 @@ int cmdInterpolationGridPeriod(PragaProject* myProject, QList argumentL variables << meteoVar; } } + else if (argumentList[i].left(3) == "-d:") + { + varString = argumentList[i].right(argumentList[i].length()-3).split(","); + foreach (var,varString) + { + meteoVar = getMeteoVar(var.toStdString()); + + if (meteoVar == noMeteoVar) { + myProject->errorString = "Unknown variable: " + var; + return PRAGA_INVALID_COMMAND; + } + + derivedVariables << meteoVar; + } + } else if (argumentList[i].left(3) == "-a:") { varString = argumentList[i].right(argumentList[i].length()-3).split(","); @@ -371,7 +386,7 @@ int cmdInterpolationGridPeriod(PragaProject* myProject, QList argumentL return PRAGA_INVALID_COMMAND; } - if (! myProject->interpolationMeteoGridPeriod(dateIni, dateFin, variables, aggrVariables, loadInterval, saveInterval)) + if (! myProject->interpolationMeteoGridPeriod(dateIni, dateFin, variables, derivedVariables, aggrVariables, loadInterval, saveInterval)) return PRAGA_ERROR; return PRAGA_OK; @@ -488,7 +503,7 @@ int cmdHourlyDerivedVariablesGrid(PragaProject* myProject, QList argume return PRAGA_INVALID_COMMAND; } - if (! myProject->derivedVariablesGrid(first, last, variables, useNetRad)) + if (! myProject->derivedVariablesMeteoGridPeriod(first, last, variables, useNetRad)) return PRAGA_ERROR; return PRAGA_OK; From 3a51861887b1c8669520129a5432379aaba67e10 Mon Sep 17 00:00:00 2001 From: Gabriele Antolini Date: Thu, 7 Nov 2024 10:17:17 +0100 Subject: [PATCH 26/34] derived variables meteogrid fixed --- pragaProject/pragaProject.cpp | 25 +++++++------------------ 1 file changed, 7 insertions(+), 18 deletions(-) diff --git a/pragaProject/pragaProject.cpp b/pragaProject/pragaProject.cpp index ac860c5da..4ecd9cc64 100644 --- a/pragaProject/pragaProject.cpp +++ b/pragaProject/pragaProject.cpp @@ -2585,6 +2585,12 @@ bool PragaProject::deriveVariableMeteoGrid(meteoVariable myVar, frequencyType my meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, myFrequency, myTime.date, myTime.getHour(), myTime.getMinutes(), &DEM, hourlyMeteoMaps->mapHourlyET0, interpolationSettings.getMeteoGridAggrMethod()); } + else if (myVar == dailyReferenceEvapotranspirationHS) + { + if (! pragaDailyMaps->computeHSET0Map(&gisSettings, myTime.date)) return false; + meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, myFrequency, myTime.date, myTime.getHour(), myTime.getMinutes(), + &DEM, pragaDailyMaps->mapDailyET0HS, interpolationSettings.getMeteoGridAggrMethod()); + } } else { @@ -2628,18 +2634,6 @@ bool PragaProject::interpolationMeteoGrid(meteoVariable myVar, frequencyType myF &DEM, getPragaMapFromVar(windVectorY), interpolationSettings.getMeteoGridAggrMethod()); meteoGridDbHandler->meteoGrid()->computeWindVectorHourly(myTime.date, myTime.getHour()); } - else if (myVar == leafWetness) - { - hourlyMeteoMaps->computeLeafWetnessMap() ; - meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, myFrequency, myTime.date, myTime.getHour(), myTime.getMinutes(), - &DEM, hourlyMeteoMaps->mapHourlyLeafW, interpolationSettings.getMeteoGridAggrMethod()); - } - else if (myVar == referenceEvapotranspiration) - { - hourlyMeteoMaps->computeET0PMMap(DEM, radiationMaps); - meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, myFrequency, myTime.date, myTime.getHour(), myTime.getMinutes(), - &DEM, hourlyMeteoMaps->mapHourlyET0, interpolationSettings.getMeteoGridAggrMethod()); - } else { if (!interpolationDemMain(myVar, myTime, getPragaMapFromVar(myVar))) return false; @@ -2649,13 +2643,8 @@ bool PragaProject::interpolationMeteoGrid(meteoVariable myVar, frequencyType myF } else if (myFrequency == daily) { - if (myVar == dailyReferenceEvapotranspirationHS) { - pragaDailyMaps->computeHSET0Map(&gisSettings, myTime.date);} - else - { - if (! interpolationDemMain(myVar, myTime, getPragaMapFromVar(myVar))) + if (! interpolationDemMain(myVar, myTime, getPragaMapFromVar(myVar))) return false; - } if (myVar == dailyAirTemperatureMax || myVar == dailyAirTemperatureMin) { if (! pragaDailyMaps->fixDailyThermalConsistency()) return false;} From a23519ca258a86bd1d2dd8823358f63b10f3080e Mon Sep 17 00:00:00 2001 From: cate-to Date: Thu, 7 Nov 2024 10:45:04 +0100 Subject: [PATCH 27/34] cleanup --- interpolation/interpolation.cpp | 16 ++++++++-------- mathFunctions/furtherMathFunctions.cpp | 6 ++---- mathFunctions/furtherMathFunctions.h | 8 ++------ 3 files changed, 12 insertions(+), 18 deletions(-) diff --git a/interpolation/interpolation.cpp b/interpolation/interpolation.cpp index 60cf172a5..0b2360b1b 100644 --- a/interpolation/interpolation.cpp +++ b/interpolation/interpolation.cpp @@ -1093,12 +1093,12 @@ void localSelection(vector &inputPoints, vector < for (unsigned long 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)); - unsigned int nrValid = 0; + int nrValid = 0; float stepRadius = 7500; // [m] float r0 = 0; // [m] float r1 = stepRadius; // [m] unsigned int i; - unsigned int nrPrimaries = 0; + int nrPrimaries = 0; int maxDistance = 0; while ((!mySettings.getUseLapseRateCode() && nrValid < minPoints) || (mySettings.getUseLapseRateCode() && nrPrimaries < minPoints)) @@ -1514,7 +1514,7 @@ void calculateFirstGuessCombinations(Crit3DProxy* myProxy) std::vector tempFirstGuess; int numSteps = 15; std::vector stepSize; - int nrParam = int(tempParam.size()/2); + unsigned nrParam = int(tempParam.size()/2); double min_,max_; for (unsigned j=0; j < nrParam; j++) @@ -1801,13 +1801,13 @@ bool multipleDetrendingElevationFitting(int elevationPos, std::vector &)>()), 400, 4, parametersMin, parametersMax, parameters, parametersDelta, - stepSize, numSteps, 1000, 0.002, 0.005, predictors, predictands, weights,firstGuessCombinations); + R2 = interpolation::bestFittingMarquardt_nDimension_singleFunction(*(myFunc.target&)>()), 4, parametersMin, parametersMax, parameters, parametersDelta, + 1000, 0.002, 0.005, predictors, predictands, weights,firstGuessCombinations); else - R2 = interpolation::bestFittingMarquardt_nDimension_singleFunction(*(myFunc.target&)>()), 400, 4, parametersMin, parametersMax, parameters, parametersDelta, - stepSize, numSteps, 1000, 0.002, 0.005, predictors, predictands,firstGuessCombinations); + R2 = interpolation::bestFittingMarquardt_nDimension_singleFunction(*(myFunc.target&)>()), 4, parametersMin, parametersMax, parameters, parametersDelta, + 1000, 0.002, 0.005, predictors, predictands,firstGuessCombinations); - mySettings->getProxy(elevationPos)->setRegressionR2(R2); + mySettings->getProxy(elevationPos)->setRegressionR2(float(R2)); std::vector> newParameters; newParameters.push_back(parameters); diff --git a/mathFunctions/furtherMathFunctions.cpp b/mathFunctions/furtherMathFunctions.cpp index bb1456450..bd13c1d19 100644 --- a/mathFunctions/furtherMathFunctions.cpp +++ b/mathFunctions/furtherMathFunctions.cpp @@ -2477,10 +2477,9 @@ namespace interpolation return norm; } double bestFittingMarquardt_nDimension_singleFunction(double (*func)(double, std::vector&), - int nrTrials, int nrMinima, + int nrMinima, std::vector & parametersMin, std::vector & parametersMax, std::vector & parameters, std::vector & parametersDelta, - std::vector & stepSize, int numSteps, int maxIterationsNr, double myEpsilon, double deltaR2, std::vector & x ,std::vector& y, std::vector& weights, std::vector> firstGuessCombinations) @@ -2547,10 +2546,9 @@ namespace interpolation } double bestFittingMarquardt_nDimension_singleFunction(double (*func)(double, std::vector&), - int nrTrials, int nrMinima, + int nrMinima, std::vector & parametersMin, std::vector & parametersMax, std::vector & parameters, std::vector & parametersDelta, - std::vector & stepSize, int numSteps, int maxIterationsNr, double myEpsilon, double deltaR2, std::vector & x ,std::vector& y, std::vector> firstGuessCombinations) diff --git a/mathFunctions/furtherMathFunctions.h b/mathFunctions/furtherMathFunctions.h index 9092fa2e1..9ed70cbfc 100644 --- a/mathFunctions/furtherMathFunctions.h +++ b/mathFunctions/furtherMathFunctions.h @@ -165,19 +165,15 @@ enum estimatedFunction {FUNCTION_CODE_SPHERICAL, FUNCTION_CODE_LINEAR, FUNCTION_ std::vector ¶meters, std::vector ¶metersDelta, int maxIterationsNr, double myEpsilon, std::vector & x, std::vector& y); - double bestFittingMarquardt_nDimension_singleFunction(double (*func)(double, std::vector&), - int nrTrials, int nrMinima, + double bestFittingMarquardt_nDimension_singleFunction(double (*func)(double, std::vector&), int nrMinima, std::vector & parametersMin, std::vector & parametersMax, std::vector & parameters, std::vector & parametersDelta, - std::vector &stepSize, int numSteps, int maxIterationsNr, double myEpsilon, double deltaR2, std::vector & x , std::vector& y, std::vector& weights, std::vector > firstGuessCombinations); - double bestFittingMarquardt_nDimension_singleFunction(double (*func)(double, std::vector&), - int nrTrials, int nrMinima, + double bestFittingMarquardt_nDimension_singleFunction(double (*func)(double, std::vector&), int nrMinima, std::vector & parametersMin, std::vector & parametersMax, std::vector & parameters, std::vector & parametersDelta, - std::vector &stepSize, int numSteps, int maxIterationsNr, double myEpsilon, double deltaR2, std::vector & x , std::vector& y, std::vector > firstGuessCombinations); From 8d4a42ec998fd127f36d84c74df70e507f7354b7 Mon Sep 17 00:00:00 2001 From: avolta Date: Fri, 8 Nov 2024 16:44:35 +0100 Subject: [PATCH 28/34] introduzione del metodo della sezione aurea nel calcolo della TAD --- interpolation/interpolation.cpp | 70 ++++++++++++++++++++++++++++++--- interpolation/interpolation.h | 11 ++++++ mathFunctions/commonConstants.h | 1 + 3 files changed, 77 insertions(+), 5 deletions(-) diff --git a/interpolation/interpolation.cpp b/interpolation/interpolation.cpp index 60cf172a5..e217e0eb8 100644 --- a/interpolation/interpolation.cpp +++ b/interpolation/interpolation.cpp @@ -42,7 +42,7 @@ #include "meteo.h" #include - +#include using namespace std; @@ -2129,6 +2129,56 @@ bool glocalDetrendingFitting(std::vector &myPoint return true; } +double goldenSectionSearch(meteoVariable myVar,Crit3DMeteoPoint* &myMeteoPoints,int nrMeteoPoints, + std::vector &interpolationPoints,Crit3DInterpolationSettings* mySettings, + Crit3DMeteoSettings* meteoSettings, double a, double b) +{ + // this function finds the minimum by means of golden section method + double tol = 1; + //const double phi = (1 + std::sqrt(5)) / 2; // golden section + double x1 = b - (b - a) / GOLDEN_SECTION; + double x2 = a + (b - a) / GOLDEN_SECTION; + int counter=0; + while (std::abs(b - a) > tol && counter<100) + { + counter++; + if (topographicDistanceInternalFunction(myVar,myMeteoPoints,nrMeteoPoints,interpolationPoints, mySettings,meteoSettings, x1) < + topographicDistanceInternalFunction(myVar,myMeteoPoints,nrMeteoPoints,interpolationPoints, mySettings,meteoSettings,x2)) + { + b = x2; + x2 = x1; + x1 = b - (b - a) / GOLDEN_SECTION; + mySettings->addToKhSeries(float(x1), topographicDistanceInternalFunction(myVar,myMeteoPoints,nrMeteoPoints,interpolationPoints, mySettings,meteoSettings, x1)); + } + else + { + a = x1; + x1 = x2; + x2 = a + (b - a) / GOLDEN_SECTION; + mySettings->addToKhSeries(float(x2), topographicDistanceInternalFunction(myVar,myMeteoPoints,nrMeteoPoints,interpolationPoints, mySettings,meteoSettings, x2)); + } + } + mySettings->addToKhSeries(float((a + b) / 2), topographicDistanceInternalFunction(myVar,myMeteoPoints,nrMeteoPoints,interpolationPoints, mySettings,meteoSettings, (a + b) / 2)); + return (a + b) / 2; // approximated minimum +} + + +double topographicDistanceInternalFunction(meteoVariable myVar, + Crit3DMeteoPoint* &myMeteoPoints, + int nrMeteoPoints, + std::vector &interpolationPoints, + Crit3DInterpolationSettings* mySettings, + Crit3DMeteoSettings* meteoSettings, double khFloat) +{ + float avgError = 0; + int kh = int(khFloat); + mySettings->setTopoDist_Kh(kh); + if (computeResiduals(myVar, myMeteoPoints, nrMeteoPoints, interpolationPoints, mySettings, meteoSettings, true, true)) + { + avgError = computeErrorCrossValidation(myMeteoPoints, nrMeteoPoints); + } + return avgError; +} void topographicDistanceOptimize(meteoVariable myVar, Crit3DMeteoPoint* &myMeteoPoints, @@ -2137,19 +2187,28 @@ void topographicDistanceOptimize(meteoVariable myVar, Crit3DInterpolationSettings* mySettings, Crit3DMeteoSettings* meteoSettings) { - float avgError; + //float avgError; mySettings->initializeKhSeries(); - int kh = 0; - int bestKh = kh; - float bestError = NODATA; + //float kh = 0; + double bestKh = 0; + //float bestError = NODATA; + std::vector avgErrorVec; + float khMin = 0; + float khMax = mySettings->getTopoDist_maxKh(); + //kh = 0.5*(khMin+khMax); + bestKh = goldenSectionSearch(myVar, myMeteoPoints, nrMeteoPoints,interpolationPoints,mySettings,meteoSettings, khMin,khMax); + mySettings->setTopoDist_Kh(int(bestKh)); + /* while (kh <= mySettings->getTopoDist_maxKh()) { mySettings->setTopoDist_Kh(kh); if (computeResiduals(myVar, myMeteoPoints, nrMeteoPoints, interpolationPoints, mySettings, meteoSettings, true, true)) { avgError = computeErrorCrossValidation(myMeteoPoints, nrMeteoPoints); + avgErrorVec.push_back(avgError); + if (isEqual(bestError, NODATA) || avgError < bestError) { bestError = avgError; @@ -2162,6 +2221,7 @@ void topographicDistanceOptimize(meteoVariable myVar, } mySettings->setTopoDist_Kh(bestKh); + */ } diff --git a/interpolation/interpolation.h b/interpolation/interpolation.h index f145ccb33..78a572917 100644 --- a/interpolation/interpolation.h +++ b/interpolation/interpolation.h @@ -51,6 +51,17 @@ void topographicDistanceOptimize(meteoVariable myVar, Crit3DMeteoPoint* &myMeteoPoints, int nrMeteoPoints, std::vector &interpolationPoints, Crit3DInterpolationSettings* mySettings, Crit3DMeteoSettings* meteoSettings); + double topographicDistanceInternalFunction(meteoVariable myVar, + Crit3DMeteoPoint* &myMeteoPoints, + int nrMeteoPoints, + std::vector &interpolationPoints, + Crit3DInterpolationSettings* mySettings, + Crit3DMeteoSettings* meteoSettings, double khFloat); + + double goldenSectionSearch(meteoVariable myVar,Crit3DMeteoPoint* &myMeteoPoints,int nrMeteoPoints, + std::vector &interpolationPoints,Crit3DInterpolationSettings* mySettings, + Crit3DMeteoSettings* meteoSettings, double a, double b); + bool krigingEstimateVariogram(float *myDist, float *mySemiVar,int sizeMyVar, int nrMyPoints,float myMaxDistance, double *mySill, double *myNugget, double *myRange, double *mySlope, TkrigingMode *myMode, int nrPointData); bool krigLinearPrep(double *mySlope, double *myNugget, int nrPointData); diff --git a/mathFunctions/commonConstants.h b/mathFunctions/commonConstants.h index f507f7b06..8c3efbf3d 100644 --- a/mathFunctions/commonConstants.h +++ b/mathFunctions/commonConstants.h @@ -227,6 +227,7 @@ #define DEG_TO_RAD 0.0174532925 #define RAD_TO_DEG 57.2957795 #define SQRT_2 1.41421356237 + #define GOLDEN_SECTION 1.6180339887498948482 #define MINIMUM_PERCENTILE_DATA 3 From 754a2e9038d178f307f8613b2e9d4d4461525413 Mon Sep 17 00:00:00 2001 From: gantolini Date: Tue, 12 Nov 2024 12:34:16 +0100 Subject: [PATCH 29/34] glocal files loading --- interpolation/interpolationSettings.cpp | 6 ++++++ interpolation/interpolationSettings.h | 1 + pragaProject/pragaProject.cpp | 13 ++++++++++--- project/project.cpp | 23 +++++++++++------------ 4 files changed, 28 insertions(+), 15 deletions(-) diff --git a/interpolation/interpolationSettings.cpp b/interpolation/interpolationSettings.cpp index 8618b983b..e017b7142 100644 --- a/interpolation/interpolationSettings.cpp +++ b/interpolation/interpolationSettings.cpp @@ -475,6 +475,7 @@ void Crit3DInterpolationSettings::initialize() initializeProxy(); } + std::string getKeyStringInterpolationMethod(TInterpolationMethod value) { std::map::const_iterator it; @@ -598,6 +599,11 @@ int Crit3DInterpolationSettings::getProxyPosFromName(TProxyVar name) return NODATA; } +bool Crit3DInterpolationSettings::isGlocalReady() +{ + return (getMacroAreasMap() != nullptr && getMacroAreas().size() > 0); +} + std::string Crit3DProxy::getName() const { return name; diff --git a/interpolation/interpolationSettings.h b/interpolation/interpolationSettings.h index 6f8072af0..531f4d8c9 100644 --- a/interpolation/interpolationSettings.h +++ b/interpolation/interpolationSettings.h @@ -230,6 +230,7 @@ void setUseLocalDetrending(bool myValue); bool getUseLocalDetrending(); + bool isGlocalReady(); void setUseGlocalDetrending(bool myValue); bool getUseGlocalDetrending(); void setMacroAreasMap(gis::Crit3DRasterGrid *value); diff --git a/pragaProject/pragaProject.cpp b/pragaProject/pragaProject.cpp index 4ecd9cc64..6f9a120ec 100644 --- a/pragaProject/pragaProject.cpp +++ b/pragaProject/pragaProject.cpp @@ -2611,6 +2611,13 @@ bool PragaProject::interpolationMeteoGrid(meteoVariable myVar, frequencyType myF return false; } + // check glocal + if (interpolationSettings.getUseGlocalDetrending() && ! interpolationSettings.isGlocalReady()) + { + if (! loadGlocalAreasMap()) return false; + if (! loadGlocalStationsAndCells(!interpolationSettings.getMeteoGridUpscaleFromDem())) return false; + } + if (interpolationSettings.getMeteoGridUpscaleFromDem()) { if (myFrequency == hourly) @@ -2665,8 +2672,8 @@ bool PragaProject::interpolationMeteoGrid(meteoVariable myVar, frequencyType myF } else if (myVar == windVectorDirection || myVar == windVectorIntensity) { - if (! interpolationMeteoGrid(windVectorX, hourly, myTime)) return false; - if (! interpolationMeteoGrid(windVectorY, hourly, myTime)) return false; + if (! interpolationGrid(windVectorX, myTime)) return false; + if (! interpolationGrid(windVectorY, myTime)) return false; meteoGridDbHandler->meteoGrid()->computeWindVectorHourly(myTime.date, myTime.getHour()); } else if (myVar == globalIrradiance) @@ -2818,7 +2825,7 @@ bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QL return false; } - if (interpolationSettings.getUseGlocalDetrending()) + if (interpolationSettings.getUseGlocalDetrending() && ! interpolationSettings.isGlocalReady()) { if (! loadGlocalAreasMap()) return false; if (! loadGlocalStationsAndCells(!interpolationSettings.getMeteoGridUpscaleFromDem())) return false; diff --git a/project/project.cpp b/project/project.cpp index 4c7bb2721..407377f37 100644 --- a/project/project.cpp +++ b/project/project.cpp @@ -235,7 +235,7 @@ bool Project::checkProxy(Crit3DProxy &myProxy, QString* error, bool isActive) if (myProxy.getFittingParametersRange().size() != nrParameters*2) { - *error = "wrong numer of parameters for proxy: " + QString::fromStdString(name_); + *error = "wrong number of parameters for proxy: " + QString::fromStdString(name_); return false; } @@ -2288,9 +2288,9 @@ bool Project::loadGlocalStationsAndCells(bool isGrid) Crit3DMacroArea myArea; //per ogni area, rappresentata da righe di areaPoints, si fa ciclo sui meteopoints - for (int j = 0; j < areaPoints.size(); j++) + for (unsigned j = 0; j < areaPoints.size(); j++) { - for (int k = 0; k < areaPoints[j].size(); k++) + for (unsigned k = 0; k < areaPoints[j].size(); k++) { //controlla se id si trova nel vettore areaPoints[j] e salva l'indice del meteopoint for (int i=0; i < nrMeteoPoints; i++) @@ -2424,11 +2424,9 @@ bool Project::loadGlocalStationsCsv(QString fileName, std::vector temp; while(!myStream.atEnd()) { - nrLine++; line = myStream.readLine().split(','); for (int i = 1; i < line.size(); i++) temp.push_back(line[i].toInt()); @@ -3105,6 +3103,13 @@ bool Project::interpolationDemMain(meteoVariable myVar, const Crit3DTime& myTime if (! checkInterpolation(myVar)) return false; + // check glocal + if (interpolationSettings.getUseGlocalDetrending() && ! interpolationSettings.isGlocalReady()) + { + if (! loadGlocalAreasMap()) return false; + if (! loadGlocalStationsAndCells(!interpolationSettings.getMeteoGridUpscaleFromDem())) return false; + } + // solar radiation model if (myVar == globalIrradiance) { @@ -3169,12 +3174,6 @@ bool Project::interpolationGrid(meteoVariable myVar, const Crit3DTime& myTime) if (interpolationSettings.getUseMultipleDetrending()) interpolationSettings.clearFitting(); - if (interpolationSettings.getUseGlocalDetrending()) - { - if (! loadGlocalAreasMap()) return false; - if (! loadGlocalStationsAndCells(!interpolationSettings.getMeteoGridUpscaleFromDem())) return false; - } - // check quality and pass data to interpolation if (! checkAndPassDataToInterpolation(quality, myVar, meteoPoints, nrMeteoPoints, myTime, &qualityInterpolationSettings, &interpolationSettings, meteoSettings, &climateParameters, interpolationPoints, @@ -4391,7 +4390,7 @@ void Project::showLocalProxyGraph(gis::Crit3DGeoPoint myPoint) myZDEM = DEM.value[row][col]; } - if (interpolationSettings.getUseGlocalDetrending()) + if (interpolationSettings.getUseGlocalDetrending() && ! interpolationSettings.isGlocalReady()) { if (loadGlocalAreasMap()) loadGlocalStationsAndCells(false); From fc7c7ed3696797711b41a9737f07d0a6a4249d3a Mon Sep 17 00:00:00 2001 From: Gabriele Antolini Date: Tue, 12 Nov 2024 17:43:47 +0100 Subject: [PATCH 30/34] fixed glocal files --- project/project.cpp | 58 ++++++++++++++++----------------------------- project/project.h | 3 ++- 2 files changed, 23 insertions(+), 38 deletions(-) diff --git a/project/project.cpp b/project/project.cpp index 407377f37..b6268e11e 100644 --- a/project/project.cpp +++ b/project/project.cpp @@ -102,7 +102,7 @@ void Project::initializeProject() outputPointsFileName = ""; currentDbOutputFileName = ""; glocalMapName = ""; - + glocalPointsName = ""; meteoPointsLoaded = false; meteoGridLoaded = false; loadGridDataAtStart = false; @@ -633,7 +633,10 @@ bool Project::loadParameters(QString parametersFileName) interpolationSettings.setUseGlocalDetrending(parametersSettings->value("glocalDetrending").toBool()); if (parametersSettings->contains("glocalMapName")) - this->glocalMapName = parametersSettings->value("glocalMapName").toString(); + glocalMapName = parametersSettings->value("glocalMapName").toString(); + + if (parametersSettings->contains("glocalPointsName")) + glocalPointsName = parametersSettings->value("glocalPointsName").toString(); if (parametersSettings->contains("meteogrid_upscalefromdem")) interpolationSettings.setMeteoGridUpscaleFromDem(parametersSettings->value("meteogrid_upscalefromdem").toBool()); @@ -2238,34 +2241,22 @@ bool Project::loadGlocalAreasMap() { //TODO: add a check for the code values? - QString mapsFolder = defaultPath + PATH_GEO; - if (! QDir(mapsFolder).exists()) - { - //errore? - } - - /*QString mapsFolder = projectPath + PATH_GLOCAL; - if (! QDir(mapsFolder).exists()) - { - QDir().mkdir(mapsFolder); - }*/ - std::string myError; - gis::Crit3DRasterGrid* macroAreasGrid = new gis::Crit3DRasterGrid(); - std::string fileNameMap = mapsFolder.toStdString() + glocalMapName.toStdString() + "_map"; + QString fileNameMap = getCompleteFileName(glocalMapName, PATH_GEO); - if (!QFile::exists(QString::fromStdString(fileNameMap + ".flt"))) + if (!QFile::exists(fileNameMap + ".flt")) { - myError = "Couldn't find " + fileNameMap + " file in " + mapsFolder.toStdString(); - logError(QString::fromStdString(myError)); + errorString = "Could not find " + fileNameMap + " file"; return false; } - if (gis::readEsriGrid(fileNameMap, macroAreasGrid, myError)) + gis::Crit3DRasterGrid* macroAreasGrid = new gis::Crit3DRasterGrid(); + + if (gis::readEsriGrid(fileNameMap.toStdString(), macroAreasGrid, myError)) interpolationSettings.setMacroAreasMap(macroAreasGrid); else { - myError = "Error loading:\n" + fileNameMap; + errorString = "Error loading:\n" + fileNameMap + "\n" + QString::fromStdString(myError); return false; } @@ -2275,12 +2266,10 @@ bool Project::loadGlocalAreasMap() bool Project::loadGlocalStationsAndCells(bool isGrid) { //leggi csv aree - QString mapsFolder = defaultPath + PATH_GEO; - std::string fileNameStations = mapsFolder.toStdString() + glocalMapName.toStdString() + "_stations"; + QString fileNameStations = getCompleteFileName(glocalPointsName, PATH_GEO); std::vector> areaPoints; - std::string myError; - loadGlocalStationsCsv(QString::fromStdString(fileNameStations), areaPoints); + if (! loadGlocalStationsCsv(fileNameStations, areaPoints)) return false; //salva in vettore gli indici dei meteopoints che appartengono a area n, infine salva quel vettore std::vector temp; @@ -2310,7 +2299,7 @@ bool Project::loadGlocalStationsAndCells(bool isGrid) //assegnazione pesi a ogni cella di ogni area if (!loadGlocalWeightMaps(myAreas, isGrid)) { - myError = "error while loading the glocal weight maps."; + errorString = "Error loading glocal weight maps"; return false; } @@ -2389,28 +2378,23 @@ bool Project::loadGlocalWeightMaps(std::vector &myAreas, bool i bool Project::loadGlocalStationsCsv(QString fileName, std::vector> &areaPoints) { - std::string myError; - if (fileName == "") { - myError = "Missing csv filename"; - logError(QString::fromStdString(myError)); + errorString = "Missing csv filename"; return false; } areaPoints.clear(); - if (! QFile(fileName + ".csv").exists() || ! QFileInfo(fileName + ".csv").isFile()) + if (! QFile(fileName).exists() || ! QFileInfo(fileName).isFile()) { - myError = "Missing csv file: " + fileName.toStdString(); - logError(QString::fromStdString(myError)); + errorString = "Missing csv file: " + fileName; return false; } - QFile myFile(fileName + ".csv"); + QFile myFile(fileName); if (!myFile.open(QIODevice::ReadOnly)) { - myError = "Open CSV failed: " + fileName.toStdString(); - logError(QString::fromStdString(myError)); + errorString = "Failed opening: " + fileName; return false; } @@ -2419,7 +2403,7 @@ bool Project::loadGlocalStationsCsv(QString fileName, std::vector Date: Thu, 14 Nov 2024 16:57:02 +0100 Subject: [PATCH 31/34] fixed daily interpolation time --- interpolation/spatialControl.cpp | 5 +++-- meteo/meteo.cpp | 1 + meteo/quality.cpp | 2 +- pragaProject/pragaProject.cpp | 6 +++--- 4 files changed, 8 insertions(+), 6 deletions(-) diff --git a/interpolation/spatialControl.cpp b/interpolation/spatialControl.cpp index c450e51da..11def4168 100644 --- a/interpolation/spatialControl.cpp +++ b/interpolation/spatialControl.cpp @@ -1,3 +1,4 @@ +#include #include #include @@ -338,7 +339,7 @@ bool checkData(Crit3DQuality* myQuality, meteoVariable myVar, Crit3DMeteoPoint* for (int i = 0; i < nrMeteoPoints; i++) { meteoPoints[i].currentValue = meteoPoints[i].elaboration; - if (int(meteoPoints[i].currentValue) != int(NODATA)) + if (isEqual(meteoPoints[i].currentValue, NODATA)) meteoPoints[i].quality = quality::accepted; else meteoPoints[i].quality = quality::missing_data; @@ -350,7 +351,7 @@ bool checkData(Crit3DQuality* myQuality, meteoVariable myVar, Crit3DMeteoPoint* for (int i = 0; i < nrMeteoPoints; i++) { meteoPoints[i].currentValue = meteoPoints[i].anomaly; - if (int(meteoPoints[i].currentValue) != int(NODATA)) + if (isEqual(meteoPoints[i].currentValue, NODATA)) meteoPoints[i].quality = quality::accepted; else meteoPoints[i].quality = quality::missing_data; diff --git a/meteo/meteo.cpp b/meteo/meteo.cpp index 514bf3749..916464924 100644 --- a/meteo/meteo.cpp +++ b/meteo/meteo.cpp @@ -23,6 +23,7 @@ ftomei@arpae.it */ +#include #include #include diff --git a/meteo/quality.cpp b/meteo/quality.cpp index df282654b..2e5ad6e42 100644 --- a/meteo/quality.cpp +++ b/meteo/quality.cpp @@ -153,7 +153,7 @@ void Crit3DQuality::syntacticQualityControl(meteoVariable myVar, Crit3DMeteoPoin for (int i = 0; i < nrMeteoPoints; i++) { - if (int(meteoPoints[i].currentValue) == int(NODATA)) + if (isEqual(meteoPoints[i].currentValue, NODATA)) meteoPoints[i].quality = quality::missing_data; else { diff --git a/pragaProject/pragaProject.cpp b/pragaProject/pragaProject.cpp index 6f9a120ec..afa1e4542 100644 --- a/pragaProject/pragaProject.cpp +++ b/pragaProject/pragaProject.cpp @@ -2851,7 +2851,7 @@ bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QL QDate loadDateFin = QDate(1800, 1, 1); while (myDate <= dateFin) - { + { countDaysSaving++; // check if load needed @@ -2911,14 +2911,14 @@ bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QL if (getVarFrequency(myVar) == daily) { logInfo(QString::fromStdString(getMeteoVarName(myVar))); - if (! interpolationMeteoGrid(myVar, daily, getCrit3DTime(myDate, myHour))) return false; + if (! interpolationMeteoGrid(myVar, daily, getCrit3DTime(myDate, 1))) return false; } } foreach (myVar, dailyDerivedVars) { logInfo(QString::fromStdString(getMeteoVarName(myVar))); - deriveVariableMeteoGrid(myVar, daily, getCrit3DTime(myDate, myHour)); + deriveVariableMeteoGrid(myVar, daily, getCrit3DTime(myDate, 1)); } } From db99c12bc2d8486d149417aa2cc7d50e35b55497 Mon Sep 17 00:00:00 2001 From: Caterina Toscano Date: Fri, 15 Nov 2024 12:54:38 +0100 Subject: [PATCH 32/34] fixed glocal detrending, added checks upon loading project, added features for multiple detrending in GUI --- interpolation/interpolationSettings.cpp | 10 ++ project/dialogInterpolation.cpp | 227 ++++++++++++++++++++++-- project/dialogInterpolation.h | 7 + project/project.cpp | 110 ++++++++---- project/project.h | 3 +- 5 files changed, 308 insertions(+), 49 deletions(-) diff --git a/interpolation/interpolationSettings.cpp b/interpolation/interpolationSettings.cpp index 8618b983b..5077f9d78 100644 --- a/interpolation/interpolationSettings.cpp +++ b/interpolation/interpolationSettings.cpp @@ -342,8 +342,16 @@ void Crit3DInterpolationSettings::setChosenElevationFunction(TFittingFunction ch int elPos = NODATA; for (int i = 0; i < getProxyNr(); i++) + { if (getProxyPragaName(getProxy(i)->getName()) == proxyHeight) elPos = i; + else //if a proxy has been checked by the user but has no ranges for its parameters, add them + { + if (getProxy(i)->getFittingParametersRange().empty() && getSelectedCombination().isProxyActive(i)) + + getProxy(i)->setFittingParametersRange({-1, -40, 1, 50}); + } + } double MIN_T = -20; double MAX_T = 40; @@ -401,6 +409,8 @@ void Crit3DInterpolationSettings::setChosenElevationFunction(TFittingFunction ch getProxy(elPos)->setFittingFunctionName(chosenFunction); } } + + } void Crit3DInterpolationSettings::setPointsRange(double min, double max) diff --git a/project/dialogInterpolation.cpp b/project/dialogInterpolation.cpp index 8e799db08..2f844b352 100644 --- a/project/dialogInterpolation.cpp +++ b/project/dialogInterpolation.cpp @@ -113,6 +113,12 @@ DialogInterpolation::DialogInterpolation(Project *myProject) optimalDetrendingEdit = new QCheckBox(tr("optimal detrending")); optimalDetrendingEdit->setChecked(_interpolationSettings->getUseBestDetrending()); + if (_interpolationSettings->getUseBestDetrending()) + { + _interpolationSettings->setUseMultipleDetrending(false); + _interpolationSettings->setUseLocalDetrending(false); + _interpolationSettings->setUseGlocalDetrending(false); + } layoutDetrending->addWidget(optimalDetrendingEdit); connect(optimalDetrendingEdit, SIGNAL(stateChanged(int)), this, SLOT(optimalDetrendingChanged(int))); @@ -234,7 +240,6 @@ void DialogInterpolation::upscaleFromDemChanged(int active) void DialogInterpolation::multipleDetrendingChanged(int active) { if (active == Qt::Checked) optimalDetrendingEdit->setChecked(Qt::Unchecked); - optimalDetrendingEdit->setEnabled(active == Qt::Unchecked); } void DialogInterpolation::localDetrendingChanged(int active) @@ -244,21 +249,21 @@ void DialogInterpolation::localDetrendingChanged(int active) optimalDetrendingEdit->setEnabled(active == Qt::Unchecked); if (active == Qt::Checked) glocalDetrendingEdit->setChecked(Qt::Unchecked); glocalDetrendingEdit->setEnabled(active == Qt::Unchecked); + if (active == Qt::Checked) multipleDetrendingEdit->setChecked(Qt::Checked); } void DialogInterpolation::glocalDetrendingChanged(int active) { - if (active == Qt::Checked) topographicDistanceEdit->setChecked(Qt::Unchecked); - topographicDistanceEdit->setEnabled(active == Qt::Unchecked); - maxTdMultiplierEdit.setEnabled(active == Qt::Unchecked); if (active == Qt::Checked) optimalDetrendingEdit->setChecked(Qt::Unchecked); optimalDetrendingEdit->setEnabled(active == Qt::Unchecked); if (active == Qt::Checked) localDetrendingEdit->setChecked(Qt::Unchecked); localDetrendingEdit->setEnabled(active == Qt::Unchecked); + if (active == Qt::Checked) multipleDetrendingEdit->setChecked(Qt::Checked); } void DialogInterpolation::optimalDetrendingChanged(int active) { + if (active == Qt::Checked) multipleDetrendingEdit->setChecked(Qt::Unchecked); if (active == Qt::Checked) localDetrendingEdit->setChecked(Qt::Unchecked); localDetrendingEdit->setEnabled(active == Qt::Unchecked); if (active == Qt::Checked) glocalDetrendingEdit->setChecked(Qt::Unchecked); @@ -387,6 +392,61 @@ void ProxyDialog::changedProxy(bool savePrevious) _field.setCurrentIndex(index); _proxyGridName.setText(QString::fromStdString(myProxy->getGridName())); _forQuality.setChecked(myProxy->getForQualityControl()); + + _param0.clear(); + _param1.clear(); + _param2.clear(); + _param3.clear(); + _param4.clear(); + _param5.clear(); + + if (_project->interpolationSettings.getUseMultipleDetrending()) + { + if (!myProxy->getFittingParametersMin().empty() && !myProxy->getFittingParametersMax().empty()) + { + _param0.setText(QString("%1").arg(myProxy->getFittingParametersMin()[0], 0, 'f', 4) + ", " + QString("%1").arg(myProxy->getFittingParametersMax()[0], 0, 'f', 4)); + _param1.setText(QString("%1").arg(myProxy->getFittingParametersMin()[1], 0, 'f', 4) + ", " + QString("%1").arg(myProxy->getFittingParametersMax()[1], 0, 'f', 4)); + if (myProxy->getFittingParametersMin().size() > 2 && myProxy->getFittingParametersMax().size() > 2) { + _param2.setText(QString("%1").arg(myProxy->getFittingParametersMin()[2], 0, 'f', 4) + ", " + QString("%1").arg(myProxy->getFittingParametersMax()[2], 0, 'f', 4)); + _param3.setText(QString("%1").arg(myProxy->getFittingParametersMin()[3], 0, 'f', 4) + ", " + QString("%1").arg(myProxy->getFittingParametersMax()[3], 0, 'f', 4)); + } + if (myProxy->getFittingParametersMin().size() > 4 && myProxy->getFittingParametersMax().size() > 4) + _param4.setText(QString("%1").arg(myProxy->getFittingParametersMin()[4], 0, 'f', 4) + ", " + QString("%1").arg(myProxy->getFittingParametersMax()[4], 0, 'f', 4)); + if (myProxy->getFittingParametersMin().size() > 5 && myProxy->getFittingParametersMax().size() > 5) + _param5.setText(QString("%1").arg(myProxy->getFittingParametersMin()[5], 0, 'f', 4) + ", " + QString("%1").arg(myProxy->getFittingParametersMax()[5], 0, 'f', 4)); + + _param0.setEnabled(true); + _param1.setEnabled(true); + _param2.setEnabled(true); + _param3.setEnabled(true); + _param4.setEnabled(true); + _param5.setEnabled(true); + } + if (getProxyPragaName(myProxy->getName()) != proxyHeight) + { + _param2.setEnabled(false); + _param3.setEnabled(false); + _param4.setEnabled(false); + _param5.setEnabled(false); + } + else + { + if (myProxy->getFittingFunctionName() != piecewiseThreeFree) + _param5.setEnabled(false); + if (myProxy->getFittingFunctionName() != piecewiseThree) + _param4.setEnabled(false); + } + } + else + { + _param0.setEnabled(false); + _param1.setEnabled(false); + _param2.setEnabled(false); + _param3.setEnabled(false); + _param4.setEnabled(false); + _param5.setEnabled(false); + } + } void ProxyDialog::selectGridFile() @@ -463,16 +523,67 @@ void ProxyDialog::saveProxy(Crit3DProxy* myProxy) if (_proxyGridName.text() != "") myProxy->setGridName(_proxyGridName.text().toStdString()); + std::vector parametersMin; + std::vector parametersMax; + QString temp; + QStringList tempList; + + temp.clear(); + tempList.clear(); + parametersMin.clear(); + parametersMax.clear(); + + temp = _param0.toPlainText(); + if (temp != "") + { + tempList.append(temp.split(QRegularExpression(","))); + temp = _param1.toPlainText(); + tempList.append(temp.split(QRegularExpression(","))); + temp = _param2.toPlainText(); + if (temp != "") + { + tempList.append(temp.split(QRegularExpression(","))); + temp = _param3.toPlainText(); + tempList.append(temp.split(QRegularExpression(","))); + temp = _param4.toPlainText(); + if (temp != "") + { + tempList.append(temp.split(QRegularExpression(","))); + temp = _param5.toPlainText(); + if (temp != "") + tempList.append(temp.split(QRegularExpression(","))); + } + } + } + + for (int j = 0; j < tempList.size(); j = j + 2) + { + parametersMin.push_back(tempList[j].toDouble()); + parametersMax.push_back(tempList[j+1].toDouble()); + } + + std::vector parameters = parametersMin; + for (int j = 0; j < parametersMax.size(); j++) + parameters.push_back(parametersMax[j]); + + myProxy->setFittingParametersRange(parameters); + myProxy->setForQualityControl(_forQuality.isChecked()); } ProxyDialog::ProxyDialog(Project *myProject) { - QVBoxLayout *layoutMain = new QVBoxLayout(); + QVBoxLayout *layoutLeft = new QVBoxLayout(); QHBoxLayout *layoutProxyCombo = new QHBoxLayout(); QVBoxLayout *layoutProxy = new QVBoxLayout(); QVBoxLayout *layoutPointValues = new QVBoxLayout(); QVBoxLayout *layoutGrid = new QVBoxLayout(); + QGroupBox *groupParameters = new QGroupBox(tr("Multiple detrending parameters")); + QHBoxLayout *layoutParametersUp = new QHBoxLayout; + QHBoxLayout *layoutParametersMiddle = new QHBoxLayout; + QHBoxLayout *layoutParametersDown = new QHBoxLayout; + QVBoxLayout *layoutRight = new QVBoxLayout; + QHBoxLayout *layoutMain = new QHBoxLayout(); this->resize(300, 100); @@ -499,7 +610,7 @@ ProxyDialog::ProxyDialog(Project *myProject) connect(_delete, &QPushButton::clicked, [=](){ this->deleteProxy(); }); layoutProxy->addLayout(layoutProxyCombo); - layoutMain->addLayout(layoutProxy); + layoutLeft->addLayout(layoutProxy); QLabel *labelTableList = new QLabel(tr("table for point values")); layoutPointValues->addWidget(labelTableList); @@ -514,7 +625,7 @@ ProxyDialog::ProxyDialog(Project *myProject) layoutPointValues->addWidget(labelFieldList); layoutPointValues->addWidget(&_field); - layoutMain->addLayout(layoutPointValues); + layoutLeft->addLayout(layoutPointValues); // grid QLabel *labelGrid = new QLabel(tr("proxy grid")); @@ -524,11 +635,104 @@ ProxyDialog::ProxyDialog(Project *myProject) _proxyGridName.setEnabled(false); layoutGrid->addWidget(&_proxyGridName); connect(_selectGrid, &QPushButton::clicked, [=](){ this->selectGridFile(); }); - layoutMain->addLayout(layoutGrid); + layoutLeft->addLayout(layoutGrid); // quality control _forQuality.setText("use for quality control"); - layoutMain->addWidget(&_forQuality); + layoutLeft->addWidget(&_forQuality); + + proxyIndex = NODATA; + + if (_proxyCombo.count() > 0) + proxyIndex = 0; + + //parameters + + const double H0_MIN = -350; //height of single inversion point (double piecewise) or first inversion point (triple piecewise) + const double H0_MAX = 2500; + 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; + + QLabel *par0Label = new QLabel(tr("par0\n(min, max)")); + QLabel *par1Label = new QLabel(tr("par1\n(min, max)")); + QLabel *par2Label = new QLabel(tr("par2\n(min, max)")); + QLabel *par3Label = new QLabel(tr("par3\n(min, max)")); + QLabel *par4Label = new QLabel(tr("par4\n(min, max)")); + QLabel *par5Label = new QLabel(tr("par5\n(min, max)")); + + std::vector parametersMin = _proxy[proxyIndex].getFittingParametersMin(); + std::vector parametersMax = _proxy[proxyIndex].getFittingParametersMax(); + + if ((parametersMin.empty() || parametersMax.empty()) && _project->interpolationSettings.getSelectedCombination().isProxyActive(proxyIndex)) + { + //se il proxy è stato attivato e non ci sono parametri caricati in precedenza, carica i default + if (_proxy[proxyIndex].getFittingFunctionName() == piecewiseTwo) + { + _proxy[proxyIndex].setFittingParametersRange({0, -20, SLOPE_MIN, INVSLOPE_MIN, + H0_MAX, 40, SLOPE_MAX, INVSLOPE_MAX}); + _proxy[proxyIndex].setFittingFirstGuess({0,1,1,1}); + } + else if (_proxy[proxyIndex].getFittingFunctionName() == piecewiseThree) + { + _proxy[proxyIndex].setFittingParametersRange({H0_MIN, -20, DELTA_MIN, SLOPE_MIN, INVSLOPE_MIN, + H0_MAX, 40, DELTA_MAX, SLOPE_MAX, INVSLOPE_MAX}); + _proxy[proxyIndex].setFittingFirstGuess({0,1,1,1,1}); + } + else if (_proxy[proxyIndex].getFittingFunctionName() == piecewiseThreeFree) + { + _proxy[proxyIndex].setFittingParametersRange({H0_MIN, -20, DELTA_MIN, SLOPE_MIN, INVSLOPE_MIN, INVSLOPE_MIN, + H0_MAX, 40, DELTA_MAX, SLOPE_MAX, INVSLOPE_MAX, INVSLOPE_MAX}); + _proxy[proxyIndex].setFittingFirstGuess({0,1,1,1,1,1}); + } + else + _proxy[proxyIndex].setFittingParametersRange({-1, 50, 1, -40}); + + parametersMin = _proxy[proxyIndex].getFittingParametersMin(); + parametersMax = _proxy[proxyIndex].getFittingParametersMax(); + } + + _param0.setText(QString("%1").arg(parametersMin[0], 0, 'f', 4) + ", " + QString("%1").arg(parametersMax[0], 0, 'f', 4)); + _param1.setText(QString("%1").arg(parametersMin[1], 0, 'f', 4) + ", " + QString("%1").arg(parametersMax[1], 0, 'f', 4)); + if (parametersMin.size() > 2 && parametersMax.size() > 2) { + _param2.setText(QString("%1").arg(parametersMin[2], 0, 'f', 4) + ", " + QString("%1").arg(parametersMax[2], 0, 'f', 4)); + _param3.setText(QString("%1").arg(parametersMin[3], 0, 'f', 4) + ", " + QString("%1").arg(parametersMax[3], 0, 'f', 4)); + } + if (parametersMin.size() > 4 && parametersMax.size() > 4) + _param4.setText(QString("%1").arg(parametersMin[4], 0, 'f', 4) + ", " + QString("%1").arg(parametersMax[4], 0, 'f', 4)); + if (parametersMin.size() > 5 && parametersMax.size() > 5) + _param5.setText(QString("%1").arg(parametersMin[5], 0, 'f', 4) + ", " + QString("%1").arg(parametersMax[5], 0, 'f', 4)); + + layoutParametersUp->addWidget(par0Label); + _param0.setFixedHeight(45); + layoutParametersUp->addWidget(&_param0); + layoutParametersUp->addWidget(par3Label); + _param3.setFixedHeight(45); + layoutParametersUp->addWidget(&_param3); + layoutParametersMiddle->addWidget(par1Label); + _param1.setFixedHeight(45); + layoutParametersMiddle->addWidget(&_param1); + layoutParametersMiddle->addWidget(par4Label); + _param4.setFixedHeight(45); + layoutParametersMiddle->addWidget(&_param4); + layoutParametersDown->addWidget(par2Label); + _param2.setFixedHeight(45); + layoutParametersDown->addWidget(&_param2); + layoutParametersDown->addWidget(par5Label); + _param5.setFixedHeight(45); + layoutParametersDown->addWidget(&_param5); + + layoutRight->addLayout(layoutParametersUp); + layoutRight->addLayout(layoutParametersMiddle); + layoutRight->addLayout(layoutParametersDown); + + groupParameters->setLayout(layoutRight); + + layoutMain->addLayout(layoutLeft); + layoutMain->addWidget(groupParameters); // buttons QDialogButtonBox *buttonBox = new QDialogButtonBox(QDialogButtonBox::Ok | QDialogButtonBox::Cancel); @@ -579,6 +783,7 @@ bool ProxyDialog::checkProxies(QString *error) void ProxyDialog::saveProxies() { + Crit3DProxyCombination myCombination = _project->interpolationSettings.getSelectedCombination(); _project->interpolationSettings.initializeProxy(); _project->qualityInterpolationSettings.initializeProxy(); @@ -587,9 +792,9 @@ void ProxyDialog::saveProxies() std::vector proxyOrder; for (unsigned i=0; i < _proxy.size(); i++) - { + { proxyList.push_back(_proxy[i]); - proxyActive.push_back(true); + proxyActive.push_back(myCombination.isProxyActive(i)); proxyOrder.push_back(i+1); } diff --git a/project/dialogInterpolation.h b/project/dialogInterpolation.h index 4a26141da..52c6c0a8d 100644 --- a/project/dialogInterpolation.h +++ b/project/dialogInterpolation.h @@ -71,6 +71,13 @@ QComboBox _table; QLineEdit _proxyGridName; QCheckBox _forQuality; + QTextEdit _param0; + QTextEdit _param1; + QTextEdit _param2; + QTextEdit _param3; + QTextEdit _param4; + QTextEdit _param5; + void changedProxy(bool savePrevious); void changedTable(); diff --git a/project/project.cpp b/project/project.cpp index 8343eedba..6a62c36d0 100644 --- a/project/project.cpp +++ b/project/project.cpp @@ -219,34 +219,63 @@ bool Project::checkProxy(Crit3DProxy &myProxy, QString* error, bool isActive) } if (interpolationSettings.getUseMultipleDetrending() && isActive) + checkProxyForMultipleDetrending(myProxy, isHeight); + + return true; +} + +void Project::checkProxyForMultipleDetrending(Crit3DProxy &myProxy, bool isHeight) +{ + const double H0_MIN = -350; //height of single inversion point (double piecewise) or first inversion point (triple piecewise) + const double H0_MAX = 2500; + 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; + + int nrParameters = 0; + if (isHeight) { - int nrParameters = 0; - if (isHeight) - { - if (myProxy.getFittingFunctionName() == piecewiseTwo) - nrParameters = 4; - else if (myProxy.getFittingFunctionName() == piecewiseThree) - nrParameters = 5; - else if (myProxy.getFittingFunctionName() == piecewiseThreeFree) - nrParameters = 6; - } - else - nrParameters = 2; + if (myProxy.getFittingFunctionName() == piecewiseTwo) + nrParameters = 4; + else if (myProxy.getFittingFunctionName() == piecewiseThree) + nrParameters = 5; + else if (myProxy.getFittingFunctionName() == piecewiseThreeFree) + nrParameters = 6; + } + else + nrParameters = 2; - if (myProxy.getFittingParametersRange().size() != nrParameters*2) - { - *error = "wrong numer of parameters for proxy: " + QString::fromStdString(name_); - return false; - } + if (myProxy.getFittingParametersRange().size() != nrParameters*2) + { + logError("wrong number of parameters for proxy: " + QString::fromStdString(myProxy.getName()) + ". You can change them in the parameters.ini file, but default parameters will be used for now."); + if (nrParameters == 2) + myProxy.setFittingParametersRange({-10, -50, 11, 50}); + else if (nrParameters == 4) + myProxy.setFittingParametersRange({0, -30, SLOPE_MIN, INVSLOPE_MIN, + H0_MAX, 50, SLOPE_MAX, INVSLOPE_MAX}); + else if (nrParameters == 5) + myProxy.setFittingParametersRange({H0_MIN, -30, DELTA_MIN, SLOPE_MIN, INVSLOPE_MIN, + H0_MAX, 50, DELTA_MAX, SLOPE_MAX, INVSLOPE_MAX}); + else if (nrParameters == 6) + myProxy.setFittingParametersRange({H0_MIN, -30, DELTA_MIN, SLOPE_MIN, INVSLOPE_MIN, INVSLOPE_MIN, + H0_MAX, 50, DELTA_MAX, SLOPE_MAX, INVSLOPE_MAX, INVSLOPE_MAX}); + } - if (isHeight && myProxy.getFittingFirstGuess().size() != nrParameters) - { - *error = "wrong number of first guess settings for proxy: " + QString::fromStdString(name_); - return false; - } + if (isHeight && myProxy.getFittingFirstGuess().size() != nrParameters) + { + logError("wrong number of first guess settings for proxy: " + QString::fromStdString(myProxy.getName()) + ". You can change them in the parameters.ini file, but default settings will be used for now."); + if (nrParameters == 4) + myProxy.setFittingFirstGuess({0, 1, 1, 1}); + else if (nrParameters == 5) + myProxy.setFittingFirstGuess({0, 1, 1, 1, 1}); + else if (nrParameters == 6) + myProxy.setFittingFirstGuess({0, 1, 1, 1, 1, 1}); } - return true; + return; } @@ -626,30 +655,37 @@ bool Project::loadParameters(QString parametersFileName) if (parametersSettings->contains("topographicDistance")) interpolationSettings.setUseTD(parametersSettings->value("topographicDistance").toBool()); + if (parametersSettings->contains("optimalDetrending")) + interpolationSettings.setUseBestDetrending(parametersSettings->value("optimalDetrending").toBool()); + + if (parametersSettings->contains("multipleDetrending")) + interpolationSettings.setUseMultipleDetrending(parametersSettings->value("multipleDetrending").toBool()); + if (parametersSettings->contains("localDetrending")) interpolationSettings.setUseLocalDetrending(parametersSettings->value("localDetrending").toBool()); if (parametersSettings->contains("glocalDetrending")) interpolationSettings.setUseGlocalDetrending(parametersSettings->value("glocalDetrending").toBool()); + if (interpolationSettings.getUseMultipleDetrending() && interpolationSettings.getUseBestDetrending()) + { + interpolationSettings.setUseMultipleDetrending(false); + interpolationSettings.setUseLocalDetrending(false); + interpolationSettings.setUseGlocalDetrending(false); + } + if (parametersSettings->contains("glocalMapName")) this->glocalMapName = parametersSettings->value("glocalMapName").toString(); if (parametersSettings->contains("meteogrid_upscalefromdem")) interpolationSettings.setMeteoGridUpscaleFromDem(parametersSettings->value("meteogrid_upscalefromdem").toBool()); - if (parametersSettings->contains("multipleDetrending")) - interpolationSettings.setUseMultipleDetrending(parametersSettings->value("multipleDetrending").toBool()); - if (parametersSettings->contains("lapseRateCode")) { interpolationSettings.setUseLapseRateCode(parametersSettings->value("lapseRateCode").toBool()); qualityInterpolationSettings.setUseLapseRateCode(parametersSettings->value("lapseRateCode").toBool()); } - if (parametersSettings->contains("optimalDetrending")) - interpolationSettings.setUseBestDetrending(parametersSettings->value("optimalDetrending").toBool()); - if (parametersSettings->contains("minRegressionR2")) { interpolationSettings.setMinRegressionR2(parametersSettings->value("minRegressionR2").toFloat()); @@ -2277,7 +2313,7 @@ bool Project::loadGlocalStationsAndCells(bool isGrid) //leggi csv aree QString mapsFolder = defaultPath + PATH_GEO; std::string fileNameStations = mapsFolder.toStdString() + glocalMapName.toStdString() + "_stations"; - std::vector> areaPoints; + std::vector> areaPoints; std::string myError; loadGlocalStationsCsv(QString::fromStdString(fileNameStations), areaPoints); @@ -2294,7 +2330,7 @@ bool Project::loadGlocalStationsAndCells(bool isGrid) { //controlla se id si trova nel vettore areaPoints[j] e salva l'indice del meteopoint for (int i=0; i < nrMeteoPoints; i++) - if (areaPoints[j][k] == std::stoi(meteoPoints[i].id)) + if (areaPoints[j][k] == meteoPoints[i].id) { temp.push_back(i); break; @@ -2338,7 +2374,7 @@ bool Project::loadGlocalWeightMaps(std::vector &myAreas, bool i if (!isGrid) { nrCols = DEM.header->nrCols; - nrRows = DEM.header->nrCols; + nrRows = DEM.header->nrRows; } else { @@ -2387,7 +2423,7 @@ bool Project::loadGlocalWeightMaps(std::vector &myAreas, bool i return true; } -bool Project::loadGlocalStationsCsv(QString fileName, std::vector> &areaPoints) +bool Project::loadGlocalStationsCsv(QString fileName, std::vector> &areaPoints) { std::string myError; @@ -2425,13 +2461,13 @@ bool Project::loadGlocalStationsCsv(QString fileName, std::vector temp; + std::vector temp; while(!myStream.atEnd()) { nrLine++; line = myStream.readLine().split(','); for (int i = 1; i < line.size(); i++) - temp.push_back(line[i].toInt()); + temp.push_back(line[i].toStdString()); areaPoints.resize(line[0].toInt() + 1); areaPoints[line[0].toInt()] = temp; @@ -3750,8 +3786,8 @@ void Project::saveProxies() if (myProxy->getProxyField() != "") parametersSettings->setValue("field", QString::fromStdString(myProxy->getProxyField())); if (myProxy->getGridName() != "") parametersSettings->setValue("raster", getRelativePath(QString::fromStdString(myProxy->getGridName()))); if (myProxy->getStdDevThreshold() != NODATA) parametersSettings->setValue("stddev_threshold", QString::number(myProxy->getStdDevThreshold())); - if (myProxy->getFittingParametersRange().size() > 0) parametersSettings->setValue("fitting_parameters_min", DoubleVectorToStringList(myProxy->getFittingParametersMin())); - if (myProxy->getFittingParametersRange().size() > 0) parametersSettings->setValue("fitting_parameters_max", DoubleVectorToStringList(myProxy->getFittingParametersMax())); + if (!myProxy->getFittingParametersMin().empty()) parametersSettings->setValue("fitting_parameters_min", DoubleVectorToStringList(myProxy->getFittingParametersMin())); + if (!myProxy->getFittingParametersMax().empty()) parametersSettings->setValue("fitting_parameters_max", DoubleVectorToStringList(myProxy->getFittingParametersMax())); if (myProxy->getFittingFirstGuess().size() > 0) parametersSettings->setValue("fitting_first_guess", IntVectorToStringList(myProxy->getFittingFirstGuess())); if (getProxyPragaName(myProxy->getName()) == proxyHeight && myProxy->getFittingFunctionName() != noFunction) parametersSettings->setValue("fitting_function", QString::fromStdString(getKeyStringElevationFunction(myProxy->getFittingFunctionName()))); parametersSettings->endGroup(); diff --git a/project/project.h b/project/project.h index 75cf4c809..5e5380bf3 100644 --- a/project/project.h +++ b/project/project.h @@ -187,6 +187,7 @@ bool checkProxy(Crit3DProxy &myProxy, QString *error, bool isActive); bool addProxyToProject(std::vector proxyList, std::deque proxyActive, std::vector proxyOrder); void addProxyGridSeries(QString name_, std::vector gridNames, std::vector gridYears); + void checkProxyForMultipleDetrending(Crit3DProxy &myProxy, bool isHeight); void setCurrentDate(QDate myDate); void setCurrentHour(int myHour); void setCurrentVariable(meteoVariable variable); @@ -264,7 +265,7 @@ bool loadGlocalAreasMap(); bool loadGlocalStationsAndCells(bool isGrid); bool loadGlocalWeightMaps(std::vector &myAreas, bool isGrid); - bool loadGlocalStationsCsv(QString fileName, std::vector> &areaPoints); + bool loadGlocalStationsCsv(QString fileName, std::vector > &areaPoints); bool groupCellsInArea(std::vector &areaPoints, unsigned index, bool isGrid); bool glocalWeightsMaps(float windowWidth); void macroAreaDetrending(Crit3DMacroArea myArea, meteoVariable myVar, std::vector interpolationPoints, std::vector &subsetInterpolationPoints, int elevationPos); From 7a3798ea1e46c3c473dbb984fa204b945ee8a396 Mon Sep 17 00:00:00 2001 From: Caterina Toscano Date: Fri, 15 Nov 2024 13:04:00 +0100 Subject: [PATCH 33/34] fix conflict --- project/project.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/project/project.cpp b/project/project.cpp index 2e078175b..251eaa6d3 100644 --- a/project/project.cpp +++ b/project/project.cpp @@ -2303,7 +2303,7 @@ bool Project::loadGlocalStationsAndCells(bool isGrid) { //leggi csv aree QString fileNameStations = getCompleteFileName(glocalPointsName, PATH_GEO); - std::vector> areaPoints; + std::vector> areaPoints; if (! loadGlocalStationsCsv(fileNameStations, areaPoints)) return false; From 3bf95c7f71988d65b3e3e2921ee777f9710a1ede Mon Sep 17 00:00:00 2001 From: Caterina Toscano Date: Mon, 18 Nov 2024 11:03:03 +0100 Subject: [PATCH 34/34] fix --- project/project.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/project/project.cpp b/project/project.cpp index 251eaa6d3..f61431ab8 100644 --- a/project/project.cpp +++ b/project/project.cpp @@ -3127,7 +3127,7 @@ bool Project::interpolationDemMain(meteoVariable myVar, const Crit3DTime& myTime if (interpolationSettings.getUseGlocalDetrending() && ! interpolationSettings.isGlocalReady()) { if (! loadGlocalAreasMap()) return false; - if (! loadGlocalStationsAndCells(!interpolationSettings.getMeteoGridUpscaleFromDem())) return false; + if (! loadGlocalStationsAndCells(false)) return false; } // solar radiation model