Skip to content

Commit

Permalink
fix interpolation output points
Browse files Browse the repository at this point in the history
  • Loading branch information
@ftomei
ftomei committed Nov 30, 2023
1 parent e8cc852 commit 1854bd7
Show file tree
Hide file tree
Showing 7 changed files with 89 additions and 62 deletions.
1 change: 0 additions & 1 deletion dbMeteoPoints/dbMeteoPointsHandler.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -224,7 +224,6 @@ QDateTime Crit3DMeteoPointsDbHandler::getFirstDate(frequencyType frequency)
}
}

QString tmp = firstDate.toString("yyyy-MM-dd");
return firstDate;
}

Expand Down
36 changes: 33 additions & 3 deletions pragaProject/pragaProject.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}
}
Expand All @@ -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
Expand Down
8 changes: 4 additions & 4 deletions project/meteoMaps.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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++)
{
Expand All @@ -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);
}


Expand Down
2 changes: 1 addition & 1 deletion project/meteoMaps.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand Down
29 changes: 20 additions & 9 deletions project/project.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2089,9 +2089,14 @@ bool Project::interpolationOutputPoints(std::vector <Crit3DInterpolationDataPoin
outputGrid->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;
}
}
Expand All @@ -2100,6 +2105,7 @@ bool Project::interpolationOutputPoints(std::vector <Crit3DInterpolationDataPoin
return true;
}


bool Project::computeStatisticsCrossValidation(Crit3DTime myTime, meteoVariable myVar, crossValidationStatistics* myStats)
{
myStats->initialize();
Expand Down Expand Up @@ -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 <Crit3DInterpolationDataPoint> interpolationPoints;

Expand All @@ -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;
}
}
Expand All @@ -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;
}

Expand All @@ -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
Expand All @@ -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)
{
Expand Down
67 changes: 26 additions & 41 deletions solarRadiation/solarRadiation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
{
Expand All @@ -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;

/*
Expand All @@ -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)
{
Expand Down Expand Up @@ -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);
Expand All @@ -941,7 +953,6 @@ bool computeRadiationRsun(Crit3DRadiationSettings* radSettings, float temperatur
radiationMaps->globalRadiationMap->setMapTime(myTime);

radiationMaps->setComputed(true);
return true;
}


Expand Down Expand Up @@ -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<gis::Crit3DOutputPoint> &outputPoints,
const Crit3DTime& myTime)
{
if (! preConditionsRadiationGrid(*radiationMaps))
return false;

if (radSettings->getAlgorithm() != RADIATION_ALGORITHM_RSUN)
return false;

Expand All @@ -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;
}

Expand Down
8 changes: 5 additions & 3 deletions solarRadiation/solarRadiation.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,16 +66,18 @@
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);

bool computeRadiationOutputPoints(Crit3DRadiationSettings *radSettings, const gis::Crit3DRasterGrid& myDEM,
Crit3DRadiationMaps *radiationMaps, std::vector<gis::Crit3DOutputPoint> &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);
Expand Down

0 comments on commit 1854bd7

Please sign in to comment.