From cbab5e06e1c98ac9baa875de52ed2645a770b31a Mon Sep 17 00:00:00 2001 From: Caterina Toscano Date: Mon, 28 Oct 2024 13:12:13 +0100 Subject: [PATCH] glocal detrending with buffer areas --- agrolib/interpolation/interpolation.cpp | 16 +- agrolib/interpolation/interpolation.h | 2 +- .../interpolation/interpolationSettings.cpp | 4 +- agrolib/interpolation/interpolationSettings.h | 6 +- agrolib/pragaProject/pragaProject.cpp | 7 - agrolib/project/project.cpp | 156 +++++++++++++++--- agrolib/project/project.h | 4 +- agrolib/proxyWidget/localProxyWidget.cpp | 11 +- agrolib/proxyWidget/proxyWidget.cpp | 2 +- src/mainWindow.cpp | 2 +- 10 files changed, 163 insertions(+), 47 deletions(-) diff --git a/agrolib/interpolation/interpolation.cpp b/agrolib/interpolation/interpolation.cpp index a3139d781..0cd8d5e86 100644 --- a/agrolib/interpolation/interpolation.cpp +++ b/agrolib/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/agrolib/interpolation/interpolation.h b/agrolib/interpolation/interpolation.h index 63b6c6a22..99fdfcf57 100644 --- a/agrolib/interpolation/interpolation.h +++ b/agrolib/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/agrolib/interpolation/interpolationSettings.cpp b/agrolib/interpolation/interpolationSettings.cpp index e6c5860e5..f286df8d0 100644 --- a/agrolib/interpolation/interpolationSettings.cpp +++ b/agrolib/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/agrolib/interpolation/interpolationSettings.h b/agrolib/interpolation/interpolationSettings.h index cc665edd8..4371d82bb 100644 --- a/agrolib/interpolation/interpolationSettings.h +++ b/agrolib/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/agrolib/pragaProject/pragaProject.cpp b/agrolib/pragaProject/pragaProject.cpp index 0ba7aa533..27b4aeecc 100644 --- a/agrolib/pragaProject/pragaProject.cpp +++ b/agrolib/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/agrolib/project/project.cpp b/agrolib/project/project.cpp index c13a0bd65..2e69f48dc 100644 --- a/agrolib/project/project.cpp +++ b/agrolib/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/agrolib/project/project.h b/agrolib/project/project.h index c39ef7a4a..46302b3d1 100644 --- a/agrolib/project/project.h +++ b/agrolib/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/agrolib/proxyWidget/localProxyWidget.cpp b/agrolib/proxyWidget/localProxyWidget.cpp index cb9600bfd..93ef34072 100644 --- a/agrolib/proxyWidget/localProxyWidget.cpp +++ b/agrolib/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/agrolib/proxyWidget/proxyWidget.cpp b/agrolib/proxyWidget/proxyWidget.cpp index 3409034b4..4140ec1c3 100644 --- a/agrolib/proxyWidget/proxyWidget.cpp +++ b/agrolib/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(); diff --git a/src/mainWindow.cpp b/src/mainWindow.cpp index d0b4d6941..8893dac8d 100644 --- a/src/mainWindow.cpp +++ b/src/mainWindow.cpp @@ -6775,7 +6775,7 @@ void MainWindow::on_actionInterpolationTopographicIndex_triggered() void MainWindow::on_actionInterpolationWriteGlocalWeightMaps_triggered() { - if (myProject.loadGlocalFiles() < 1) + if (!myProject.loadGlocalAreasMap()) { myProject.logError("Error loading zone grid " + myProject.glocalMapName); return;