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

…e datasets with irregular vertical levels
qgis · Jan 17, 2020 · 830e080 · 830e080
1 parent f780b59
commit 830e080
Show file tree

Hide file tree

Showing 3 changed files with 141 additions and 67 deletions.
diff --git a/src/core/mesh/qgsmesh3daveraging.cpp b/src/core/mesh/qgsmesh3daveraging.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 )

diff --git a/src/core/mesh/qgsmesh3daveraging.h b/src/core/mesh/qgsmesh3daveraging.h
@@ -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
      */

diff --git a/tests/src/core/testqgsmesh3daveraging.cpp b/tests/src/core/testqgsmesh3daveraging.cpp
@@ -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"