From 58756b0adad68be602819d071b461667b14ab686 Mon Sep 17 00:00:00 2001 From: ftomei Date: Wed, 29 Nov 2023 16:32:22 +0100 Subject: [PATCH 1/6] improve db dates --- project/project.cpp | 4 ++++ project/project.h | 1 + 2 files changed, 5 insertions(+) diff --git a/project/project.cpp b/project/project.cpp index 81d0edce6..94b62b815 100644 --- a/project/project.cpp +++ b/project/project.cpp @@ -63,6 +63,8 @@ void Project::initializeProject() outputPointsDbHandler = nullptr; meteoGridDbHandler = nullptr; aggregationDbHandler = nullptr; + meteoPointsDbFirstTime.setTime_t(0); + meteoPointsDbLastTime.setTime_t(0); meteoSettings->initialize(); quality->initialize(); @@ -1091,6 +1093,7 @@ bool Project::loadMeteoPointsDB(QString fileName) // find dates meteoPointsDbLastTime = findDbPointLastTime(); + meteoPointsDbFirstTime.setTime_t(0); if (! meteoPointsDbLastTime.isNull()) { @@ -1187,6 +1190,7 @@ bool Project::loadAggregationDBAsMeteoPoints(QString fileName) // find dates meteoPointsDbLastTime = findDbPointLastTime(); + meteoPointsDbFirstTime.setTime_t(0); if (! meteoPointsDbLastTime.isNull()) { diff --git a/project/project.h b/project/project.h index 9c4e94701..5f312fe16 100644 --- a/project/project.h +++ b/project/project.h @@ -110,6 +110,7 @@ Crit3DMeteoPointsDbHandler* meteoPointsDbHandler; Crit3DOutputPointsDbHandler* outputPointsDbHandler; Crit3DAggregationsDbHandler* aggregationDbHandler; + QDateTime meteoPointsDbFirstTime; QDateTime meteoPointsDbLastTime; Crit3DColorScale* meteoPointsColorScale; From f38a055af85db96936a85a7f14e238ac0a4923cf Mon Sep 17 00:00:00 2001 From: ftomei Date: Wed, 29 Nov 2023 17:18:20 +0100 Subject: [PATCH 2/6] update dbmeteohandler --- dbMeteoPoints/dbMeteoPointsHandler.cpp | 50 ++++++++++++++------------ dbMeteoPoints/dbMeteoPointsHandler.h | 6 ++-- pragaProject/pragaProject.cpp | 33 ++++++++++++++--- pragaProject/pragaProject.h | 2 +- 4 files changed, 61 insertions(+), 30 deletions(-) diff --git a/dbMeteoPoints/dbMeteoPointsHandler.cpp b/dbMeteoPoints/dbMeteoPointsHandler.cpp index da8753834..ca9414369 100644 --- a/dbMeteoPoints/dbMeteoPointsHandler.cpp +++ b/dbMeteoPoints/dbMeteoPointsHandler.cpp @@ -502,12 +502,12 @@ bool Crit3DMeteoPointsDbHandler::deleteAllPointsFromDataset(QList datas } -bool Crit3DMeteoPointsDbHandler::loadDailyData(Crit3DDate firstDate, Crit3DDate lastDate, Crit3DMeteoPoint *meteoPoint) +bool Crit3DMeteoPointsDbHandler::loadDailyData(const Crit3DDate &firstDate, const Crit3DDate &lastDate, Crit3DMeteoPoint *meteoPoint) { // check dates if (firstDate > lastDate) { - this->errorStr = "wrong dates: first > last"; + errorStr = "wrong dates: first > last"; return false; } @@ -545,57 +545,61 @@ bool Crit3DMeteoPointsDbHandler::loadDailyData(Crit3DDate firstDate, Crit3DDate } -bool Crit3DMeteoPointsDbHandler::loadHourlyData(Crit3DDate dateStart, Crit3DDate dateEnd, Crit3DMeteoPoint *meteoPoint) +bool Crit3DMeteoPointsDbHandler::loadHourlyData(const Crit3DDate &firstDate, const Crit3DDate &lastDate, Crit3DMeteoPoint *meteoPoint) { - meteoVariable variable; - int idVar; - float value; + // check dates + if (firstDate > lastDate) + { + errorStr = "wrong dates: first > last"; + return false; + } - int numberOfDays = difference(dateStart, dateEnd)+1; + // initialize obs data + int numberOfDays = difference(firstDate, lastDate) + 1; int myHourlyFraction = 1; - QString startDate = QString::fromStdString(dateStart.toStdString()); - QString endDate = QString::fromStdString(dateEnd.toStdString()); - - QSqlQuery qry(_db); - - meteoPoint->initializeObsDataH(myHourlyFraction, numberOfDays, dateStart); + meteoPoint->initializeObsDataH(myHourlyFraction, numberOfDays, firstDate); + QString startDateStr = QString::fromStdString(firstDate.toStdString()); + QString endDateStr = QString::fromStdString(lastDate.toStdString()); QString tableName = QString::fromStdString(meteoPoint->id) + "_H"; QString statement = QString( "SELECT * FROM `%1` WHERE date_time >= DATETIME('%2 01:00:00') AND date_time <= DATETIME('%3 00:00:00', '+1 day')") - .arg(tableName, startDate, endDate); - if( !qry.exec(statement) ) + .arg(tableName, startDateStr, endDateStr); + QSqlQuery qry(_db); + if(! qry.exec(statement) ) { - qDebug() << qry.lastError(); + errorStr = qry.lastError().text(); return false; } else { + meteoVariable variable; while (qry.next()) { QDateTime d = qry.value(0).toDateTime(); Crit3DDate myDate = Crit3DDate(d.date().day(), d.date().month(), d.date().year()); - //myDate = QDate::fromString(dateStr.mid(0,10), "yyyy-MM-dd"); - //myTime = QTime::fromString(dateStr.mid(11,8), "HH:mm:ss"); - //QDateTime d(QDateTime(myDate, myTime, Qt::UTC)); - idVar = qry.value(1).toInt(); - try { + int idVar = qry.value(1).toInt(); + try + { variable = _mapIdMeteoVar.at(idVar); } - catch (const std::out_of_range& ) { + catch (const std::out_of_range& ) + { variable = noMeteoVar; } if (variable != noMeteoVar) { - value = qry.value(2).toFloat(); + float value = qry.value(2).toFloat(); meteoPoint->setMeteoPointValueH(myDate, d.time().hour(), d.time().minute(), variable, value); // copy scalar intensity to vector intensity (instantaneous values are equivalent, following WMO) // should be removed when hourly averages are available if (variable == windScalarIntensity) + { meteoPoint->setMeteoPointValueH(myDate, d.time().hour(), d.time().minute(), windVectorIntensity, value); + } } } } diff --git a/dbMeteoPoints/dbMeteoPointsHandler.h b/dbMeteoPoints/dbMeteoPointsHandler.h index 08940e901..d6d78c1fd 100644 --- a/dbMeteoPoints/dbMeteoPointsHandler.h +++ b/dbMeteoPoints/dbMeteoPointsHandler.h @@ -64,11 +64,13 @@ const gis::Crit3DGisSettings& gisSettings, QString& errorString); bool getPropertiesGivenId(QString id, Crit3DMeteoPoint* meteoPoint, const gis::Crit3DGisSettings& gisSettings, QString& errorString); - bool loadDailyData(Crit3DDate firstDate, Crit3DDate lastDate, Crit3DMeteoPoint *meteoPoint); + bool loadDailyData(const Crit3DDate &firstDate, const Crit3DDate &lastDate, Crit3DMeteoPoint *meteoPoint); std::vector loadDailyVar(QString *myError, meteoVariable variable, Crit3DDate dateStart, Crit3DDate dateEnd, QDate* firstDateDB, Crit3DMeteoPoint *meteoPoint); - bool loadHourlyData(Crit3DDate dateStart, Crit3DDate dateEnd, Crit3DMeteoPoint *meteoPoint); + + bool loadHourlyData(const Crit3DDate &firstDate, const Crit3DDate &lastDate, Crit3DMeteoPoint *meteoPoint); + std::vector loadHourlyVar(QString *myError, meteoVariable variable, Crit3DDate dateStart, Crit3DDate dateEnd, QDateTime* firstDateDB, Crit3DMeteoPoint *meteoPoint); diff --git a/pragaProject/pragaProject.cpp b/pragaProject/pragaProject.cpp index 0b43c0921..83145397a 100644 --- a/pragaProject/pragaProject.cpp +++ b/pragaProject/pragaProject.cpp @@ -1887,7 +1887,7 @@ bool PragaProject::hourlyDerivedVariablesGrid(QDate first, QDate last, bool load } -bool PragaProject::interpolationOutputPointsPeriod(QDate dateIni, QDate dateFin, QList variables) +bool PragaProject::interpolationOutputPointsPeriod(QDate firstDate, QDate lastDate, QList variables) { // check if (variables.size() == 0) @@ -1909,9 +1909,9 @@ bool PragaProject::interpolationOutputPointsPeriod(QDate dateIni, QDate dateFin, } // check dates - if (dateIni.isNull() || dateFin.isNull() || dateIni > dateFin) + if (firstDate.isNull() || lastDate.isNull() || firstDate > lastDate) { - errorString = "Wrong period."; + errorString = "Wrong dates"; return false; } @@ -1944,7 +1944,32 @@ bool PragaProject::interpolationOutputPointsPeriod(QDate dateIni, QDate dateFin, varToSave.push_back(windVectorIntensity); } - errorString = "TODO"; + int nrDays = firstDate.daysTo(lastDate) + 1; + int nrDaysLoading = std::min(nrDays, 30); + int countDaysSaving = 0; + + QDate myDate = firstDate; + QDate lastLoadingDate = firstDate.addDays(nrDaysLoading - 1); + while (myDate <= lastDate) + { + countDaysSaving++; + + // check if load needed + if (myDate == firstDate || myDate > lastLoadingDate) + { + lastLoadingDate = myDate.addDays(nrDaysLoading - 1); + if (lastLoadingDate > lastDate) + lastLoadingDate = lastDate; + + // load one day before (for transmissivity) + logInfoGUI("Loading meteo points data from " + myDate.addDays(-1).toString("yyyy-MM-dd") + " to " + lastLoadingDate.toString("yyyy-MM-dd")); + if (! loadMeteoPointsData(myDate.addDays(-1), lastLoadingDate, isHourly, isDaily, false)) + return false; + } + } + + closeLogInfo(); + errorString = "TODO save data"; return false; } diff --git a/pragaProject/pragaProject.h b/pragaProject/pragaProject.h index c485cb83c..b5dba359f 100644 --- a/pragaProject/pragaProject.h +++ b/pragaProject/pragaProject.h @@ -107,7 +107,7 @@ bool downloadDailyDataArkimet(QList variables, bool prec0024, QDate startDate, QDate endDate, bool showInfo); bool downloadHourlyDataArkimet(QList variables, QDate startDate, QDate endDate, bool showInfo); - bool interpolationOutputPointsPeriod(QDate dateIni, QDate dateFin, QList variables); + bool interpolationOutputPointsPeriod(QDate dateIni, QDate lastDate, QList variables); bool interpolationMeteoGrid(meteoVariable myVar, frequencyType myFrequency, const Crit3DTime& myTime); bool interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QList variables, From c2c6ad9fb2291cd84b7af7b5fcdc55e103416401 Mon Sep 17 00:00:00 2001 From: ftomei Date: Wed, 29 Nov 2023 17:22:02 +0100 Subject: [PATCH 3/6] update --- pragaProject/pragaProject.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pragaProject/pragaProject.cpp b/pragaProject/pragaProject.cpp index 83145397a..15380fde0 100644 --- a/pragaProject/pragaProject.cpp +++ b/pragaProject/pragaProject.cpp @@ -1945,7 +1945,7 @@ bool PragaProject::interpolationOutputPointsPeriod(QDate firstDate, QDate lastDa } int nrDays = firstDate.daysTo(lastDate) + 1; - int nrDaysLoading = std::min(nrDays, 30); + int nrDaysLoading = std::min(nrDays, 10); int countDaysSaving = 0; QDate myDate = firstDate; From 9beebe9ac2a8fa6e317dd3832b63860a9495b048 Mon Sep 17 00:00:00 2001 From: ftomei Date: Wed, 29 Nov 2023 17:24:37 +0100 Subject: [PATCH 4/6] fix epoch --- project/project.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/project/project.cpp b/project/project.cpp index 94b62b815..ed775dea6 100644 --- a/project/project.cpp +++ b/project/project.cpp @@ -63,8 +63,8 @@ void Project::initializeProject() outputPointsDbHandler = nullptr; meteoGridDbHandler = nullptr; aggregationDbHandler = nullptr; - meteoPointsDbFirstTime.setTime_t(0); - meteoPointsDbLastTime.setTime_t(0); + meteoPointsDbFirstTime.setSecsSinceEpoch(0); + meteoPointsDbLastTime.setSecsSinceEpoch(0); meteoSettings->initialize(); quality->initialize(); @@ -1093,7 +1093,7 @@ bool Project::loadMeteoPointsDB(QString fileName) // find dates meteoPointsDbLastTime = findDbPointLastTime(); - meteoPointsDbFirstTime.setTime_t(0); + meteoPointsDbFirstTime.setSecsSinceEpoch(0); if (! meteoPointsDbLastTime.isNull()) { @@ -1190,7 +1190,7 @@ bool Project::loadAggregationDBAsMeteoPoints(QString fileName) // find dates meteoPointsDbLastTime = findDbPointLastTime(); - meteoPointsDbFirstTime.setTime_t(0); + meteoPointsDbFirstTime.setSecsSinceEpoch(0); if (! meteoPointsDbLastTime.isNull()) { From e8cc852099193147882010d13d276e04c6f9f840 Mon Sep 17 00:00:00 2001 From: ftomei Date: Thu, 30 Nov 2023 17:59:07 +0100 Subject: [PATCH 5/6] update interpolation output points --- dbMeteoGrid/dbMeteoGrid.cpp | 2 +- dbMeteoPoints/dbMeteoPointsHandler.cpp | 1 + gis/gis.cpp | 2 + interpolation/interpolation.cpp | 38 +++++------ interpolation/interpolation.h | 4 +- pragaProject/pragaMeteoMaps.cpp | 9 +++ pragaProject/pragaMeteoMaps.h | 2 + pragaProject/pragaProject.cpp | 94 ++++++++++++++++++++++++-- project/meteoMaps.cpp | 30 +++++--- project/meteoMaps.h | 2 + project/project.cpp | 4 ++ 11 files changed, 151 insertions(+), 37 deletions(-) diff --git a/dbMeteoGrid/dbMeteoGrid.cpp b/dbMeteoGrid/dbMeteoGrid.cpp index c2ac990b0..1f93361f5 100644 --- a/dbMeteoGrid/dbMeteoGrid.cpp +++ b/dbMeteoGrid/dbMeteoGrid.cpp @@ -2374,7 +2374,7 @@ std::vector Crit3DMeteoGridDbHandler::loadGridHourlyVar(QString *myError, QString tableH = _tableHourly.prefix + meteoPoint + _tableHourly.postFix; QDateTime dateTime, previousDateTime; dateTime.setTimeSpec(Qt::UTC); - previousDateTime.setTimeSpec((Qt::UTC)); + previousDateTime.setTimeSpec(Qt::UTC); std::vector hourlyVarList; diff --git a/dbMeteoPoints/dbMeteoPointsHandler.cpp b/dbMeteoPoints/dbMeteoPointsHandler.cpp index ca9414369..569c6b33a 100644 --- a/dbMeteoPoints/dbMeteoPointsHandler.cpp +++ b/dbMeteoPoints/dbMeteoPointsHandler.cpp @@ -224,6 +224,7 @@ QDateTime Crit3DMeteoPointsDbHandler::getFirstDate(frequencyType frequency) } } + QString tmp = firstDate.toString("yyyy-MM-dd"); return firstDate; } diff --git a/gis/gis.cpp b/gis/gis.cpp index 37aae1b10..b0beb0d32 100644 --- a/gis/gis.cpp +++ b/gis/gis.cpp @@ -371,6 +371,7 @@ namespace gis } + // clean the grid (all NO DATA) void Crit3DRasterGrid::emptyGrid() { for (int row = 0; row < header->nrRows; row++) @@ -378,6 +379,7 @@ namespace gis value[row][col] = header->flag; } + Crit3DRasterGrid::~Crit3DRasterGrid() { clear(); diff --git a/interpolation/interpolation.cpp b/interpolation/interpolation.cpp index 34f43cd4f..2e565091a 100644 --- a/interpolation/interpolation.cpp +++ b/interpolation/interpolation.cpp @@ -36,6 +36,7 @@ #include "meteoPoint.h" #include "gis.h" #include "spatialControl.h" +#include "interpolationPoint.h" #include "interpolation.h" #include @@ -1048,25 +1049,18 @@ void localSelection(vector &inputPoints, vector < mySettings.setLocalRadius(r1); } -bool checkPrecipitationZero(std::vector &myPoints, float precThreshold, int* nrPrecNotNull, bool* flatPrecipitation) + +bool checkPrecipitationZero(const std::vector &myPoints, float precThreshold, int &nrNotNull) { - *flatPrecipitation = true; - *nrPrecNotNull = 0; - float myValue = NODATA; + nrNotNull = 0; for (unsigned int i = 0; i < myPoints.size(); i++) if (myPoints[i].isActive) - if (int(myPoints[i].value) != int(NODATA)) - if (myPoints[i].value >= float(precThreshold)) - { - if (*nrPrecNotNull > 0 && myPoints[i].value != myValue) - *flatPrecipitation = false; + if (! isEqual(myPoints[i].value, NODATA)) + if (myPoints[i].value >= precThreshold) + nrNotNull++; - myValue = myPoints[i].value; - (*nrPrecNotNull)++; - } - - return (*nrPrecNotNull == 0); + return (nrNotNull == 0); } @@ -1643,8 +1637,7 @@ bool preInterpolation(std::vector &myPoints, Crit if (myVar == precipitation || myVar == dailyPrecipitation) { int nrPrecNotNull; - bool isFlatPrecipitation; - if (checkPrecipitationZero(myPoints, meteoSettings->getRainfallThreshold(), &nrPrecNotNull, &isFlatPrecipitation)) + if (checkPrecipitationZero(myPoints, meteoSettings->getRainfallThreshold(), nrPrecNotNull)) { mySettings->setPrecipitationAllZero(true); return true; @@ -1727,15 +1720,16 @@ float interpolate(vector &myPoints, Crit3DInterpo myResult = MAXVALUE(myResult, 0); return myResult; - } -bool getActiveProxyValues(Crit3DInterpolationSettings* mySettings, std::vector & allProxyValues, std::vector & activeProxyValues) + +bool getActiveProxyValues(Crit3DInterpolationSettings *mySettings, const std::vector &allProxyValues, std::vector &activeProxyValues) { Crit3DProxyCombination myCombination = mySettings->getCurrentCombination(); - std::vector myValues; - if (allProxyValues.size() != mySettings->getProxyNr()) return false; + if (allProxyValues.size() != mySettings->getProxyNr()) + return false; + activeProxyValues.clear(); bool isComplete = true; @@ -1744,12 +1738,14 @@ bool getActiveProxyValues(Crit3DInterpolationSettings* mySettings, std::vector < if (myCombination.getValue(i) && mySettings->getProxy(i)->getIsSignificant()) { activeProxyValues.push_back(allProxyValues[i]); - if (allProxyValues[i] == NODATA) isComplete = false; + if (allProxyValues[i] == NODATA) + isComplete = false; } return isComplete; } + void getProxyValuesXY(float x, float y, Crit3DInterpolationSettings* mySettings, std::vector &myValues) { float myValue; diff --git a/interpolation/interpolation.h b/interpolation/interpolation.h index 875354f18..4d1f9097b 100644 --- a/interpolation/interpolation.h +++ b/interpolation/interpolation.h @@ -53,13 +53,15 @@ bool krigLinearPrep(double *mySlope, double *myNugget, int nrPointData); void clearInterpolationPoints(); + bool checkPrecipitationZero(const std::vector &myPoints, float precThreshold, int &nrNotNull); bool neighbourhoodVariability(meteoVariable myVar, std::vector &myInterpolationPoints, Crit3DInterpolationSettings *mySettings, float x, float y, float z, int nMax, float* devSt, float* avgDeltaZ, float* minDistance); float interpolate(std::vector &myPoints, Crit3DInterpolationSettings *mySettings, Crit3DMeteoSettings *meteoSettings, meteoVariable myVar, float myX, float myY, float myZ, std::vector myProxyValues, bool excludeSupplemental); void getProxyValuesXY(float x, float y, Crit3DInterpolationSettings* mySettings, std::vector &myValues); - bool getActiveProxyValues(Crit3DInterpolationSettings* mySettings, std::vector &allProxyValues, std::vector &activeProxyValues); + + bool getActiveProxyValues(Crit3DInterpolationSettings *mySettings, const std::vector &allProxyValues, std::vector &activeProxyValues); void detrending(std::vector &myPoints, Crit3DProxyCombination myCombination, Crit3DInterpolationSettings *mySettings, Crit3DClimateParameters *myClimate, diff --git a/pragaProject/pragaMeteoMaps.cpp b/pragaProject/pragaMeteoMaps.cpp index 36a7bb91f..46298f9c9 100644 --- a/pragaProject/pragaMeteoMaps.cpp +++ b/pragaProject/pragaMeteoMaps.cpp @@ -38,6 +38,15 @@ void PragaHourlyMeteoMaps::clear() } +void PragaHourlyMeteoMaps::initialize() +{ + mapHourlyWindVectorInt->emptyGrid(); + mapHourlyWindVectorDir->emptyGrid(); + mapHourlyWindVectorX->emptyGrid(); + mapHourlyWindVectorY->emptyGrid(); +} + + gis::Crit3DRasterGrid* PragaHourlyMeteoMaps::getMapFromVar(meteoVariable myVar) { if (myVar == windVectorIntensity) diff --git a/pragaProject/pragaMeteoMaps.h b/pragaProject/pragaMeteoMaps.h index 36fde1986..cd54771e1 100644 --- a/pragaProject/pragaMeteoMaps.h +++ b/pragaProject/pragaMeteoMaps.h @@ -22,6 +22,8 @@ ~PragaHourlyMeteoMaps(); void clear(); + void initialize(); + gis::Crit3DRasterGrid* getMapFromVar(meteoVariable myVar); bool computeWindVector(); }; diff --git a/pragaProject/pragaProject.cpp b/pragaProject/pragaProject.cpp index 15380fde0..35069f500 100644 --- a/pragaProject/pragaProject.cpp +++ b/pragaProject/pragaProject.cpp @@ -1699,13 +1699,17 @@ QString getMapFileOutName(meteoVariable myVar, QDate myDate, int myHour) return name; } + gis::Crit3DRasterGrid* PragaProject::getPragaMapFromVar(meteoVariable myVar) { gis::Crit3DRasterGrid* myGrid = nullptr; + // search in project maps (hourlyMeteoMaps and radiationMaps) myGrid = getHourlyMeteoRaster(myVar); - if (myGrid == nullptr) myGrid = pragaHourlyMaps->getMapFromVar(myVar); + // search in pragaDailyMaps if (myGrid == nullptr) myGrid = pragaDailyMaps->getMapFromVar(myVar); + // saerch in pragaHourlyMaps + if (myGrid == nullptr) myGrid = pragaHourlyMaps->getMapFromVar(myVar); return myGrid; } @@ -1908,6 +1912,12 @@ bool PragaProject::interpolationOutputPointsPeriod(QDate firstDate, QDate lastDa return false; } + if (! DEM.isLoaded) + { + errorString = "Load a Digital Elevation Model before."; + return false; + } + // check dates if (firstDate.isNull() || lastDate.isNull() || firstDate > lastDate) { @@ -1944,8 +1954,34 @@ bool PragaProject::interpolationOutputPointsPeriod(QDate firstDate, QDate lastDa varToSave.push_back(windVectorIntensity); } + // initialize maps + if (pragaDailyMaps == nullptr) pragaDailyMaps = new Crit3DDailyMeteoMaps(DEM); + if (pragaHourlyMaps == nullptr) pragaHourlyMaps = new PragaHourlyMeteoMaps(DEM); + + // check if topographic distance is needed + bool useTd = false; + if (interpolationSettings.getUseTD()) + { + foreach (myVar, variables) + { + if (getUseTdVar(myVar)) + { + useTd = true; + break; + } + } + } + + // load topographic distance + if (useTd) + { + logInfoGUI("Loading topographic distance maps..."); + if (! loadTopographicDistanceMaps(true, false)) + return false; + } + int nrDays = firstDate.daysTo(lastDate) + 1; - int nrDaysLoading = std::min(nrDays, 10); + int nrDaysLoading = std::min(nrDays, 30); int countDaysSaving = 0; QDate myDate = firstDate; @@ -1966,6 +2002,57 @@ bool PragaProject::interpolationOutputPointsPeriod(QDate firstDate, QDate lastDa if (! loadMeteoPointsData(myDate.addDays(-1), lastLoadingDate, isHourly, isDaily, false)) return false; } + + if (isHourly) + { + // initialize + hourlyMeteoMaps->initialize(); + radiationMaps->initialize(); + pragaHourlyMaps->initialize(); + + for (int hour = 1; hour <= 24; hour++) + { + logInfoGUI("Interpolating hourly variables for " + myDate.toString("yyyy-MM-dd") + " " + QString::number(hour) + ":00"); + foreach (myVar, variables) + { + if (getVarFrequency(myVar) == hourly) + { + // TODO RH and other special variables + + if (! interpolationDemMain(myVar, getCrit3DTime(myDate, hour), getPragaMapFromVar(myVar))) + return false; + } + } + } + } + + if (isDaily) + { + logInfoGUI("Interpolating daily variables for " + myDate.toString("yyyy-MM-dd")); + + // initialize + pragaDailyMaps->initialize(); + + foreach (myVar, variables) + { + if (getVarFrequency(myVar) == daily) + { + // TODO special variables + + if (! interpolationDemMain(myVar, getCrit3DTime(myDate, 0), getPragaMapFromVar(myVar))) + return false; + + // fix daily temperatures consistency + if (myVar == dailyAirTemperatureMax || myVar == dailyAirTemperatureMin) + { + if (! pragaDailyMaps->fixDailyThermalConsistency()) + return false; + } + } + } + } + + myDate = myDate.addDays(1); } closeLogInfo(); @@ -2214,7 +2301,6 @@ bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QL } meteoGridDbHandler->meteoGrid()->spatialAggregateMeteoGrid(myVar, daily, getCrit3DDate(myDate), 0, 0, &DEM, myGrid, interpolationSettings.getMeteoGridAggrMethod()); - } } } @@ -2246,9 +2332,9 @@ bool PragaProject::interpolationMeteoGridPeriod(QDate dateIni, QDate dateFin, QL return false; return true; - } + bool PragaProject::interpolationMeteoGrid(meteoVariable myVar, frequencyType myFrequency, const Crit3DTime& myTime) { if (meteoGridDbHandler != nullptr) diff --git a/project/meteoMaps.cpp b/project/meteoMaps.cpp index c0366e6b5..d6a4b411d 100644 --- a/project/meteoMaps.cpp +++ b/project/meteoMaps.cpp @@ -58,6 +58,20 @@ void Crit3DDailyMeteoMaps::clear() } +void Crit3DDailyMeteoMaps::initialize() +{ + mapDailyTAvg->emptyGrid(); + mapDailyTMax->emptyGrid(); + mapDailyTMin->emptyGrid(); + mapDailyPrec->emptyGrid(); + mapDailyRHAvg->emptyGrid(); + mapDailyRHMin->emptyGrid(); + mapDailyRHMax->emptyGrid(); + mapDailyLeafW->emptyGrid(); + mapDailyET0HS->emptyGrid(); +} + + bool Crit3DDailyMeteoMaps::computeHSET0Map(gis::Crit3DGisSettings* gisSettings, Crit3DDate myDate) { float airTmin, airTmax; @@ -84,30 +98,26 @@ bool Crit3DDailyMeteoMaps::computeHSET0Map(gis::Crit3DGisSettings* gisSettings, return gis::updateMinMaxRasterGrid(mapDailyET0HS); } + bool Crit3DDailyMeteoMaps::fixDailyThermalConsistency() { if (! mapDailyTMax->isLoaded || ! mapDailyTMin->isLoaded) return true; if ( mapDailyTMax->getMapTime() != mapDailyTMin->getMapTime()) return true; - float TRange = NODATA; - unsigned row, col; - - for (row = 0; row < unsigned(mapDailyTMax->header->nrRows); row++) - for (col = 0; col < unsigned(mapDailyTMax->header->nrCols); col++) + for (unsigned row = 0; row < unsigned(mapDailyTMax->header->nrRows); row++) + { + for (unsigned col = 0; col < unsigned(mapDailyTMax->header->nrCols); col++) { if (! isEqual(mapDailyTMax->value[row][col], mapDailyTMax->header->flag) && ! isEqual(mapDailyTMin->value[row][col], mapDailyTMin->header->flag)) { if (mapDailyTMin->value[row][col] > mapDailyTMax->value[row][col]) { - //TRange = findNeighbourTRangeRaster(grdTmax, grdTmin, myRow, myCol) - if (! isEqual(TRange, NODATA)) - mapDailyTMin->value[row][col] = mapDailyTMax->value[row][col] - TRange; - else - mapDailyTMin->value[row][col] = mapDailyTMax->value[row][col] - 0.1f; + mapDailyTMin->value[row][col] = mapDailyTMax->value[row][col] - 0.1f; } } } + } return true; } diff --git a/project/meteoMaps.h b/project/meteoMaps.h index 4fd13786f..3ca912239 100644 --- a/project/meteoMaps.h +++ b/project/meteoMaps.h @@ -32,8 +32,10 @@ ~Crit3DDailyMeteoMaps(); void clear(); + void initialize(); gis::Crit3DRasterGrid* getMapFromVar(meteoVariable myVar); + bool computeHSET0Map(gis::Crit3DGisSettings *gisSettings, Crit3DDate myDate); bool fixDailyThermalConsistency(); }; diff --git a/project/project.cpp b/project/project.cpp index ed775dea6..6326cda33 100644 --- a/project/project.cpp +++ b/project/project.cpp @@ -63,6 +63,9 @@ void Project::initializeProject() outputPointsDbHandler = nullptr; meteoGridDbHandler = nullptr; aggregationDbHandler = nullptr; + + meteoPointsDbFirstTime.setTimeSpec(Qt::UTC); + meteoPointsDbLastTime.setTimeSpec(Qt::UTC); meteoPointsDbFirstTime.setSecsSinceEpoch(0); meteoPointsDbLastTime.setSecsSinceEpoch(0); @@ -2228,6 +2231,7 @@ bool Project::interpolationDem(meteoVariable myVar, const Crit3DTime& myTime, gi } myRaster->setMapTime(myTime); + return true; } From 1854bd7f978043c7c2c7b8e4dbd7456bdd57447f Mon Sep 17 00:00:00 2001 From: ftomei Date: Thu, 30 Nov 2023 20:00:08 +0100 Subject: [PATCH 6/6] fix interpolation output points --- dbMeteoPoints/dbMeteoPointsHandler.cpp | 1 - pragaProject/pragaProject.cpp | 36 ++++++++++++-- project/meteoMaps.cpp | 8 +-- project/meteoMaps.h | 2 +- project/project.cpp | 29 +++++++---- solarRadiation/solarRadiation.cpp | 67 ++++++++++---------------- solarRadiation/solarRadiation.h | 8 +-- 7 files changed, 89 insertions(+), 62 deletions(-) diff --git a/dbMeteoPoints/dbMeteoPointsHandler.cpp b/dbMeteoPoints/dbMeteoPointsHandler.cpp index 569c6b33a..ca9414369 100644 --- a/dbMeteoPoints/dbMeteoPointsHandler.cpp +++ b/dbMeteoPoints/dbMeteoPointsHandler.cpp @@ -224,7 +224,6 @@ QDateTime Crit3DMeteoPointsDbHandler::getFirstDate(frequencyType frequency) } } - QString tmp = firstDate.toString("yyyy-MM-dd"); return firstDate; } diff --git a/pragaProject/pragaProject.cpp b/pragaProject/pragaProject.cpp index 35069f500..5a7cccd0d 100644 --- a/pragaProject/pragaProject.cpp +++ b/pragaProject/pragaProject.cpp @@ -2017,9 +2017,33 @@ bool PragaProject::interpolationOutputPointsPeriod(QDate firstDate, QDate lastDa { if (getVarFrequency(myVar) == hourly) { - // TODO RH and other special variables + setComputeOnlyPoints(true); - if (! interpolationDemMain(myVar, getCrit3DTime(myDate, hour), getPragaMapFromVar(myVar))) + // TODO other special variables + + bool isOk; + if (myVar == airRelHumidity && interpolationSettings.getUseDewPoint()) + { + if (interpolationSettings.getUseInterpolatedTForRH()) + { + passInterpolatedTemperatureToHumidityPoints(getCrit3DTime(myDate, hour), meteoSettings); + } + + isOk = interpolationDemMain(airDewTemperature, getCrit3DTime(myDate, hour), hourlyMeteoMaps->mapHourlyTdew); + + if (isOk) + { + hourlyMeteoMaps->computeRelativeHumidityMap(hourlyMeteoMaps->mapHourlyRelHum); + } + } + else + { + isOk = interpolationDemMain(myVar, getCrit3DTime(myDate, hour), getPragaMapFromVar(myVar)); + } + + setComputeOnlyPoints(false); + + if (! isOk) return false; } } @@ -2037,9 +2061,15 @@ bool PragaProject::interpolationOutputPointsPeriod(QDate firstDate, QDate lastDa { if (getVarFrequency(myVar) == daily) { + setComputeOnlyPoints(true); + // TODO special variables - if (! interpolationDemMain(myVar, getCrit3DTime(myDate, 0), getPragaMapFromVar(myVar))) + bool isOk = interpolationDemMain(myVar, getCrit3DTime(myDate, 0), getPragaMapFromVar(myVar)); + + setComputeOnlyPoints(false); + + if (! isOk) return false; // fix daily temperatures consistency diff --git a/project/meteoMaps.cpp b/project/meteoMaps.cpp index d6a4b411d..f790088ac 100644 --- a/project/meteoMaps.cpp +++ b/project/meteoMaps.cpp @@ -295,10 +295,10 @@ bool Crit3DHourlyMeteoMaps::computeLeafWetnessMap() -bool Crit3DHourlyMeteoMaps::computeRelativeHumidityMap(gis::Crit3DRasterGrid* myGrid) +bool Crit3DHourlyMeteoMaps::computeRelativeHumidityMap(gis::Crit3DRasterGrid* myRaster) { float airT, dewT; - myGrid->emptyGrid(); + myRaster->emptyGrid(); for (long row = 0; row < mapHourlyRelHum->header->nrRows ; row++) { @@ -309,12 +309,12 @@ bool Crit3DHourlyMeteoMaps::computeRelativeHumidityMap(gis::Crit3DRasterGrid* my if (! isEqual(airT, mapHourlyTair->header->flag) && ! isEqual(dewT, mapHourlyTdew->header->flag)) { - myGrid->value[row][col] = relHumFromTdew(dewT, airT); + myRaster->value[row][col] = relHumFromTdew(dewT, airT); } } } - return gis::updateMinMaxRasterGrid(myGrid); + return gis::updateMinMaxRasterGrid(myRaster); } diff --git a/project/meteoMaps.h b/project/meteoMaps.h index 3ca912239..a8e3150a7 100644 --- a/project/meteoMaps.h +++ b/project/meteoMaps.h @@ -63,7 +63,7 @@ gis::Crit3DRasterGrid* getMapFromVar(meteoVariable myVar); bool computeET0PMMap(const gis::Crit3DRasterGrid &DEM, Crit3DRadiationMaps *radMaps); - bool computeRelativeHumidityMap(gis::Crit3DRasterGrid* myGrid); + bool computeRelativeHumidityMap(gis::Crit3DRasterGrid* myRaster); bool computeLeafWetnessMap(); void setComputed(bool value); bool getComputed(); diff --git a/project/project.cpp b/project/project.cpp index 6326cda33..deed22571 100644 --- a/project/project.cpp +++ b/project/project.cpp @@ -2089,9 +2089,14 @@ bool Project::interpolationOutputPoints(std::vector getRowCol(x, y, row, col); if (! gis::isOutOfGridRowCol(row, col, *outputGrid)) { - if (getUseDetrendingVar(myVar)) getProxyValuesXY(x, y, &interpolationSettings, proxyValues); + if (getUseDetrendingVar(myVar)) + { + getProxyValuesXY(x, y, &interpolationSettings, proxyValues); + } + outputPoints[i].currentValue = interpolate(interpolationPoints, &interpolationSettings, - meteoSettings, myVar, x, y, z, proxyValues, true); + meteoSettings, myVar, x, y, z, proxyValues, true); + outputGrid->value[row][col] = outputPoints[i].currentValue; } } @@ -2100,6 +2105,7 @@ bool Project::interpolationOutputPoints(std::vector initialize(); @@ -2318,7 +2324,13 @@ bool Project::interpolationDemLocalDetrending(meteoVariable myVar, const Crit3DT bool Project::interpolateDemRadiation(const Crit3DTime& myTime, gis::Crit3DRasterGrid *myRaster) { - this->radiationMaps->initialize(); + if (! radiationMaps->latMap->isLoaded || ! radiationMaps->lonMap->isLoaded) + return false; + + if (! radiationMaps->slopeMap->isLoaded || ! radiationMaps->aspectMap->isLoaded) + return false; + + radiationMaps->initialize(); std::vector interpolationPoints; @@ -2334,7 +2346,7 @@ bool Project::interpolateDemRadiation(const Crit3DTime& myTime, gis::Crit3DRaste // TODO: add flag to parameters. Could be NOT wanted if (! computeTransmissivityFromTRange(meteoPoints, nrMeteoPoints, myTime)) { - logError("Function interpolateDemRadiation: cannot compute transmissivity."); + logError("Error in function interpolateDemRadiation: cannot compute transmissivity."); return false; } } @@ -2345,7 +2357,7 @@ bool Project::interpolateDemRadiation(const Crit3DTime& myTime, gis::Crit3DRaste meteoSettings, &climateParameters, interpolationPoints, checkSpatialQuality); if (! result) { - logError("Function interpolateDemRadiation: not enough transmissivity data."); + logError("Error in function interpolateDemRadiation: not enough transmissivity data."); return false; } @@ -2355,7 +2367,6 @@ bool Project::interpolateDemRadiation(const Crit3DTime& myTime, gis::Crit3DRaste // interpolate transmissivity if (getComputeOnlyPoints()) { - radiationMaps->transmissivityMap->initializeGrid(DEM); result = interpolationOutputPoints(interpolationPoints, radiationMaps->transmissivityMap, atmTransmissivity); } else @@ -2365,18 +2376,18 @@ bool Project::interpolateDemRadiation(const Crit3DTime& myTime, gis::Crit3DRaste } if (! result) { - logError("Function interpolateDemRadiation: error interpolating transmissivity."); + logError("Function interpolateDemRadiation: error in interpolating transmissivity."); return false; } // compute radiation if (getComputeOnlyPoints()) { - result = radiation::computeRadiationOutputPoints(&radSettings, this->DEM, this->radiationMaps, outputPoints, myTime); + result = radiation::computeRadiationOutputPoints(&radSettings, DEM, radiationMaps, outputPoints, myTime); } else { - result = radiation::computeRadiationGrid(&radSettings, this->DEM, this->radiationMaps, myTime); + result = radiation::computeRadiationDEM(&radSettings, DEM, radiationMaps, myTime); } if (! result) { diff --git a/solarRadiation/solarRadiation.cpp b/solarRadiation/solarRadiation.cpp index 3c59b5683..6f3209526 100644 --- a/solarRadiation/solarRadiation.cpp +++ b/solarRadiation/solarRadiation.cpp @@ -824,7 +824,7 @@ bool computeRadiationRsun(Crit3DRadiationSettings* radSettings, float temperatur } - bool computeRadiationRSunGridPoint(Crit3DRadiationSettings* radSettings, const gis::Crit3DRasterGrid& dem, + bool computeRadiationDemPoint(Crit3DRadiationSettings* radSettings, const gis::Crit3DRasterGrid& myDem, Crit3DRadiationMaps* radiationMaps, TradPoint radPoint, int row, int col, const Crit3DTime& myTime) { @@ -839,8 +839,8 @@ bool computeRadiationRsun(Crit3DRadiationSettings* radSettings, float temperatur float transmissivity = radiationMaps->transmissivityMap->value[row][col]; TsunPosition sunPosition; - if (!computeRadiationRsun(radSettings, TEMPERATURE_DEFAULT, PRESSURE_SEALEVEL, myTime, - linke, albedo, radSettings->getClearSky(), transmissivity, &sunPosition, &radPoint, dem)) + if (! computeRadiationRsun(radSettings, TEMPERATURE_DEFAULT, PRESSURE_SEALEVEL, myTime, + linke, albedo, radSettings->getClearSky(), transmissivity, &sunPosition, &radPoint, myDem)) return false; /* @@ -857,6 +857,7 @@ bool computeRadiationRsun(Crit3DRadiationSettings* radSettings, float temperatur return true; } + bool computeRadiationPotentialRSunMeteoPoint(Crit3DRadiationSettings* radSettings, const gis::Crit3DRasterGrid& dem, Crit3DMeteoPoint* myMeteoPoint, float slope, float aspect, const Crit3DTime& myTime, TradPoint* radPoint) { @@ -905,27 +906,38 @@ bool computeRadiationRsun(Crit3DRadiationSettings* radSettings, float temperatur } - bool computeRadiationGridRsun(Crit3DRadiationSettings* radSettings, const gis::Crit3DRasterGrid& dem, - Crit3DRadiationMaps* radiationMaps, const Crit3DTime& myTime) + bool computeRadiationDEM(Crit3DRadiationSettings* radSettings, const gis::Crit3DRasterGrid& myDem, + Crit3DRadiationMaps* radiationMaps, const Crit3DTime& myTime) { + if (radSettings->getAlgorithm() != RADIATION_ALGORITHM_RSUN) + return false; + int row, col; TradPoint radPoint; - for (row = 0; row < dem.header->nrRows; row++ ) + for (row = 0; row < myDem.header->nrRows; row++ ) { - for (col = 0; col < dem.header->nrCols; col++) + for (col = 0; col < myDem.header->nrCols; col++) { - if(isGridPointComputable(radSettings, row, col, dem, radiationMaps)) + if(isGridPointComputable(radSettings, row, col, myDem, radiationMaps)) { - dem.getXY(row, col, radPoint.x, radPoint.y); - radPoint.height = dem.value[row][col]; + myDem.getXY(row, col, radPoint.x, radPoint.y); + radPoint.height = myDem.value[row][col]; - if (! computeRadiationRSunGridPoint(radSettings, dem, radiationMaps, radPoint, row, col, myTime)) + if (! computeRadiationDemPoint(radSettings, myDem, radiationMaps, radPoint, row, col, myTime)) return false; } } } + updateRadiationMaps(radiationMaps, myTime); + + return true; + } + + + void updateRadiationMaps(Crit3DRadiationMaps* radiationMaps, const Crit3DTime &myTime) + { gis::updateMinMaxRasterGrid(radiationMaps->sunElevationMap); gis::updateMinMaxRasterGrid(radiationMaps->transmissivityMap); gis::updateMinMaxRasterGrid(radiationMaps->beamRadiationMap); @@ -941,7 +953,6 @@ bool computeRadiationRsun(Crit3DRadiationSettings* radSettings, float temperatur radiationMaps->globalRadiationMap->setMapTime(myTime); radiationMaps->setComputed(true); - return true; } @@ -1080,37 +1091,10 @@ bool computeRadiationRsun(Crit3DRadiationSettings* radSettings, float temperatur } - bool preConditionsRadiationGrid(const Crit3DRadiationMaps& radiationMaps) - { - if (! radiationMaps.latMap->isLoaded || ! radiationMaps.lonMap->isLoaded) return false; - if (! radiationMaps.slopeMap->isLoaded || ! radiationMaps.aspectMap->isLoaded) return false; - - return true; - } - - - bool computeRadiationGrid(Crit3DRadiationSettings* radSettings, const gis::Crit3DRasterGrid& dem, - Crit3DRadiationMaps* radiationMaps, const Crit3DTime& myTime) - { - if (! preConditionsRadiationGrid(*radiationMaps)) - return false; - - if (radSettings->getAlgorithm() == RADIATION_ALGORITHM_RSUN) - { - return computeRadiationGridRsun(radSettings, dem, radiationMaps, myTime); - } - else - return false; - } - - bool computeRadiationOutputPoints(Crit3DRadiationSettings *radSettings, const gis::Crit3DRasterGrid& dem, Crit3DRadiationMaps* radiationMaps, std::vector &outputPoints, const Crit3DTime& myTime) { - if (! preConditionsRadiationGrid(*radiationMaps)) - return false; - if (radSettings->getAlgorithm() != RADIATION_ALGORITHM_RSUN) return false; @@ -1127,13 +1111,14 @@ bool computeRadiationRsun(Crit3DRadiationSettings* radSettings, float temperatur if(isGridPointComputable(radSettings, row, col, dem, radiationMaps)) { - if (! computeRadiationRSunGridPoint(radSettings, dem, radiationMaps, radPoint, row, col, myTime)) + if (! computeRadiationDemPoint(radSettings, dem, radiationMaps, radPoint, row, col, myTime)) return false; } } } - radiationMaps->setComputed(true); + updateRadiationMaps(radiationMaps, myTime); + return true; } diff --git a/solarRadiation/solarRadiation.h b/solarRadiation/solarRadiation.h index 75dc79cc2..e7b9c7af7 100644 --- a/solarRadiation/solarRadiation.h +++ b/solarRadiation/solarRadiation.h @@ -66,7 +66,10 @@ float myLinke, float myAlbedo, float myClearSkyTransmissivity, float transmissivity, TsunPosition* mySunPosition, TradPoint* myPoint, const gis::Crit3DRasterGrid& myDEM); - bool computeRadiationRSunGridPoint(Crit3DRadiationSettings* mySettings, const gis::Crit3DRasterGrid& myDEM, + bool computeRadiationDEM(Crit3DRadiationSettings *radSettings, const gis::Crit3DRasterGrid& myDem, + Crit3DRadiationMaps* radiationMaps, const Crit3DTime& myTime); + + bool computeRadiationDemPoint(Crit3DRadiationSettings* mySettings, const gis::Crit3DRasterGrid& myDem, Crit3DRadiationMaps* radiationMaps, TradPoint radPoint, int row, int col, const Crit3DTime& myTime); @@ -74,8 +77,7 @@ Crit3DRadiationMaps *radiationMaps, std::vector &outputPoints, const Crit3DTime& myCrit3DTime); - bool computeRadiationGrid(Crit3DRadiationSettings *radSettings, const gis::Crit3DRasterGrid& myDEM, - Crit3DRadiationMaps* radiationMaps, const Crit3DTime& myCrit3DTime); + void updateRadiationMaps(Crit3DRadiationMaps* radiationMaps, const Crit3DTime &myTime); float computePointTransmissivity(Crit3DRadiationSettings *mySettings, const gis::Crit3DPoint& myPoint, Crit3DTime myTime, float* measuredRad, int windowWidth, int timeStepSecond, const gis::Crit3DRasterGrid& myDEM);