Skip to content
Permalink
Browse files

fix bug when using elevation/depth/height averaging 3d method for som…

…e datasets with irregular vertical levels
  • Loading branch information
PeterPetrik committed Jan 17, 2020
1 parent f780b59 commit 830e080865d499bf354f31e737a258e36041cccf
Showing with 141 additions and 67 deletions.
  1. +90 −67 src/core/mesh/qgsmesh3daveraging.cpp
  2. +15 −0 src/core/mesh/qgsmesh3daveraging.h
  3. +36 −0 tests/src/core/testqgsmesh3daveraging.cpp
@@ -94,85 +94,35 @@ QgsMeshDataBlock QgsMesh3dAveragingMethod::calculate( const QgsMesh3dDataBlock &
std::swap( faceLevelTop, faceLevelBottom );
}

double totalAveragedHeight = 0;
double nSumX = 0.0;
double nSumY = 0.0;

double methodLevelTop = std::numeric_limits<double>::quiet_NaN();
double methodLevelBottom = std::numeric_limits<double>::quiet_NaN();

volumeRangeForFace( methodLevelTop,
methodLevelBottom,
verticalLevelsForFace );

if ( std::isnan( methodLevelTop ) || std::isnan( methodLevelBottom ) )
continue;

// the level is value below surface, so top level (-0.1m) is usually higher number than bottom level (e.g. -1.2m)
if ( methodLevelTop < methodLevelBottom )
{
std::swap( methodLevelTop, methodLevelBottom );
}

// check if we are completely outside the limits
if ( methodLevelTop < faceLevelBottom )
continue;
if ( methodLevelBottom > faceLevelTop )
continue;

// Now go through all volumes below face and check if we need to take that volume into consideration
for ( int relativeVolumeIndex = 0; relativeVolumeIndex < volumesBelowFaceCount; ++relativeVolumeIndex )
if ( !std::isnan( methodLevelTop ) && !std::isnan( methodLevelBottom ) )
{
const int volumeIndex = startVolumeIndex + relativeVolumeIndex;
double volumeLevelTop = verticalLevelsForFace[relativeVolumeIndex];
double volumeLevelBottom = verticalLevelsForFace[relativeVolumeIndex + 1];
if ( volumeLevelTop < volumeLevelBottom )
{
std::swap( volumeLevelTop, volumeLevelBottom );
}

const double intersectionLevelTop = std::min( methodLevelTop, volumeLevelTop );
const double intersectionLevelBottom = std::max( methodLevelBottom, volumeLevelBottom );
const double effectiveInterval = intersectionLevelTop - intersectionLevelBottom;

if ( effectiveInterval > eps )
// the level is value below surface, so top level (-0.1m) is usually higher number than bottom level (e.g. -1.2m)
if ( methodLevelTop < methodLevelBottom )
{
if ( isVector )
{
const double x = volumeValues[2 * volumeIndex ];
const double y = volumeValues[ 2 * volumeIndex + 1 ];
if ( ! std::isnan( x ) &&
! std::isnan( y )
)
{
nSumX += x * effectiveInterval;
nSumY += y * effectiveInterval;
totalAveragedHeight += effectiveInterval;
}
}
else
{
const double x = volumeValues[ volumeIndex ];
if ( ! std::isnan( x ) )
{
nSumX += x * effectiveInterval;
totalAveragedHeight += effectiveInterval;
}
}
std::swap( methodLevelTop, methodLevelBottom );
}
}

// calculate average
if ( totalAveragedHeight > eps )
{
if ( isVector )
// check if we are completely outside the limits
if ( ( methodLevelTop >= faceLevelBottom ) && ( methodLevelBottom <= faceLevelTop ) )
{
valuesFaces[2 * faceIndex] = nSumX / totalAveragedHeight;
valuesFaces[2 * faceIndex + 1 ] = nSumY / totalAveragedHeight;
}
else
{
valuesFaces[faceIndex] = nSumX / totalAveragedHeight;
averageVolumeValuesForFace(
faceIndex,
volumesBelowFaceCount,
startVolumeIndex,
methodLevelTop,
methodLevelBottom,
isVector,
verticalLevelsForFace,
volumeValues,
valuesFaces
);
}
}

@@ -188,6 +138,79 @@ QgsMesh3dAveragingMethod::Method QgsMesh3dAveragingMethod::method() const
return mMethod;
}

void QgsMesh3dAveragingMethod::averageVolumeValuesForFace(
int faceIndex,
int volumesBelowFaceCount,
int startVolumeIndex,
double methodLevelTop,
double methodLevelBottom,
bool isVector,
const QVector<double> &verticalLevelsForFace,
const QVector<double> &volumeValues,
QVector<double> &valuesFaces
) const
{
double totalAveragedHeight = 0;
double nSumX = 0.0;
double nSumY = 0.0;

// Now go through all volumes below face and check if we need to take that volume into consideration
for ( int relativeVolumeIndex = 0; relativeVolumeIndex < volumesBelowFaceCount; ++relativeVolumeIndex )
{
const int volumeIndex = startVolumeIndex + relativeVolumeIndex;
double volumeLevelTop = verticalLevelsForFace[relativeVolumeIndex];
double volumeLevelBottom = verticalLevelsForFace[relativeVolumeIndex + 1];
if ( volumeLevelTop < volumeLevelBottom )
{
std::swap( volumeLevelTop, volumeLevelBottom );
}

const double intersectionLevelTop = std::min( methodLevelTop, volumeLevelTop );
const double intersectionLevelBottom = std::max( methodLevelBottom, volumeLevelBottom );
const double effectiveInterval = intersectionLevelTop - intersectionLevelBottom;

if ( effectiveInterval > eps )
{
if ( isVector )
{
const double x = volumeValues[2 * volumeIndex ];
const double y = volumeValues[ 2 * volumeIndex + 1 ];
if ( ! std::isnan( x ) &&
! std::isnan( y )
)
{
nSumX += x * effectiveInterval;
nSumY += y * effectiveInterval;
totalAveragedHeight += effectiveInterval;
}
}
else
{
const double x = volumeValues[ volumeIndex ];
if ( ! std::isnan( x ) )
{
nSumX += x * effectiveInterval;
totalAveragedHeight += effectiveInterval;
}
}
}
}

// calculate average
if ( totalAveragedHeight > eps )
{
if ( isVector )
{
valuesFaces[2 * faceIndex] = nSumX / totalAveragedHeight;
valuesFaces[2 * faceIndex + 1 ] = nSumY / totalAveragedHeight;
}
else
{
valuesFaces[faceIndex] = nSumX / totalAveragedHeight;
}
}
}

bool QgsMesh3dAveragingMethod::equals( const QgsMesh3dAveragingMethod *a, const QgsMesh3dAveragingMethod *b )
{
if ( a )
@@ -119,6 +119,21 @@ class CORE_EXPORT QgsMesh3dAveragingMethod SIP_ABSTRACT
//! Returns whether the method is correctly initialized
virtual bool hasValidInputs() const = 0;

/**
* For one face, Calculates average of volume values
*/
void averageVolumeValuesForFace(
int faceIndex,
int volumesBelowFaceCount,
int startVolumeIndex,
double methodLevelTop,
double methodLevelBottom,
bool isVector,
const QVector<double> &verticalLevelsForFace,
const QVector<double> &volumeValues,
QVector<double> &valuesFaces
) const;

/**
* For one face, calculates absolute vertical levels between which we need to average
*/
@@ -61,6 +61,7 @@ class TestQgsMesh3dAveraging: public QObject

void testMeshElevationAveragingMethod_data();
void testMeshElevationAveragingMethod();
void testMeshElevationAveragingMethodVariableMesh();

private:
void compare( const QgsMesh3dAveragingMethod *method, double expected, bool valid );
@@ -328,5 +329,40 @@ void TestQgsMesh3dAveraging::testMeshElevationAveragingMethod()
compare( &method, expected, startParam <= 0 );
}

void TestQgsMesh3dAveraging::testMeshElevationAveragingMethodVariableMesh()
{
// Test the situation when the number of vertical levels is different for
// each face and also that for face 1 the vertical levels are outside of
// requested elevation range
QVector<int> faceToVolumeIndex = { 0, 4, 7 };
QVector<double> verticalLevels = { -1.0, -2.0, -3.0, -4.0, -5.0,
-1.0, -1.1, -1.3, -1.5,
-1.0, -2.0, -3.0, -4.0, -5.0, -6.0
};

QVector<int> verticalLevelsCount = { 4, 3, 5 };

QgsMesh3dDataBlock scalarBlock2( 3, false );
QVector<double> values = { 1, 2, 3, 4,
0, 0, 0,
100, 200, 300, 400, 500
};

scalarBlock2.setFaceToVolumeIndex( faceToVolumeIndex );
scalarBlock2.setVerticalLevelsCount( verticalLevelsCount );
scalarBlock2.setVerticalLevels( verticalLevels );
scalarBlock2.setValues( values );
scalarBlock2.setValid( true );

QgsMeshElevationAveragingMethod method( -2.0, -6.0 );
QgsMeshDataBlock block = method.calculate( scalarBlock2 );
QVERIFY( block.isValid() );
QVERIFY( block.count() == 3 );

QCOMPARE( block.value( 0 ).scalar(), ( 2 + 3 + 4 ) / 3.0 );
QVERIFY( std::isnan( block.value( 1 ).scalar() ) );
QCOMPARE( block.value( 2 ).scalar(), ( 200 + 300 + 400 + 500 ) / 4.0 );
}

QGSTEST_MAIN( TestQgsMesh3dAveraging )
#include "testqgsmesh3daveraging.moc"

0 comments on commit 830e080

Please sign in to comment.
You can’t perform that action at this time.