Merge pull request #7582 from PeterPetrik/mesh_plots

[mesh] [feature] function to get value for the point on map
qgis · Aug 10, 2018 · ad4ddb1 · ad4ddb1
2 parents af075bf + 271bab1
commit ad4ddb1
Show file tree

Hide file tree

Showing 14 changed files with 229 additions and 35 deletions.
diff --git a/python/core/auto_generated/mesh/qgsmeshdataprovider.sip.in b/python/core/auto_generated/mesh/qgsmeshdataprovider.sip.in
@@ -140,6 +140,13 @@ such as whether the data is vector or scalar, name
 #include "qgsmeshdataprovider.h"
 %End
   public:
+
+    enum DataType
+    {
+      DataOnFaces,
+      DataOnVertices
+    };
+
     QgsMeshDatasetGroupMetadata();
 %Docstring
 Constructs an empty metadata object
@@ -178,9 +185,9 @@ Returns whether dataset group has vector data
 Returns whether dataset group has scalar data
 %End
 
-    bool isOnVertices() const;
+    DataType dataType() const;
 %Docstring
-Returns whether dataset group data is defined on vertices
+Returns whether dataset group data is defined on vertices or faces
 %End
 
 };

diff --git a/python/core/auto_generated/mesh/qgsmeshlayer.sip.in b/python/core/auto_generated/mesh/qgsmeshlayer.sip.in
@@ -189,6 +189,43 @@ If dataset is not vector based, do nothing. Triggers repaint
     QgsMeshDatasetIndex activeVectorDataset() const;
 %Docstring
 Returns active vector dataset
+%End
+
+    QgsMeshDatasetValue datasetValue( const QgsMeshDatasetIndex &index, const QgsPointXY &point ) const;
+%Docstring
+Interpolates the value on the given point from given dataset.
+
+.. note::
+
+   It uses previously cached and indexed triangular mesh
+   and so if the layer has not been rendered previously
+   (e.g. when used in a script) it returns NaN value
+
+:param index: dataset index specifying group and dataset to extract value from
+:param point: point to query in map coordinates
+
+:return: interpolated value at the point. Returns NaN values for values
+         outside the mesh layer, nodata values and in case triangular mesh was not
+         previously used for rendering
+
+
+.. versionadded:: 3.4
+%End
+
+  signals:
+
+    void activeScalarDatasetChanged( const QgsMeshDatasetIndex &index );
+%Docstring
+Emitted when active scalar dataset is changed
+
+.. versionadded:: 3.4
+%End
+
+    void activeVectorDatasetChanged( const QgsMeshDatasetIndex &index );
+%Docstring
+Emitted when active vector dataset is changed
+
+.. versionadded:: 3.4
 %End
 
   private: // Private methods

diff --git a/src/app/mesh/qgsmeshrendereractivedatasetwidget.cpp b/src/app/mesh/qgsmeshrendereractivedatasetwidget.cpp
@@ -126,8 +126,8 @@ void QgsMeshRendererActiveDatasetWidget::updateMetadata( QgsMeshDatasetIndex dat
 
     const QgsMeshDatasetGroupMetadata gmeta = mMeshLayer->dataProvider()->datasetGroupMetadata( datasetIndex );
     msg += QStringLiteral( "<tr><td>%1</td><td>%2</td></tr>" )
-           .arg( tr( "Is on vertices" ) )
-           .arg( gmeta.isOnVertices() ? tr( "Yes" ) : tr( "No" ) );
+           .arg( tr( "Data Type" ) )
+           .arg( gmeta.dataType() == QgsMeshDatasetGroupMetadata::DataOnVertices ? tr( "Defined on vertices" ) : tr( "Defined on faces" ) );
 
     msg += QStringLiteral( "<tr><td>%1</td><td>%2</td></tr>" )
            .arg( tr( "Is vector" ) )

diff --git a/src/app/mesh/qgsmeshrendererscalarsettingswidget.cpp b/src/app/mesh/qgsmeshrendererscalarsettingswidget.cpp
@@ -112,7 +112,7 @@ void QgsMeshRendererScalarSettingsWidget::calcMinMax( QgsMeshDatasetIndex datase
     return;
 
   const QgsMeshDatasetGroupMetadata metadata = mMeshLayer->dataProvider()->datasetGroupMetadata( datasetIndex );
-  bool scalarDataOnVertices = metadata.isOnVertices();
+  bool scalarDataOnVertices = metadata.dataType() == QgsMeshDatasetGroupMetadata::DataOnVertices;
   int count;
   if ( scalarDataOnVertices )
     count = mMeshLayer->dataProvider()->vertexCount();

diff --git a/src/core/mesh/qgsmeshdataprovider.cpp b/src/core/mesh/qgsmeshdataprovider.cpp
@@ -161,15 +161,16 @@ bool QgsMeshDatasetGroupMetadata::isScalar() const
   return mIsScalar;
 }
 
+
+
 QString QgsMeshDatasetGroupMetadata::name() const
 {
   return mName;
 }
 
-
-bool QgsMeshDatasetGroupMetadata::isOnVertices() const
+QgsMeshDatasetGroupMetadata::DataType QgsMeshDatasetGroupMetadata::dataType() const
 {
-  return mIsOnVertices;
+  return ( mIsOnVertices ) ? DataType::DataOnVertices : DataType::DataOnFaces;
 }
 
 int QgsMeshDatasetSourceInterface::datasetCount( QgsMeshDatasetIndex index ) const

diff --git a/src/core/mesh/qgsmeshdataprovider.h b/src/core/mesh/qgsmeshdataprovider.h
@@ -127,6 +127,14 @@ class CORE_EXPORT QgsMeshDatasetValue
 class CORE_EXPORT QgsMeshDatasetGroupMetadata
 {
   public:
+
+    //! Location of where data is specified for datasets in the dataset group
+    enum DataType
+    {
+      DataOnFaces, //!< Data is defined on faces
+      DataOnVertices //!< Data is defined on vertices
+    };
+
     //! Constructs an empty metadata object
     QgsMeshDatasetGroupMetadata() = default;
 
@@ -164,9 +172,9 @@ class CORE_EXPORT QgsMeshDatasetGroupMetadata
     bool isScalar() const;
 
     /**
-     * \brief Returns whether dataset group data is defined on vertices
+     * \brief Returns whether dataset group data is defined on vertices or faces
      */
-    bool isOnVertices() const;
+    DataType dataType() const;
 
   private:
     QString mName;

diff --git a/src/core/mesh/qgsmeshlayer.cpp b/src/core/mesh/qgsmeshlayer.cpp
@@ -16,6 +16,7 @@
  ***************************************************************************/
 
 #include <cstddef>
+#include <limits>
 
 #include <QUuid>
 
@@ -27,6 +28,7 @@
 #include "qgsproviderregistry.h"
 #include "qgsreadwritecontext.h"
 #include "qgstriangularmesh.h"
+#include "qgsmeshlayerinterpolator.h"
 
 QgsMeshLayer::QgsMeshLayer( const QString &meshLayerPath,
                             const QString &baseName,
@@ -154,6 +156,8 @@ void QgsMeshLayer::setActiveScalarDataset( QgsMeshDatasetIndex index )
     mActiveScalarDataset = QgsMeshDatasetIndex();
 
   triggerRepaint();
+
+  emit activeScalarDatasetChanged( mActiveScalarDataset );
 }
 
 void QgsMeshLayer::setActiveVectorDataset( QgsMeshDatasetIndex index )
@@ -175,6 +179,43 @@ void QgsMeshLayer::setActiveVectorDataset( QgsMeshDatasetIndex index )
   }
 
   triggerRepaint();
+
+  emit activeVectorDatasetChanged( mActiveVectorDataset );
+}
+
+QgsMeshDatasetValue QgsMeshLayer::datasetValue( const QgsMeshDatasetIndex &index, const QgsPointXY &point ) const
+{
+  QgsMeshDatasetValue value;
+
+  if ( mTriangularMesh && dataProvider() && dataProvider()->isValid() && index.isValid() )
+  {
+    int faceIndex = mTriangularMesh->faceIndexForPoint( point ) ;
+    if ( faceIndex >= 0 )
+    {
+      if ( dataProvider()->datasetGroupMetadata( index ).dataType() == QgsMeshDatasetGroupMetadata::DataOnFaces )
+      {
+        int nativeFaceIndex = mTriangularMesh->trianglesToNativeFaces().at( faceIndex );
+        return dataProvider()->datasetValue( index, nativeFaceIndex );
+      }
+      else
+      {
+        const QgsMeshFace &face = mTriangularMesh->triangles()[faceIndex];
+        const int v1 = face[0], v2 = face[1], v3 = face[2];
+        const QgsPoint p1 = mTriangularMesh->vertices()[v1], p2 = mTriangularMesh->vertices()[v2], p3 = mTriangularMesh->vertices()[v3];
+        const QgsMeshDatasetValue val1 = dataProvider()->datasetValue( index, v1 );
+        const QgsMeshDatasetValue val2 = dataProvider()->datasetValue( index, v2 );
+        const QgsMeshDatasetValue val3 = dataProvider()->datasetValue( index, v3 );
+        const double x = QgsMeshLayerInterpolator::interpolateFromVerticesData( p1, p2, p3, val1.x(), val2.x(), val3.x(), point );
+        double y = std::numeric_limits<double>::quiet_NaN();
+        bool isVector = dataProvider()->datasetGroupMetadata( index ).isVector();
+        if ( isVector )
+          y = QgsMeshLayerInterpolator::interpolateFromVerticesData( p1, p2, p3, val1.y(), val2.y(), val3.y(), point );
+
+        return QgsMeshDatasetValue( x, y );
+      }
+    }
+  }
+  return value;
 }
 
 void QgsMeshLayer::fillNativeMesh()

diff --git a/src/core/mesh/qgsmeshlayer.h b/src/core/mesh/qgsmeshlayer.h
@@ -183,6 +183,39 @@ class CORE_EXPORT QgsMeshLayer : public QgsMapLayer
     //! Returns active vector dataset
     QgsMeshDatasetIndex activeVectorDataset() const { return mActiveVectorDataset; }
 
+    /**
+      * Interpolates the value on the given point from given dataset.
+      *
+      * \note It uses previously cached and indexed triangular mesh
+      * and so if the layer has not been rendered previously
+      * (e.g. when used in a script) it returns NaN value
+      *
+      * \param index dataset index specifying group and dataset to extract value from
+      * \param point point to query in map coordinates
+      * \returns interpolated value at the point. Returns NaN values for values
+      * outside the mesh layer, nodata values and in case triangular mesh was not
+      * previously used for rendering
+      *
+      * \since QGIS 3.4
+      */
+    QgsMeshDatasetValue datasetValue( const QgsMeshDatasetIndex &index, const QgsPointXY &point ) const;
+
+  signals:
+
+    /**
+     * Emitted when active scalar dataset is changed
+     *
+     * \since QGIS 3.4
+     */
+    void activeScalarDatasetChanged( const QgsMeshDatasetIndex &index );
+
+    /**
+     * Emitted when active vector dataset is changed
+     *
+     * \since QGIS 3.4
+     */
+    void activeVectorDatasetChanged( const QgsMeshDatasetIndex &index );
+
   private: // Private methods
 
     /**

diff --git a/src/core/mesh/qgsmeshlayerinterpolator.cpp b/src/core/mesh/qgsmeshlayerinterpolator.cpp
@@ -88,9 +88,8 @@ static bool E3T_physicalToBarycentric( const QgsPointXY &pA, const QgsPointXY &p
   return true;
 }
 
-
-double interpolateFromVerticesData( const QgsPointXY &p1, const QgsPointXY &p2, const QgsPointXY &p3,
-                                    double val1, double val2, double val3, const QgsPointXY &pt )
+double QgsMeshLayerInterpolator::interpolateFromVerticesData( const QgsPointXY &p1, const QgsPointXY &p2, const QgsPointXY &p3,
+    double val1, double val2, double val3, const QgsPointXY &pt )
 {
   double lam1, lam2, lam3;
   if ( !E3T_physicalToBarycentric( p1, p2, p3, pt, lam1, lam2, lam3 ) )

diff --git a/src/core/mesh/qgsmeshlayerinterpolator.h b/src/core/mesh/qgsmeshlayerinterpolator.h
@@ -59,6 +59,22 @@ class QgsMeshLayerInterpolator : public QgsRasterInterface
     int bandCount() const override;
     QgsRasterBlock *block( int, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback = nullptr ) override;
 
+    /*
+    * Interpolates value based on known values on the vertices of a triangle
+    * \param p1 first vertex of the triangle
+    * \param p2 second vertex of the triangle
+    * \param p3 third vertex of the triangle
+    * \param val1 value on p1 of the triangle
+    * \param val2 value on p2 of the triangle
+    * \param val3 value on p3 of the triangle
+    * \param pt point where to calculate value
+    * \returns value on the point pt or NaN in case the point is outside the triangle
+    */
+    static double interpolateFromVerticesData(
+      const QgsPointXY &p1, const QgsPointXY &p2, const QgsPointXY &p3,
+      double val1, double val2, double val3, const QgsPointXY &pt
+    );
+
   private:
     const QgsTriangularMesh &mTriangularMesh;
     const QVector<double> &mDatasetValues;

diff --git a/src/core/mesh/qgsmeshlayerrenderer.cpp b/src/core/mesh/qgsmeshlayerrenderer.cpp
@@ -140,7 +140,7 @@ void QgsMeshLayerRenderer::copyScalarDatasetValues( QgsMeshLayer *layer )
   if ( datasetIndex.isValid() )
   {
     const QgsMeshDatasetGroupMetadata metadata = layer->dataProvider()->datasetGroupMetadata( datasetIndex );
-    mScalarDataOnVertices = metadata.isOnVertices();
+    mScalarDataOnVertices = metadata.dataType() == QgsMeshDatasetGroupMetadata::DataOnVertices;
     int count;
     if ( mScalarDataOnVertices )
       count = mNativeMesh.vertices.count();
@@ -172,7 +172,7 @@ void QgsMeshLayerRenderer::copyVectorDatasetValues( QgsMeshLayer *layer )
     }
     else
     {
-      mVectorDataOnVertices = metadata.isOnVertices();
+      mVectorDataOnVertices = metadata.dataType() == QgsMeshDatasetGroupMetadata::DataOnVertices;
       int count;
       if ( mVectorDataOnVertices )
         count = mNativeMesh.vertices.count();
@@ -229,17 +229,7 @@ void QgsMeshLayerRenderer::renderMesh( const std::unique_ptr<QgsSymbol> &symbol,
     const QgsMeshFace &face = faces[i];
     QgsFeature feat;
     feat.setFields( fields );
-    QVector<QgsPointXY> ring;
-    for ( int j = 0; j < face.size(); ++j )
-    {
-      int vertex_id = face[j];
-      Q_ASSERT( vertex_id < mTriangularMesh.vertices().size() ); //Triangular mesh vertices contains also native mesh vertices
-      const QgsPoint &vertex = mTriangularMesh.vertices()[vertex_id];
-      ring.append( vertex );
-    }
-    QgsPolygonXY polygon;
-    polygon.append( ring );
-    QgsGeometry geom = QgsGeometry::fromPolygonXY( polygon );
+    QgsGeometry geom = QgsMeshUtils::toGeometry( face, mTriangularMesh.vertices() ); //Triangular mesh vertices contains also native mesh vertices
     feat.setGeometry( geom );
     renderer.renderFeature( feat, mContext );
   }

diff --git a/src/core/mesh/qgstriangularmesh.cpp b/src/core/mesh/qgstriangularmesh.cpp
@@ -15,6 +15,11 @@
  *                                                                         *
  ***************************************************************************/
 
+#include <memory>
+#include <QList>
+#include "qgsfeature.h"
+#include "qgspolygon.h"
+#include "qgslinestring.h"
 #include "qgstriangularmesh.h"
 #include "qgsrendercontext.h"
 #include "qgscoordinatetransform.h"
@@ -63,12 +68,12 @@ static void ENP_centroid( const QPolygonF &pX, double &cx, double &cy )
   cy /= ( 6.0 * signedArea );
 }
 
-
 void QgsTriangularMesh::update( QgsMesh *nativeMesh, QgsRenderContext *context )
 {
   Q_ASSERT( nativeMesh );
   Q_ASSERT( context );
 
+  mSpatialIndex = QgsSpatialIndex();
   mTriangularMesh.vertices.clear();
   mTriangularMesh.faces.clear();
   mTrianglesToNativeFaces.clear();
@@ -143,6 +148,14 @@ void QgsTriangularMesh::update( QgsMesh *nativeMesh, QgsRenderContext *context )
     ENP_centroid( poly, cx, cy );
     mNativeMeshFaceCentroids[i] = QgsMeshVertex( cx, cy );
   }
+
+  // CALCULATE SPATIAL INDEX
+  for ( int i = 0; i < mTriangularMesh.faces.size(); ++i )
+  {
+    const QgsMeshFace &face = mTriangularMesh.faces.at( i ) ;
+    QgsGeometry geom = QgsMeshUtils::toGeometry( face, mTriangularMesh.vertices );
+    ( void )mSpatialIndex.insertFeature( i, geom.boundingBox() );
+  }
 }
 
 const QVector<QgsMeshVertex> &QgsTriangularMesh::vertices() const
@@ -165,3 +178,31 @@ const QVector<int> &QgsTriangularMesh::trianglesToNativeFaces() const
   return mTrianglesToNativeFaces;
 }
 
+int QgsTriangularMesh::faceIndexForPoint( const QgsPointXY &point ) const
+{
+  const QList<QgsFeatureId> faceIndexes = mSpatialIndex.intersects( QgsRectangle( point, point ) );
+  for ( const QgsFeatureId fid : faceIndexes )
+  {
+    int faceIndex = static_cast<int>( fid );
+    const QgsMeshFace &face = mTriangularMesh.faces.at( faceIndex );
+    const QgsGeometry geom = QgsMeshUtils::toGeometry( face, mTriangularMesh.vertices );
+    if ( geom.contains( &point ) )
+      return faceIndex;
+  }
+  return -1;
+}
+
+QgsGeometry QgsMeshUtils::toGeometry( const QgsMeshFace &face, const QVector<QgsMeshVertex> &vertices )
+{
+  QVector<QgsPoint> ring;
+  for ( int j = 0; j < face.size(); ++j )
+  {
+    int vertexId = face[j];
+    Q_ASSERT( vertexId < vertices.size() );
+    const QgsPoint &vertex = vertices[vertexId];
+    ring.append( vertex );
+  }
+  std::unique_ptr< QgsPolygon > polygon = qgis::make_unique< QgsPolygon >();
+  polygon->setExteriorRing( new QgsLineString( ring ) );
+  return QgsGeometry( std::move( polygon ) );
+}