[feature] [mesh] fix #31550 Expose new methods in analysis library for exporting contours for mesh layers

The mesh styling dialog now has new combo box for selection of resampling method for datesets defined on faces.
Also the analysis library has new class QgsMeshContours with export function for contour lines and contour polygons.

python/analysis/analysis_auto.sip

@@ -15,6 +15,7 @@
 %Include auto_generated/raster/qgsrastercalcnode.sip
 %Include auto_generated/raster/qgstotalcurvaturefilter.sip
 %Include auto_generated/mesh/qgsmeshcalculator.sip
+%Include auto_generated/mesh/qgsmeshcontours.sip
 %Include auto_generated/vector/qgsgeometrysnapper.sip
 %Include auto_generated/vector/qgsgeometrysnappersinglesource.sip
 %Include auto_generated/vector/qgszonalstatistics.sip

python/analysis/auto_generated/mesh/qgsmeshcontours.sip.in

@@ -0,0 +1,77 @@
+/************************************************************************
+ * This file has been generated automatically from                      *
+ *                                                                      *
+ * src/analysis/mesh/qgsmeshcontours.h                                  *
+ *                                                                      *
+ * Do not edit manually ! Edit header and run scripts/sipify.pl again   *
+ ************************************************************************/
+
+
+
+
+
+
+class QgsMeshContours
+{
+%Docstring
+
+Exporter of contours lines or polygons from a mesh layer.
+
+.. versionadded:: 3.12
+%End
+
+%TypeHeaderCode
+#include "qgsmeshcontours.h"
+%End
+  public:
+
+    QgsMeshContours( QgsMeshLayer *layer );
+%Docstring
+Constructs the mesh contours exporter.
+Caches the native and triangular mesh from data provider
+
+:param layer: mesh layer to be associated with this exporter
+%End
+    ~QgsMeshContours();
+
+    QgsGeometry exportLines( const QgsMeshDatasetIndex &index,
+                             double value,
+                             QgsMeshRendererScalarSettings::DataInterpolationMethod method,
+                             QgsFeedback *feedback = 0 );
+%Docstring
+Exports multi line string containing the contour line for particular dataset and value
+
+:param index: dataset index
+:param value: value of the contour line
+:param method: for datasets defined on faces, the method will be used to convert data to vertices
+:param feedback: optional feedback object for progress and cancellation support
+
+:return: MultiLineString geometry containing contour lines
+%End
+
+    QgsGeometry exportPolygons( const QgsMeshDatasetIndex &index,
+                                double min_value,
+                                double max_value,
+                                QgsMeshRendererScalarSettings::DataInterpolationMethod method,
+                                QgsFeedback *feedback = 0 );
+%Docstring
+Exports multi polygons representing the areas with values in range for particular dataset
+
+:param index: dataset index
+:param min_value: minimum of the value interval for contour polygon
+:param max_value: maximum of the value interval for contour polygon
+:param method: for datasets defined on faces, the method will be used to convert data to vertices
+:param feedback: optional feedback object for progress and cancellation support
+
+:return: MultiPolygon geometry containing contour polygons
+%End
+
+};
+
+/************************************************************************
+ * This file has been generated automatically from                      *
+ *                                                                      *
+ * src/analysis/mesh/qgsmeshcontours.h                                  *
+ *                                                                      *
+ * Do not edit manually ! Edit header and run scripts/sipify.pl again   *
+ ************************************************************************/

python/core/auto_generated/mesh/qgsmeshrenderersettings.sip.in

@@ -83,6 +83,14 @@ Represents a mesh renderer settings for scalar datasets
 #include "qgsmeshrenderersettings.h"
 %End
   public:
+    enum DataInterpolationMethod
+    {
+
+      None,
+
+      NeighbourAverage,
+    };
+
     QgsColorRampShader colorRampShader() const;
 %Docstring
 Returns color ramp shader function
@@ -112,6 +120,22 @@ Returns opacity
     void setOpacity( double opacity );
 %Docstring
 Sets opacity
+%End
+
+    DataInterpolationMethod dataInterpolationMethod() const;
+%Docstring
+Returns the type of interpolation to use to
+convert face defined datasets to
+values on vertices
+
+.. versionadded:: 3.12
+%End
+
+    void setDataInterpolationMethod( const DataInterpolationMethod &dataInterpolationMethod );
+%Docstring
+Sets data interpolation method
+
+.. versionadded:: 3.12
 %End
 
     QDomElement writeXml( QDomDocument &doc ) const;

src/analysis/CMakeLists.txt

@@ -155,6 +155,7 @@ SET(QGIS_ANALYSIS_SRCS
   mesh/qgsmeshcalcnode.cpp
   mesh/qgsmeshcalculator.cpp
   mesh/qgsmeshcalcutils.cpp
+  mesh/qgsmeshcontours.cpp
 
   network/qgsgraph.cpp
   network/qgsgraphbuilder.cpp
@@ -286,6 +287,7 @@ SET(QGIS_ANALYSIS_HDRS
   mesh/qgsmeshcalcnode.h
   mesh/qgsmeshcalculator.h
   mesh/qgsmeshcalcutils.h
+  mesh/qgsmeshcontours.h
 
   vector/qgsgeometrysnapper.h
   vector/qgsgeometrysnappersinglesource.h

src/analysis/mesh/qgsmeshcontours.cpp

@@ -0,0 +1,391 @@
+/***************************************************************************
+                          qgsmeshcontours.cpp
+                          -------------------
+    begin                : September 25th, 2019
+    copyright            : (C) 2018 by Peter Petrik
+    email                : zilolv at gmail dot com
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+#include "qgsmeshcontours.h"
+#include "qgspoint.h"
+#include "qgsmultilinestring.h"
+#include "qgslinestring.h"
+#include "qgspolygon.h"
+#include "qgsmapsettings.h"
+#include "qgsmeshlayer.h"
+#include "qgstriangularmesh.h"
+#include "qgsgeometryutils.h"
+#include "qgsfeature.h"
+#include "qgsgeometry.h"
+#include "qgsmeshlayerutils.h"
+#include "qgsmeshdataprovider.h"
+#include "qgsfeedback.h"
+
+#include <limits>
+
+#include <QSet>
+#include <QPair>
+
+QgsMeshContours::QgsMeshContours( QgsMeshLayer *layer )
+  : mMeshLayer( layer )
+{
+  if ( !mMeshLayer ||  !mMeshLayer->dataProvider() || !mMeshLayer->dataProvider()->isValid() )
+    return;
+
+  mNativeMesh.reset( new QgsMesh() );
+  mMeshLayer->dataProvider()->populateMesh( mNativeMesh.get() );
+
+  QgsMapSettings mapSettings;
+  mapSettings.setExtent( mMeshLayer->extent() );
+  mapSettings.setDestinationCrs( mMeshLayer->crs() );
+  mapSettings.setOutputDpi( 96 );
+  QgsRenderContext context = QgsRenderContext::fromMapSettings( mapSettings );
+  mTriangularMesh.reset( new QgsTriangularMesh() );
+  mTriangularMesh->update( mNativeMesh.get(), &context );
+}
+
+QgsMeshContours::~QgsMeshContours() = default;
+
+QgsGeometry QgsMeshContours::exportPolygons(
+  const QgsMeshDatasetIndex &index,
+  double min_value,
+  double max_value,
+  QgsMeshRendererScalarSettings::DataInterpolationMethod method,
+  QgsFeedback *feedback
+)
+{
+  // Check if the layer/mesh is valid
+  if ( !mTriangularMesh )
+    return QgsGeometry();
+
+  if ( min_value > max_value )
+  {
+    double tmp = max_value;
+    max_value = min_value;
+    min_value = tmp;
+  }
+
+  QVector<QgsGeometry> multiPolygon;
+
+  // STEP 1: Get Data
+  populateCache( index, method );
+  const QVector<QgsMeshVertex> vertices = mTriangularMesh->vertices();
+  const QVector<int> &trianglesToNativeFaces = mTriangularMesh->trianglesToNativeFaces();
+
+  // STEP 2: For each triangle get the contour line
+  for ( int i = 0; i < mTriangularMesh->triangles().size(); ++i )
+  {
+    if ( feedback && feedback->isCanceled() )
+      break;
+
+    int nativeIndex = trianglesToNativeFaces.at( i );
+    if ( !mScalarActiveFaceFlagValues.active( nativeIndex ) )
+      continue;
+
+    const auto &triangle = mTriangularMesh->triangles().at( i );
+    const int indices[3] =
+    {
+      triangle.at( 0 ),
+      triangle.at( 1 ),
+      triangle.at( 2 )
+    };
+
+    const QVector<QgsMeshVertex> coords =
+    {
+      vertices.at( indices[0] ),
+      vertices.at( indices[1] ),
+      vertices.at( indices[2] )
+    };
+
+    const double values[3] =
+    {
+      mDatasetValues.at( indices[0] ),
+      mDatasetValues.at( indices[1] ),
+      mDatasetValues.at( indices[2] )
+    };
+
+    // any value is NaN
+    if ( std::isnan( values[0] ) || std::isnan( values[1] ) || std::isnan( values[2] ) )
+      continue;
+
+    // all values on vertices are outside the range
+    if ( ( ( min_value > values[0] ) && ( min_value > values[1] ) && ( min_value > values[2] ) )  ||
+         ( ( max_value < values[0] ) && ( max_value < values[1] ) && ( max_value < values[2] ) ) )
+      continue;
+
+    const bool valueInRange[3] =
+    {
+      ( min_value <= values[0] ) &&( max_value >= values[0] ),
+      ( min_value <= values[1] ) &&( max_value >= values[1] ),
+      ( min_value <= values[2] ) &&( max_value >= values[2] )
+    };
+
+    // all values are inside the range == take whole triangle
+    if ( valueInRange[0] && valueInRange[1] && valueInRange[2] )
+    {
+      QVector<QgsMeshVertex> ring = coords;
+      ring.push_back( coords[0] );
+      std::unique_ptr< QgsLineString > ext = qgis::make_unique< QgsLineString> ( coords );
+      std::unique_ptr< QgsPolygon > poly = qgis::make_unique< QgsPolygon >();
+      poly->setExteriorRing( ext.release() );
+      multiPolygon.push_back( QgsGeometry( std::move( poly ) ) );
+      continue;
+    }
+
+    // go through all edges
+    QVector<QgsMeshVertex> ring;
+    for ( int i = 0; i < 3; ++i )
+    {
+      const int j = ( i + 1 ) % 3;
+
+      if ( valueInRange[i] )
+      {
+        if ( valueInRange[j] )
+        {
+          // whole edge is part of resulting contour polygon edge
+          if ( !ring.contains( coords[i] ) )
+            ring.push_back( coords[i] );
+          if ( !ring.contains( coords[j] ) )
+            ring.push_back( coords[j] );
+        }
+        else
+        {
+          // i is part or the resulting edge
+          if ( !ring.contains( coords[i] ) )
+            ring.push_back( coords[i] );
+          // we need to find the other point
+          double value = max_value;
+          if ( values[i] > values[j] )
+          {
+            value = min_value;
+          }
+          const double fraction = ( value - values[i] ) / ( values[j] - values[i] );
+          const QgsPoint xy = QgsGeometryUtils::interpolatePointOnLine( coords[i], coords[j], fraction );
+          if ( !ring.contains( xy ) )
+            ring.push_back( xy );
+        }
+      }
+      else
+      {
+        if ( valueInRange[j] )
+        {
+          // we need to find the other point
+          double value = max_value;
+          if ( values[i] < values[j] )
+          {
+            value = min_value;
+          }
+
+          const double fraction = ( value - values[i] ) / ( values[j] - values[i] );
+          const QgsPoint xy = QgsGeometryUtils::interpolatePointOnLine( coords[i], coords[j], fraction );
+          if ( !ring.contains( xy ) )
+            ring.push_back( xy );
+
+          // j is part
+          if ( !ring.contains( coords[j] ) )
+            ring.push_back( coords[j] );
+
+        }
+        else
+        {
+          // nothing here, we are outside the range
+          continue;
+        }
+      }
+    }
+
+    // add if the polygon is not degraded
+    if ( ring.size() > 2 )
+    {
+      ring.push_back( coords[0] );
+      std::unique_ptr< QgsLineString > ext = qgis::make_unique< QgsLineString> ( coords );
+      std::unique_ptr< QgsPolygon > poly = qgis::make_unique< QgsPolygon >();
+      poly->setExteriorRing( ext.release() );
+      multiPolygon.push_back( QgsGeometry( std::move( poly ) ) );
+    }
+  }
+
+  // STEP 3: dissolve the individual polygons from triangles if possible
+  if ( multiPolygon.isEmpty() )
+  {
+    return QgsGeometry();
+  }
+  else
+  {
+    QgsGeometry res = QgsGeometry::unaryUnion( multiPolygon );
+    return res;
+  }
+}
+
+QgsGeometry QgsMeshContours::exportLines( const QgsMeshDatasetIndex &index,
+    double value,
+    QgsMeshRendererScalarSettings::DataInterpolationMethod method,
+    QgsFeedback *feedback )
+{
+  // Check if the layer/mesh is valid
+  if ( !mTriangularMesh )
+    return QgsGeometry();
+
+  std::unique_ptr<QgsMultiLineString> multiLineString( new QgsMultiLineString() );
+  QSet<QPair<int, int>> exactEdges;
+
+  // STEP 1: Get Data
+  populateCache( index, method );
+  QVector<QgsMeshVertex> vertices = mTriangularMesh->vertices();
+  const QVector<int> &trianglesToNativeFaces = mTriangularMesh->trianglesToNativeFaces();
+
+  // STEP 2: For each triangle get the contour line
+  for ( int i = 0; i < mTriangularMesh->triangles().size(); ++i )
+  {
+    if ( feedback && feedback->isCanceled() )
+      break;
+
+    int nativeIndex = trianglesToNativeFaces.at( i );
+    if ( !mScalarActiveFaceFlagValues.active( nativeIndex ) )
+      continue;
+
+    const auto &triangle = mTriangularMesh->triangles().at( i );
+
+    const int indices[3] =
+    {
+      triangle.at( 0 ),
+      triangle.at( 1 ),
+      triangle.at( 2 )
+    };
+
+    const QVector<QgsMeshVertex> coords =
+    {
+      vertices.at( indices[0] ),
+      vertices.at( indices[1] ),
+      vertices.at( indices[2] )
+    };
+
+    const double values[3] =
+    {
+      mDatasetValues.at( indices[0] ),
+      mDatasetValues.at( indices[1] ),
+      mDatasetValues.at( indices[2] )
+    };
+
+    // any value is NaN
+    if ( std::isnan( values[0] ) || std::isnan( values[1] ) || std::isnan( values[2] ) )
+      continue;
+
+    // value is outside the range
+    if ( ( ( value > values[0] ) && ( value > values[1] ) && ( value > values[2] ) )  ||
+         ( ( value > values[0] ) && ( value > values[1] ) && ( value > values[2] ) ) )
+      continue;
+
+    // all values are the same
+    if ( qgsDoubleNear( values[0], values[1] ) && qgsDoubleNear( values[1], values[2] ) )
+      continue;
+
+    // go through all edges
+    QgsPoint tmp;
+
+    for ( int i = 0; i < 3; ++i )
+    {
+      const int j = ( i + 1 ) % 3;
+      // value is outside the range
+      if ( ( ( value > values[i] ) && ( value > values[j] ) ) ||
+           ( ( value < values[i] ) && ( value < values[j] ) ) )
+        continue;
+
+      // the whole edge is result and we are done
+      if ( qgsDoubleNear( values[i], values[j] ) && qgsDoubleNear( values[i], values[j] ) )
+      {
+        if ( exactEdges.contains( { indices[i], indices[j] } ) || exactEdges.contains( { indices[j], indices[i] } ) )
+        {
+          break;
+        }
+        else
+        {
+          exactEdges.insert( { indices[i], indices[j] } );
+          std::unique_ptr<QgsLineString> line( new QgsLineString( coords[i], coords[j] ) );
+          multiLineString->addGeometry( line.release() );
+          break;
+        }
+      }
+
+      // only one point matches, we are not interested in this
+      if ( qgsDoubleNear( values[i], value ) || qgsDoubleNear( values[j], value ) )
+      {
+        continue;
+      }
+
+      // ok part of the result contour line is one point on this edge
+      const double fraction = ( value - values[i] ) / ( values[j] - values[i] );
+      const QgsPoint xy = QgsGeometryUtils::interpolatePointOnLine( coords[i], coords[j], fraction );
+
+      if ( std::isnan( tmp.x() ) )
+      {
+        // ok we have found start point of the contour line
+        tmp = xy;
+      }
+      else
+      {
+        // we have found the end point of the contour line, we are done
+        std::unique_ptr<QgsLineString> line( new QgsLineString( tmp, xy ) );
+        multiLineString->addGeometry( line.release() );
+        break;
+      }
+    }
+  }
+
+  // STEP 3: merge the contour segments to linestrings
+  if ( multiLineString->isEmpty() )
+  {
+    return QgsGeometry();
+  }
+  else
+  {
+    const QgsGeometry in( multiLineString.release() );
+    const QgsGeometry res = in.mergeLines();
+    return res;
+  }
+}
+
+void QgsMeshContours::populateCache( const QgsMeshDatasetIndex &index, QgsMeshRendererScalarSettings::DataInterpolationMethod method )
+{
+  if ( mCachedIndex != index )
+  {
+    QVector<QgsMeshVertex> vertices = mTriangularMesh->vertices();
+    bool scalarDataOnVertices = mMeshLayer->dataProvider()->datasetGroupMetadata( index ).dataType() != QgsMeshDatasetGroupMetadata::DataOnFaces;
+    int count =  scalarDataOnVertices ? mNativeMesh->vertices.count() : mNativeMesh->faces.count();
+
+    // populate scalar values
+    QgsMeshDataBlock vals = mMeshLayer->dataProvider()->datasetValues(
+                              index,
+                              0,
+                              count );
+
+    // vals could be scalar or vectors, for contour rendering we want always magnitude
+    mDatasetValues = QgsMeshLayerUtils::calculateMagnitudes( vals );
+
+    // populate face active flag, always defined on faces
+    mScalarActiveFaceFlagValues = mMeshLayer->dataProvider()->areFacesActive(
+                                    index,
+                                    0,
+                                    mNativeMesh->faces.count() );
+
+    // for data on faces, there could be request to interpolate the data to vertices
+    if ( ( !scalarDataOnVertices ) )
+    {
+      mDatasetValues = QgsMeshLayerUtils::interpolateFromFacesData(
+                         mDatasetValues,
+                         mNativeMesh.get(),
+                         mTriangularMesh.get(),
+                         &mScalarActiveFaceFlagValues,
+                         method
+                       );
+    }
+  }
+}

src/analysis/mesh/qgsmeshcontours.h

@@ -0,0 +1,105 @@
+/***************************************************************************
+                          qgsmeshcontours.h
+                          -----------------
+    begin                : September 25th, 2019
+    copyright            : (C) 2018 by Peter Petrik
+    email                : zilolv at gmail dot com
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#ifndef QGSMESHCONTOURS_H
+#define QGSMESHCONTOURS_H
+
+#include <QMap>
+#include <QString>
+#include <QStringList>
+#include <QVector>
+#include <QPair>
+#include <QList>
+#include <memory>
+#include <limits>
+
+#include "qgsmeshdataprovider.h"
+#include "qgsfeature.h"
+#include "qgis_analysis.h"
+#include "qgis_sip.h"
+#include "qgstriangularmesh.h"
+#include "qgsmeshrenderersettings.h"
+
+class QgsMeshLayer;
+class QgsPoint;
+class QgsGeometry;
+class QgsFeedback;
+
+/**
+ * \ingroup analysis
+ * \class QgsMeshContours
+ *
+ * Exporter of contours lines or polygons from a mesh layer.
+ *
+ * \since QGIS 3.12
+ */
+class ANALYSIS_EXPORT QgsMeshContours
+{
+  public:
+
+    /**
+     * Constructs the mesh contours exporter.
+     * Caches the native and triangular mesh from data provider
+     * \param layer mesh layer to be associated with this exporter
+     */
+    QgsMeshContours( QgsMeshLayer *layer );
+    ~QgsMeshContours();
+
+    /**
+     * Exports multi line string containing the contour line for particular dataset and value
+     * \param index dataset index
+     * \param value value of the contour line
+     * \param method for datasets defined on faces, the method will be used to convert data to vertices
+     * \param feedback optional feedback object for progress and cancellation support
+     * \returns MultiLineString geometry containing contour lines
+     */
+    QgsGeometry exportLines( const QgsMeshDatasetIndex &index,
+                             double value,
+                             QgsMeshRendererScalarSettings::DataInterpolationMethod method,
+                             QgsFeedback *feedback = nullptr );
+
+    /**
+     * Exports multi polygons representing the areas with values in range for particular dataset
+     * \param index dataset index
+     * \param min_value minimum of the value interval for contour polygon
+     * \param max_value maximum of the value interval for contour polygon
+     * \param method for datasets defined on faces, the method will be used to convert data to vertices
+     * \param feedback optional feedback object for progress and cancellation support
+     * \returns MultiPolygon geometry containing contour polygons
+     */
+    QgsGeometry exportPolygons( const QgsMeshDatasetIndex &index,
+                                double min_value,
+                                double max_value,
+                                QgsMeshRendererScalarSettings::DataInterpolationMethod method,
+                                QgsFeedback *feedback = nullptr );
+
+  private:
+    void populateCache(
+      const QgsMeshDatasetIndex &index,
+      QgsMeshRendererScalarSettings::DataInterpolationMethod method );
+
+    QgsMeshLayer *mMeshLayer = nullptr;
+
+    // cached values for particular index & interpolation method & layer
+    QgsMeshDatasetIndex mCachedIndex;
+    std::shared_ptr<QgsTriangularMesh> mTriangularMesh;
+    std::shared_ptr<QgsMesh> mNativeMesh;
+    QgsMeshDataBlock mScalarActiveFaceFlagValues;
+    QVector<double> mDatasetValues;
+};
+
+#endif // QGSMESHCONTOURS_H

src/app/mesh/qgsmeshrendererscalarsettingswidget.cpp

@@ -27,6 +27,12 @@ QgsMeshRendererScalarSettingsWidget::QgsMeshRendererScalarSettingsWidget( QWidge
 {
   setupUi( this );
 
+  // add items to data interpolation combo box
+  mScalarInterpolationTypeComboBox->addItem( tr( "None" ), QgsMeshRendererScalarSettings::None );
+  mScalarInterpolationTypeComboBox->addItem( tr( "Neighbour Average" ), QgsMeshRendererScalarSettings::NeighbourAverage );
+  mScalarInterpolationTypeComboBox->setCurrentIndex( 0 );
+
+  // connect
   connect( mScalarRecalculateMinMaxButton, &QPushButton::clicked, this, &QgsMeshRendererScalarSettingsWidget::recalculateMinMaxButtonClicked );
   connect( mScalarMinLineEdit, &QLineEdit::textChanged, this, &QgsMeshRendererScalarSettingsWidget::minMaxChanged );
   connect( mScalarMaxLineEdit, &QLineEdit::textChanged, this, &QgsMeshRendererScalarSettingsWidget::minMaxChanged );
@@ -34,11 +40,19 @@ QgsMeshRendererScalarSettingsWidget::QgsMeshRendererScalarSettingsWidget( QWidge
   connect( mScalarMaxLineEdit, &QLineEdit::textEdited, this, &QgsMeshRendererScalarSettingsWidget::minMaxEdited );
   connect( mScalarColorRampShaderWidget, &QgsColorRampShaderWidget::widgetChanged, this, &QgsMeshRendererScalarSettingsWidget::widgetChanged );
   connect( mOpacityWidget, &QgsOpacityWidget::opacityChanged, this, &QgsMeshRendererScalarSettingsWidget::widgetChanged );
+  connect( mScalarInterpolationTypeComboBox, qgis::overload<int>::of( &QComboBox::currentIndexChanged ), this, &QgsMeshRendererScalarSettingsWidget::widgetChanged );
 }
 
 void QgsMeshRendererScalarSettingsWidget::setLayer( QgsMeshLayer *layer )
 {
   mMeshLayer = layer;
+  mScalarInterpolationTypeComboBox->setEnabled( dataIsDefinedOnFaces() );
+}
+
+void QgsMeshRendererScalarSettingsWidget::setActiveDatasetGroup( int groupIndex )
+{
+  mActiveDatasetGroup = groupIndex;
+  mScalarInterpolationTypeComboBox->setEnabled( dataIsDefinedOnFaces() );
 }
 
 QgsMeshRendererScalarSettings QgsMeshRendererScalarSettingsWidget::settings() const
@@ -47,6 +61,7 @@ QgsMeshRendererScalarSettings QgsMeshRendererScalarSettingsWidget::settings() co
   settings.setColorRampShader( mScalarColorRampShaderWidget->shader() );
   settings.setClassificationMinimumMaximum( lineEditValue( mScalarMinLineEdit ), lineEditValue( mScalarMaxLineEdit ) );
   settings.setOpacity( mOpacityWidget->opacity() );
+  settings.setDataInterpolationMethod( dataIntepolationMethod() );
   return settings;
 }
 
@@ -68,6 +83,8 @@ void QgsMeshRendererScalarSettingsWidget::syncToLayer( )
   whileBlocking( mScalarColorRampShaderWidget )->setFromShader( shader );
   whileBlocking( mScalarColorRampShaderWidget )->setMinimumMaximum( min, max );
   whileBlocking( mOpacityWidget )->setOpacity( settings.opacity() );
+  int index = mScalarInterpolationTypeComboBox->findData( settings.dataInterpolationMethod() );
+  whileBlocking( mScalarInterpolationTypeComboBox )->setCurrentIndex( index );
 }
 
 double QgsMeshRendererScalarSettingsWidget::lineEditValue( const QLineEdit *lineEdit ) const
@@ -103,3 +120,32 @@ void QgsMeshRendererScalarSettingsWidget::recalculateMinMaxButtonClicked()
   whileBlocking( mScalarMaxLineEdit )->setText( QString::number( max ) );
   mScalarColorRampShaderWidget->setMinimumMaximumAndClassify( min, max );
 }
+
+QgsMeshRendererScalarSettings::DataInterpolationMethod QgsMeshRendererScalarSettingsWidget::dataIntepolationMethod() const
+{
+  if ( dataIsDefinedOnFaces() )
+  {
+    const int data = mScalarInterpolationTypeComboBox->currentData().toInt();
+    const QgsMeshRendererScalarSettings::DataInterpolationMethod method = static_cast<QgsMeshRendererScalarSettings::DataInterpolationMethod>( data );
+    return method;
+  }
+  else
+  {
+    return QgsMeshRendererScalarSettings::None;
+  }
+}
+
+bool QgsMeshRendererScalarSettingsWidget::dataIsDefinedOnFaces() const
+{
+  if ( !mMeshLayer || !mMeshLayer->dataProvider() || !mMeshLayer->dataProvider()->isValid() )
+    return false;
+
+  if ( mActiveDatasetGroup < 0 )
+    return false;
+
+  QgsMeshDatasetGroupMetadata meta = mMeshLayer->dataProvider()->datasetGroupMetadata( mActiveDatasetGroup );
+  const bool onFaces = ( meta.dataType() == QgsMeshDatasetGroupMetadata::DataOnFaces );
+  return onFaces;
+}
+
+

src/app/mesh/qgsmeshrendererscalarsettingswidget.h

@@ -46,7 +46,7 @@ class APP_EXPORT QgsMeshRendererScalarSettingsWidget : public QWidget, private U
     void setLayer( QgsMeshLayer *layer );
 
     //! Associates a dataset group with the widget (should be set before syncToLayer())
-    void setActiveDatasetGroup( int groupIndex ) { mActiveDatasetGroup = groupIndex; }
+    void setActiveDatasetGroup( int groupIndex );
 
     //! Returns scalar settings
     QgsMeshRendererScalarSettings settings() const;
@@ -65,6 +65,9 @@ class APP_EXPORT QgsMeshRendererScalarSettingsWidget : public QWidget, private U
 
   private:
     double lineEditValue( const QLineEdit *lineEdit ) const;
+    QgsMeshRendererScalarSettings::DataInterpolationMethod dataIntepolationMethod() const;
+
+    bool dataIsDefinedOnFaces() const;
 
     QgsMeshLayer *mMeshLayer = nullptr; // not owned
     int mActiveDatasetGroup = -1;

src/core/mesh/qgsmeshlayerrenderer.cpp

@@ -80,9 +80,11 @@ void QgsMeshLayerRenderer::copyScalarDatasetValues( QgsMeshLayer *layer )
 
   // Find out if we can use cache up to date. If yes, use it and return
   const int datasetGroupCount = layer->dataProvider()->datasetGroupCount();
+  const QgsMeshRendererScalarSettings::DataInterpolationMethod method = mRendererSettings.scalarSettings( datasetIndex.group() ).dataInterpolationMethod();
   QgsMeshLayerRendererCache *cache = layer->rendererCache();
   if ( ( cache->mDatasetGroupsCount == datasetGroupCount ) &&
-       ( cache->mActiveScalarDatasetIndex == datasetIndex ) )
+       ( cache->mActiveScalarDatasetIndex == datasetIndex ) &&
+       ( cache->mDataInterpolationMethod ==  method ) )
   {
     mScalarDatasetValues = cache->mScalarDatasetValues;
     mScalarActiveFaceFlagValues = cache->mScalarActiveFaceFlagValues;
@@ -113,6 +115,19 @@ void QgsMeshLayerRenderer::copyScalarDatasetValues( QgsMeshLayer *layer )
                                     0,
                                     mNativeMesh.faces.count() );
 
+    // for data on faces, there could be request to interpolate the data to vertices
+    if ( ( !mScalarDataOnVertices ) && ( method != QgsMeshRendererScalarSettings::None ) )
+    {
+      mScalarDataOnVertices = true;
+      mScalarDatasetValues = QgsMeshLayerUtils::interpolateFromFacesData(
+                               mScalarDatasetValues,
+                               &mNativeMesh,
+                               &mTriangularMesh,
+                               &mScalarActiveFaceFlagValues,
+                               method
+                             );
+    }
+
     const QgsMeshDatasetMetadata datasetMetadata = layer->dataProvider()->datasetMetadata( datasetIndex );
     mScalarDatasetMinimum = datasetMetadata.minimum();
     mScalarDatasetMaximum = datasetMetadata.maximum();
@@ -121,6 +136,7 @@ void QgsMeshLayerRenderer::copyScalarDatasetValues( QgsMeshLayer *layer )
   // update cache
   cache->mDatasetGroupsCount = datasetGroupCount;
   cache->mActiveScalarDatasetIndex = datasetIndex;
+  cache->mDataInterpolationMethod = method;
   cache->mScalarDatasetValues = mScalarDatasetValues;
   cache->mScalarActiveFaceFlagValues = mScalarActiveFaceFlagValues;
   cache->mScalarDataOnVertices = mScalarDataOnVertices;

src/core/mesh/qgsmeshlayerrenderer.h

@@ -61,6 +61,7 @@ struct CORE_NO_EXPORT QgsMeshLayerRendererCache
   bool mScalarDataOnVertices = true;
   double mScalarDatasetMinimum = std::numeric_limits<double>::quiet_NaN();
   double mScalarDatasetMaximum = std::numeric_limits<double>::quiet_NaN();
+  QgsMeshRendererScalarSettings::DataInterpolationMethod mDataInterpolationMethod = QgsMeshRendererScalarSettings::None;
 
   // vector dataset
   QgsMeshDatasetIndex mActiveVectorDatasetIndex;

src/core/mesh/qgsmeshlayerutils.cpp

@@ -17,6 +17,7 @@
 
 #include "qgsmeshlayerutils.h"
 #include "qgsmeshtimesettings.h"
+#include "qgstriangularmesh.h"
 
 #include <limits>
 #include <QTime>
@@ -121,6 +122,59 @@ double QgsMeshLayerUtils::interpolateFromFacesData( const QgsPointXY &p1, const
   return val;
 }
 
+QVector<double> QgsMeshLayerUtils::interpolateFromFacesData( QVector<double> valuesOnFaces, QgsMesh *nativeMesh,
+    QgsTriangularMesh *triangularMesh,
+    QgsMeshDataBlock *active,
+    QgsMeshRendererScalarSettings::DataInterpolationMethod method )
+{
+  Q_UNUSED( triangularMesh )
+
+  assert( nativeMesh );
+  assert( active );
+  assert( method == QgsMeshRendererScalarSettings::NeighbourAverage );
+
+  // assuming that native vertex count = triangular vertex count
+  assert( nativeMesh->vertices.size() == triangularMesh->vertices().size() );
+  int vertexCount = triangularMesh->vertices().size();
+
+  QVector<double> res( vertexCount, 0.0 );
+  // for face datasets do simple average of the valid values of all faces that contains this vertex
+  QVector<int> count( vertexCount, 0 );
+
+  for ( int i = 0; i < nativeMesh->faces.size(); ++i )
+  {
+    if ( active->active( i ) )
+    {
+      double val = valuesOnFaces[ i ];
+      if ( !std::isnan( val ) )
+      {
+        // assign for all vertices
+        const auto &face = nativeMesh->faces.at( i );
+        for ( int j = 0; j < face.size(); ++j )
+        {
+          int vertexIndex = face[j];
+          res[vertexIndex] += val;
+          count[vertexIndex] += 1;
+        }
+      }
+    }
+  }
+
+  for ( int i = 0; i < vertexCount; ++i )
+  {
+    if ( count.at( i ) > 0 )
+    {
+      res[i] = res[i] / double( count.at( i ) );
+    }
+    else
+    {
+      res[i] = std::numeric_limits<double>::quiet_NaN();
+    }
+  }
+
+  return res;
+}
+
 QgsRectangle QgsMeshLayerUtils::triangleBoundingBox( const QgsPointXY &p1, const QgsPointXY &p2, const QgsPointXY &p3 )
 {
   // p1

src/core/mesh/qgsmeshlayerutils.h

@@ -24,6 +24,7 @@
 #include "qgsrectangle.h"
 #include "qgsmaptopixel.h"
 #include "qgsmeshdataprovider.h"
+#include "qgsmeshrenderersettings.h"
 
 #include <QVector>
 #include <QSize>
@@ -31,6 +32,8 @@
 ///@cond PRIVATE
 
 class QgsMeshTimeSettings;
+class QgsTriangularMesh;
+class QgsMeshDataBlock;
 
 /**
  * \ingroup core
@@ -95,6 +98,19 @@ class CORE_EXPORT QgsMeshLayerUtils
       const QgsPointXY &p1, const QgsPointXY &p2, const QgsPointXY &p3,
       double val, const QgsPointXY &pt );
 
+    /**
+    * Interpolate values on vertices from values on faces
+    *
+    * \since QGIS 3.12
+    */
+    static QVector<double> interpolateFromFacesData(
+      QVector<double> valuesOnFaces,
+      QgsMesh *nativeMesh,
+      QgsTriangularMesh *triangularMesh,
+      QgsMeshDataBlock *active,
+      QgsMeshRendererScalarSettings::DataInterpolationMethod method
+    );
+
     /**
      * Calculates the bounding box of the triangle
      * \param p1 first vertex of the triangle

src/core/mesh/qgsmeshrenderersettings.cpp

@@ -91,12 +91,34 @@ double QgsMeshRendererScalarSettings::opacity() const { return mOpacity; }
 
 void QgsMeshRendererScalarSettings::setOpacity( double opacity ) { mOpacity = opacity; }
 
+QgsMeshRendererScalarSettings::DataInterpolationMethod QgsMeshRendererScalarSettings::dataInterpolationMethod() const
+{
+  return mDataInterpolationMethod;
+}
+
+void QgsMeshRendererScalarSettings::setDataInterpolationMethod( const QgsMeshRendererScalarSettings::DataInterpolationMethod &dataInterpolationMethod )
+{
+  mDataInterpolationMethod = dataInterpolationMethod;
+}
+
 QDomElement QgsMeshRendererScalarSettings::writeXml( QDomDocument &doc ) const
 {
   QDomElement elem = doc.createElement( QStringLiteral( "scalar-settings" ) );
   elem.setAttribute( QStringLiteral( "min-val" ), mClassificationMinimum );
   elem.setAttribute( QStringLiteral( "max-val" ), mClassificationMaximum );
   elem.setAttribute( QStringLiteral( "opacity" ), mOpacity );
+
+  QString methodTxt;
+  switch ( mDataInterpolationMethod )
+  {
+    case None:
+      methodTxt = QStringLiteral( "none" );
+      break;
+    case NeighbourAverage:
+      methodTxt = QStringLiteral( "neighbour-average" );
+      break;
+  }
+  elem.setAttribute( QStringLiteral( "interpolation-method" ), methodTxt );
   QDomElement elemShader = mColorRampShader.writeXml( doc );
   elem.appendChild( elemShader );
   return elem;
@@ -107,6 +129,15 @@ void QgsMeshRendererScalarSettings::readXml( const QDomElement &elem )
   mClassificationMinimum = elem.attribute( QStringLiteral( "min-val" ) ).toDouble();
   mClassificationMaximum = elem.attribute( QStringLiteral( "max-val" ) ).toDouble();
   mOpacity = elem.attribute( QStringLiteral( "opacity" ) ).toDouble();
+  QString methodTxt = elem.attribute( QStringLiteral( "interpolation-method" ) );
+  if ( QStringLiteral( "neighbour-average" ) == methodTxt )
+  {
+    mDataInterpolationMethod = DataInterpolationMethod::NeighbourAverage;
+  }
+  else
+  {
+    mDataInterpolationMethod = DataInterpolationMethod::None;
+  }
   QDomElement elemShader = elem.firstChildElement( QStringLiteral( "colorrampshader" ) );
   mColorRampShader.readXml( elemShader );
 }

src/core/mesh/qgsmeshrenderersettings.h

@@ -77,6 +77,22 @@ class CORE_EXPORT QgsMeshRendererMeshSettings
 class CORE_EXPORT QgsMeshRendererScalarSettings
 {
   public:
+    //! Interpolation of value defined on vertices from datasets with data defined on faces
+    enum DataInterpolationMethod
+    {
+
+      /**
+       * Use data defined on face centers, do not interpolate to vertices
+       */
+      None = 0,
+
+      /**
+       * For each vertex does a simple average of values defined for all faces that contains
+       * given vertex
+       */
+      NeighbourAverage,
+    };
+
     //! Returns color ramp shader function
     QgsColorRampShader colorRampShader() const;
     //! Sets color ramp shader function
@@ -94,6 +110,22 @@ class CORE_EXPORT QgsMeshRendererScalarSettings
     //! Sets opacity
     void setOpacity( double opacity );
 
+    /**
+     * Returns the type of interpolation to use to
+     * convert face defined datasets to
+     * values on vertices
+     *
+     * \since QGIS 3.12
+     */
+    DataInterpolationMethod dataInterpolationMethod() const;
+
+    /**
+     * Sets data interpolation method
+     *
+     * \since QGIS 3.12
+     */
+    void setDataInterpolationMethod( const DataInterpolationMethod &dataInterpolationMethod );
+
     //! Writes configuration to a new DOM element
     QDomElement writeXml( QDomDocument &doc ) const;
     //! Reads configuration from the given DOM element
@@ -101,9 +133,11 @@ class CORE_EXPORT QgsMeshRendererScalarSettings
 
   private:
     QgsColorRampShader mColorRampShader;
+    DataInterpolationMethod mDataInterpolationMethod = DataInterpolationMethod::None;
     double mClassificationMinimum = 0;
     double mClassificationMaximum = 0;
     double mOpacity = 1;
+
 };
 
 /**

src/ui/mesh/qgsmeshrendererscalarsettingswidgetbase.ui

@@ -7,7 +7,7 @@
     <x>0</x>
     <y>0</y>
     <width>344</width>
-    <height>267</height>
+    <height>377</height>
    </rect>
   </property>
   <property name="windowTitle">
@@ -71,6 +71,20 @@
      </item>
     </layout>
    </item>
+   <item>
+    <layout class="QHBoxLayout" name="horizontalLayout_3">
+     <item>
+      <widget class="QLabel" name="mLabelDataInterpolation">
+       <property name="text">
+        <string>Resampling method</string>
+       </property>
+      </widget>
+     </item>
+     <item>
+      <widget class="QComboBox" name="mScalarInterpolationTypeComboBox"/>
+     </item>
+    </layout>
+   </item>
    <item>
     <widget class="QgsColorRampShaderWidget" name="mScalarColorRampShaderWidget" native="true">
      <property name="sizePolicy">

tests/src/analysis/CMakeLists.txt

@@ -88,6 +88,7 @@ SET(TESTS
  testqgsnetworkanalysis.cpp
  testqgsninecellfilters.cpp 
  testqgsmeshcalculator.cpp
+ testqgsmeshcontours.cpp
 )
 
 FOREACH(TESTSRC ${TESTS})

tests/src/analysis/testqgsmeshcontours.cpp

@@ -0,0 +1,225 @@
+/***************************************************************************
+  testqgsmeshcontours.cpp
+  --------------------------------------
+Date                 : September 2019
+Copyright            : (C) 2019 by Peter Petrik
+Email                : zilolv at gmail dot com
+ ***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+#include "qgstest.h"
+#include <limits>
+#include <cmath>
+#include <qdebug.h>
+#include <iostream>
+
+#include "qgsmeshcontours.h"
+#include "qgsmeshdataprovider.h"
+#include "qgsmeshlayer.h"
+#include "qgsapplication.h"
+#include "qgsproject.h"
+#include "qgsmeshmemorydataprovider.h"
+#include "qgslinesegment.h"
+#include "qgsmultilinestring.h"
+#include "qgslinestring.h"
+#include "qgsgeometryfactory.h"
+
+class TestQgsMeshContours : public QObject
+{
+    Q_OBJECT
+
+  public:
+    TestQgsMeshContours() = default;
+
+  private slots:
+    void initTestCase();// will be called before the first testfunction is executed.
+    void cleanupTestCase();// will be called after the last testfunction was executed.
+    void init() ;// will be called before each testfunction is executed.
+    void cleanup() ;// will be called after every testfunction
+
+    void equals( QgsGeometry geom, QgsGeometry expected );
+
+    void testQuadAndTriangleVertexScalarLine_data();
+    void testQuadAndTriangleVertexScalarLine();
+
+    void testQuadAndTriangleFaceScalarLine_data();
+    void testQuadAndTriangleFaceScalarLine();
+
+    void testQuadAndTriangleVertexScalarPoly_data();
+    void testQuadAndTriangleVertexScalarPoly();
+
+    void testQuadAndTriangleFaceScalarPoly_data();
+    void testQuadAndTriangleFaceScalarPoly();
+
+  private:
+    QgsMeshLayer *mpMeshLayer = nullptr;
+};
+
+void  TestQgsMeshContours::initTestCase()
+{
+  //
+  // Runs once before any tests are run
+  //
+  // init QGIS's paths - true means that all path will be inited from prefix
+  QgsApplication::init();
+  QgsApplication::initQgis();
+
+  QString testDataDir = QStringLiteral( TEST_DATA_DIR ) + QStringLiteral( "/mesh/" );
+  QString uri( testDataDir + "quad_and_triangle.2dm" );
+  mpMeshLayer = new QgsMeshLayer( uri, "Triangle and Quad MDAL", "mdal" );
+  mpMeshLayer->dataProvider()->addDataset( testDataDir + "/quad_and_triangle_vertex_scalar.dat" );
+  mpMeshLayer->dataProvider()->addDataset( testDataDir + "/quad_and_triangle_els_face_scalar.dat" );
+  QVERIFY( mpMeshLayer->isValid() );
+}
+
+void  TestQgsMeshContours::cleanupTestCase()
+{
+  QgsApplication::exitQgis();
+}
+
+void  TestQgsMeshContours::init()
+{
+}
+void  TestQgsMeshContours::cleanup()
+{
+}
+
+void TestQgsMeshContours::equals( QgsGeometry geom, QgsGeometry expected )
+{
+  if ( expected.isNull() )
+  {
+    QVERIFY( geom.isNull() );
+    return;
+  }
+
+  QVERIFY( !geom.isNull() );
+
+  QVERIFY( geom.type() == expected.type() );
+  QVERIFY( geom.isMultipart() == expected.isMultipart() );
+
+  const QString gWkt = geom.asWkt();
+  const QString eWkt = expected.asWkt();
+
+  QCOMPARE( gWkt.trimmed(), eWkt.trimmed() );
+}
+
+void TestQgsMeshContours::testQuadAndTriangleVertexScalarLine_data()
+{
+  QTest::addColumn< double >( "value" );
+  QTest::addColumn< QgsGeometry >( "expected" );
+
+  QTest::newRow( "left_edge" ) << 1.0 << QgsGeometry( QgsGeometryFactory::geomFromWkt( "LineStringZ (1000 3000 10, 1000 2000 20)" ) );
+  QTest::newRow( "two_intersections" ) <<  1.5 << QgsGeometry( QgsGeometryFactory::geomFromWkt( "LineStringZ (1500 3000 30, 1500 2500 35, 1500 2000 25)" ) );
+  QTest::newRow( "between two triangles" ) <<  2.0 << QgsGeometry( QgsGeometryFactory::geomFromWkt( "LineStringZ (2000 2000 30, 2000 3000 50)" ) );
+
+  QTest::newRow( "only one point" ) <<  3.0 << QgsGeometry();
+  QTest::newRow( "outside lower" ) << -3.0 << QgsGeometry();
+  QTest::newRow( "outside upper" ) <<  4.0 << QgsGeometry();
+}
+
+void TestQgsMeshContours::testQuadAndTriangleVertexScalarLine()
+{
+  QFETCH( double, value );
+  QFETCH( QgsGeometry, expected );
+
+  QgsMeshDatasetIndex datasetIndex( 1, 0 );
+
+  QgsMeshContours contours( mpMeshLayer );
+
+  QgsGeometry res = contours.exportLines( datasetIndex, value, QgsMeshRendererScalarSettings::None );
+  equals( res, expected );
+}
+
+void TestQgsMeshContours::testQuadAndTriangleFaceScalarLine_data()
+{
+  QTest::addColumn< double >( "value" );
+  QTest::addColumn< QgsGeometry >( "expected" );
+
+  QTest::newRow( "left_edge" ) <<  1.0 << QgsGeometry( QgsGeometryFactory::geomFromWkt( "LineStringZ (1000 3000 10, 1000 2000 20)" ) );
+  QTest::newRow( "two_intersections" ) << 1.25 << QgsGeometry( QgsGeometryFactory::geomFromWkt( "LineStringZ (1500 3000 30, 1500 2500 35, 1500 2000 25)" ) );
+  QTest::newRow( "between two triangles" ) <<  1.5 << QgsGeometry( QgsGeometryFactory::geomFromWkt( "LineStringZ (2000 2000 30, 2000 3000 50)" ) );
+
+  QTest::newRow( "only one point" ) <<  2.0 << QgsGeometry();
+  QTest::newRow( "outside lower" ) <<  -3.0 << QgsGeometry();
+  QTest::newRow( "outside upper" ) <<  3.0 << QgsGeometry();
+}
+
+void TestQgsMeshContours::testQuadAndTriangleFaceScalarLine()
+{
+  QFETCH( double, value );
+  QFETCH( QgsGeometry, expected );
+
+  QgsMeshDatasetIndex datasetIndex( 2, 0 );
+
+  QgsMeshContours contours( mpMeshLayer );
+
+  QgsGeometry res = contours.exportLines( datasetIndex, value, QgsMeshRendererScalarSettings::NeighbourAverage );
+  equals( res, expected );
+}
+
+void TestQgsMeshContours::testQuadAndTriangleVertexScalarPoly_data()
+{
+  QTest::addColumn< double >( "min_value" );
+  QTest::addColumn< double >( "max_value" );
+  QTest::addColumn< QgsGeometry >( "expected" );
+
+  QTest::newRow( "left" ) <<  1.0 << 1.5 << QgsGeometry( QgsGeometryFactory::geomFromWkt( "PolygonZ ((1000 2000 20, 1000 3000 10, 2000 3000 50, 2000 2000 30, 1000 2000 20))" ) );
+  QTest::newRow( "middle" ) <<  1.5 << 2.0 << QgsGeometry( QgsGeometryFactory::geomFromWkt( "PolygonZ ((1000 2000 20, 1000 3000 10, 2000 3000 50, 2000 2000 30, 1000 2000 20))" ) );
+  QTest::newRow( "right" ) <<  2.0 << 2.5 << QgsGeometry( QgsGeometryFactory::geomFromWkt( "PolygonZ ((3000 2000 40, 2000 3000 50, 2000 2000 30, 3000 2000 40))" ) );
+
+  QTest::newRow( "outside_left_edge" ) <<  0.0 << 1.0 << QgsGeometry();
+  QTest::newRow( "only one point" ) << 3.0 << 3.0 << QgsGeometry();
+  QTest::newRow( "outside lower" ) <<  -3.0 << -4.0 << QgsGeometry();
+  QTest::newRow( "outside upper" ) << 4.0 << 5.0 << QgsGeometry();
+}
+
+void TestQgsMeshContours::testQuadAndTriangleVertexScalarPoly()
+{
+  QFETCH( double, min_value );
+  QFETCH( double, max_value );
+  QFETCH( QgsGeometry, expected );
+
+  QgsMeshDatasetIndex datasetIndex( 1, 0 );
+
+  QgsMeshContours contours( mpMeshLayer );
+
+  QgsGeometry res = contours.exportPolygons( datasetIndex, min_value, max_value, QgsMeshRendererScalarSettings::None );
+  equals( res, expected );
+}
+
+void TestQgsMeshContours::testQuadAndTriangleFaceScalarPoly_data()
+{
+  QTest::addColumn< double >( "min_value" );
+  QTest::addColumn< double >( "max_value" );
+  QTest::addColumn< QgsGeometry >( "expected" );
+
+  QTest::newRow( "left" ) << 1.0 << 1.25 << QgsGeometry( QgsGeometryFactory::geomFromWkt( "PolygonZ ((1000 2000 20, 1000 3000 10, 2000 3000 50, 2000 2000 30, 1000 2000 20))" ) );
+  QTest::newRow( "middle" ) << 1.25 << 1.5 <<  QgsGeometry( QgsGeometryFactory::geomFromWkt( "PolygonZ ((1000 2000 20, 1000 3000 10, 2000 3000 50, 2000 2000 30, 1000 2000 20))" ) );
+  QTest::newRow( "right" ) <<  1.5 << 1.75 << QgsGeometry( QgsGeometryFactory::geomFromWkt( "PolygonZ ((3000 2000 40, 2000 3000 50, 2000 2000 30, 3000 2000 40))" ) );
+
+  QTest::newRow( "only one point" ) << 2.0 << 2.0 << QgsGeometry();
+  QTest::newRow( "outside lower" ) << -3.0 << -4.0 << QgsGeometry();
+  QTest::newRow( "outside upper" ) << 3.0 << 4.0 << QgsGeometry();
+}
+
+void TestQgsMeshContours::testQuadAndTriangleFaceScalarPoly()
+{
+  QFETCH( double, min_value );
+  QFETCH( double, max_value );
+  QFETCH( QgsGeometry, expected );
+
+  QgsMeshDatasetIndex datasetIndex( 2, 0 );
+
+  QgsMeshContours contours( mpMeshLayer );
+
+  QgsGeometry res = contours.exportPolygons( datasetIndex, min_value, max_value, QgsMeshRendererScalarSettings::NeighbourAverage );
+  equals( res, expected );
+}
+
+QGSTEST_MAIN( TestQgsMeshContours )
+#include "testqgsmeshcontours.moc"

tests/src/core/testqgsmeshlayerrenderer.cpp

@@ -71,6 +71,7 @@ class TestQgsMeshRenderer : public QObject
     void test_vertex_scalar_dataset_rendering();
     void test_vertex_vector_dataset_rendering();
     void test_face_scalar_dataset_rendering();
+    void test_face_scalar_dataset_interpolated_neighbour_average_rendering();
     void test_face_vector_dataset_rendering();
     void test_vertex_scalar_dataset_with_inactive_face_rendering();
     void test_face_vector_on_user_grid();
@@ -242,6 +243,23 @@ void TestQgsMeshRenderer::test_face_scalar_dataset_rendering()
   QVERIFY( imageCheck( "quad_and_triangle_face_scalar_dataset", mMemoryLayer ) );
 }
 
+void TestQgsMeshRenderer::test_face_scalar_dataset_interpolated_neighbour_average_rendering()
+{
+  QgsMeshDatasetIndex ds( 2, 0 );
+  const QgsMeshDatasetGroupMetadata metadata = mMemoryLayer->dataProvider()->datasetGroupMetadata( ds );
+  QVERIFY( metadata.name() == "FaceScalarDataset" );
+
+  QgsMeshRendererSettings rendererSettings = mMemoryLayer->rendererSettings();
+  rendererSettings.setActiveScalarDataset( ds );
+  auto scalarRendererSettings = rendererSettings.scalarSettings( 2 );
+  scalarRendererSettings.setDataInterpolationMethod( QgsMeshRendererScalarSettings::NeighbourAverage );
+  rendererSettings.setScalarSettings( 2, scalarRendererSettings );
+  mMemoryLayer->setRendererSettings( rendererSettings );
+
+  QVERIFY( imageCheck( "quad_and_triangle_face_scalar_interpolated_neighbour_average_dataset", mMemoryLayer ) );
+}
+
+
 void TestQgsMeshRenderer::test_face_vector_dataset_rendering()
 {
   QgsMeshDatasetIndex ds( 3, 0 );