identify tool for 1d mesh

qgis · Apr 28, 2020 · a6ef994 · a6ef994
1 parent f7b001f
commit a6ef994
Show file tree

Hide file tree

Showing 6 changed files with 145 additions and 7 deletions.
diff --git a/src/core/mesh/qgsmeshlayer.cpp b/src/core/mesh/qgsmeshlayer.cpp
@@ -259,18 +259,17 @@ QString QgsMeshLayer::formatTime( double hours )
     return QgsMeshLayerUtils::formatTime( hours, QDateTime(), mTimeSettings );
 }
 
-QgsMeshDatasetValue QgsMeshLayer::datasetValue( const QgsMeshDatasetIndex &index, const QgsPointXY &point ) const
+QgsMeshDatasetValue QgsMeshLayer::datasetValue( const QgsMeshDatasetIndex &index, const QgsPointXY &point, double searchRadius ) const
 {
   QgsMeshDatasetValue value;
-
   const QgsTriangularMesh *mesh = triangularMesh();
 
   if ( mesh && dataProvider() && dataProvider()->isValid() && index.isValid() )
   {
     if ( dataProvider()->contains( QgsMesh::ElementType::Edge ) )
     {
-      // Identify for edges not implemented yet
-      return value;
+      QgsRectangle searchRectangle( point.x() - searchRadius, point.y() - searchRadius, point.x() + searchRadius, point.y() + searchRadius );
+      return dataset1dValue( index, point, searchRadius );
     }
     int faceIndex = mesh->faceIndexForPoint( point ) ;
     if ( faceIndex >= 0 )
@@ -352,6 +351,68 @@ QgsMesh3dDataBlock QgsMeshLayer::dataset3dValue( const QgsMeshDatasetIndex &inde
   return block3d;
 }
 
+QgsMeshDatasetValue QgsMeshLayer::dataset1dValue( const QgsMeshDatasetIndex &index, const QgsPointXY &point, double searchRadius ) const
+{
+  QgsMeshDatasetValue value;
+  QgsRectangle searchRectangle( point.x() - searchRadius, point.y() - searchRadius, point.x() + searchRadius, point.y() + searchRadius );
+  const QgsTriangularMesh *mesh = triangularMesh();
+  if ( mesh &&
+       mesh->contains( QgsMesh::Edge ) &&
+       mDataProvider->isValid() &&
+       index.isValid() )
+  {
+    // search for the closest edge in rectangle from point
+    const QList<int> edgeIndexes = mesh->edgeIndexesForRectangle( searchRectangle );
+    int selectedIndex = -1;
+    double sqrMaxDistFromPoint = pow( searchRadius, 2 );
+    QgsPointXY projectedPoint;
+    for ( const int edgeIndex : edgeIndexes )
+    {
+      const QgsMeshEdge &edge = mesh->edges().at( edgeIndex );
+      const QgsMeshVertex &vertex1 = mesh->vertices()[edge.first];
+      const QgsMeshVertex &vertex2 = mesh->vertices()[edge.second];
+      QgsPointXY projPoint;
+      double sqrDist = point.sqrDistToSegment( vertex1.x(), vertex1.y(), vertex2.x(), vertex2.y(), projPoint );
+      if ( sqrDist < sqrMaxDistFromPoint )
+      {
+        selectedIndex = edgeIndex;
+        projectedPoint = projPoint;
+        sqrMaxDistFromPoint = sqrDist;
+      }
+    }
+
+    if ( selectedIndex >= 0 )
+    {
+      const QgsMeshDatasetGroupMetadata::DataType dataType = dataProvider()->datasetGroupMetadata( index ).dataType();
+      switch ( dataType )
+      {
+        case QgsMeshDatasetGroupMetadata::DataOnEdges:
+        {
+          value = dataProvider()->datasetValue( index, selectedIndex );
+        }
+        break;
+
+        case QgsMeshDatasetGroupMetadata::DataOnVertices:
+        {
+          const QgsMeshEdge &edge = mesh->edges()[selectedIndex];
+          const int v1 = edge.first, v2 = edge.second;
+          const QgsPoint p1 = mesh->vertices()[v1], p2 = mesh->vertices()[v2];
+          const QgsMeshDatasetValue val1 = dataProvider()->datasetValue( index, v1 );
+          const QgsMeshDatasetValue val2 = dataProvider()->datasetValue( index, v2 );
+          double edgeLength = p1.distance( p2 );
+          double dist1 = p1.distance( projectedPoint.x(), projectedPoint.y() );
+          value = QgsMeshLayerUtils::interpolateFromVerticesData( dist1 / edgeLength, val1, val2 );
+        }
+        break;
+        default:
+          break;
+      }
+    }
+  }
+
+  return value;
+}
+
 void QgsMeshLayer::setTransformContext( const QgsCoordinateTransformContext &transformContext )
 {
   if ( mDataProvider )

diff --git a/src/core/mesh/qgsmeshlayer.h b/src/core/mesh/qgsmeshlayer.h
@@ -204,6 +204,8 @@ class CORE_EXPORT QgsMeshLayer : public QgsMapLayer
      * Gets native mesh and updates (creates if it doesn't exist) the base triangular mesh
      *
      * \param transform Transformation from layer CRS to destination (e.g. map) CRS. With invalid transform, it keeps the native mesh CRS
+     *
+     * \since QGIS 3.14
      */
     void updateTriangularMesh( const QgsCoordinateTransform &transform = QgsCoordinateTransform() );
 
@@ -257,27 +259,32 @@ class CORE_EXPORT QgsMeshLayer : public QgsMapLayer
 
     /**
       * Interpolates the value on the given point from given dataset.
+      * For 3D datasets, it uses #dataset3dValue \n
+      * For 1D datasets, it uses #dataset1dValue with \a searchRadius
       *
       * \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
+      * \see updateTriangularMesh
       *
       * \param index dataset index specifying group and dataset to extract value from
       * \param point point to query in map coordinates
+      * \param searchRadius the radius of the search area in map unit
       * \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;
+    QgsMeshDatasetValue datasetValue( const QgsMeshDatasetIndex &index, const QgsPointXY &point, double searchRadius = 0 ) const;
 
     /**
       * Returns the 3d values of stacked 3d mesh defined by the given point
       *
       * \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
+      * \see updateTriangularMesh
       *
       * \param index dataset index specifying group and dataset to extract value from
       * \param point point to query in map coordinates
@@ -289,6 +296,24 @@ class CORE_EXPORT QgsMeshLayer : public QgsMapLayer
       */
     QgsMesh3dDataBlock dataset3dValue( const QgsMeshDatasetIndex &index, const QgsPointXY &point ) const;
 
+    /**
+      * Returns the value of 1D mesh dataset defined on edge that are in the search area defined by point ans searchRadius
+      *
+      * \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
+      * \see updateTriangularMesh
+      *
+      * \param index dataset index specifying group and dataset to extract value from
+      * \param point the center point of the search area
+      * \param searchRadius the radius of the searc area in map unit
+      * \returns interpolated value at the projected point. Returns NaN values for values
+      * outside the mesh layer and in case triangular mesh was not previously used for rendering
+      *
+      * \since QGIS 3.14
+      */
+    QgsMeshDatasetValue dataset1dValue( const QgsMeshDatasetIndex &index, const QgsPointXY &point, double searchRadius ) const;
+
     /**
       * Returns dataset index from active scalar group depending on the time range.
       * If the temporal properties is not active, return the static dataset

diff --git a/src/core/mesh/qgsmeshlayerutils.cpp b/src/core/mesh/qgsmeshlayerutils.cpp
@@ -269,6 +269,12 @@ double QgsMeshLayerUtils::interpolateFromVerticesData( double fraction, double v
   return val1 + ( val2 - val1 ) * fraction;
 }
 
+QgsMeshDatasetValue QgsMeshLayerUtils::interpolateFromVerticesData( double fraction, const QgsMeshDatasetValue &val1, const QgsMeshDatasetValue &val2 )
+{
+  return QgsMeshDatasetValue( interpolateFromVerticesData( fraction, val1.x(), val2.x() ),
+                              interpolateFromVerticesData( fraction, val1.y(), val2.y() ) );
+}
+
 
 QgsVector QgsMeshLayerUtils::interpolateVectorFromVerticesData( const QgsPointXY &p1, const QgsPointXY &p2, const QgsPointXY &p3, QgsVector vect1, QgsVector vect2, QgsVector vect3, const QgsPointXY &pt )
 {

diff --git a/src/core/mesh/qgsmeshlayerutils.h b/src/core/mesh/qgsmeshlayerutils.h
@@ -137,6 +137,16 @@ class CORE_EXPORT QgsMeshLayerUtils
       double val1, double val2
     );
 
+    /**
+    * Interpolates value based on known values on the vertices of a edge
+    * \returns value on the point pt a or NaN
+    *
+    * \since QGIS 3.14
+    */
+    static QgsMeshDatasetValue interpolateFromVerticesData( double fraction,
+        const QgsMeshDatasetValue &val1, const QgsMeshDatasetValue &val2
+                                                          );
+
     /**
     * Interpolates value based on known values on the vertices of a triangle
     * \param p1 first vertex of the triangle

diff --git a/src/gui/qgsmaptoolidentify.cpp b/src/gui/qgsmaptoolidentify.cpp
@@ -259,12 +259,14 @@ bool QgsMapToolIdentify::identifyMeshLayer( QList<QgsMapToolIdentify::IdentifyRe
 
   QMap< QString, QString > scalarAttributes, vectorAttributes, raw3dAttributes;
 
+  double searchRadius = mOverrideCanvasSearchRadius < 0 ? searchRadiusMU( mCanvas ) : mOverrideCanvasSearchRadius;
+
   QString scalarGroup;
   if ( scalarDatasetIndex.isValid() )
   {
     scalarGroup = layer->dataProvider()->datasetGroupMetadata( scalarDatasetIndex.group() ).name();
 
-    const QgsMeshDatasetValue scalarValue = layer->datasetValue( scalarDatasetIndex, point );
+    const QgsMeshDatasetValue scalarValue = layer->datasetValue( scalarDatasetIndex, point, searchRadius );
     const double scalar = scalarValue.scalar();
     if ( std::isnan( scalar ) )
       scalarAttributes.insert( tr( "Scalar Value" ), tr( "no data" ) );
@@ -277,7 +279,7 @@ bool QgsMapToolIdentify::identifyMeshLayer( QList<QgsMapToolIdentify::IdentifyRe
   {
     vectorGroup = layer->dataProvider()->datasetGroupMetadata( vectorDatasetIndex.group() ).name();
 
-    const QgsMeshDatasetValue vectorValue = layer->datasetValue( vectorDatasetIndex, point );
+    const QgsMeshDatasetValue vectorValue = layer->datasetValue( vectorDatasetIndex, point, searchRadius );
     const double vectorX = vectorValue.x();
     const double vectorY = vectorValue.y();
 

diff --git a/tests/src/core/testqgsmeshlayer.cpp b/tests/src/core/testqgsmeshlayer.cpp
@@ -82,6 +82,8 @@ class TestQgsMeshLayer : public QObject
     void test_reload_extra_dataset();
 
     void test_mesh_simplification();
+
+    void test_dataset_value_from_layer();
 };
 
 QString TestQgsMeshLayer::readFile( const QString &fname ) const
@@ -909,6 +911,38 @@ void TestQgsMeshLayer::test_mesh_simplification()
     delete m;
 }
 
+void TestQgsMeshLayer::test_dataset_value_from_layer()
+{
+  QgsMeshDatasetValue value;
+
+  //1D mesh
+  mMdal1DLayer->updateTriangularMesh();
+  value = mMdal1DLayer->datasetValue( QgsMeshDatasetIndex( 0, 0 ), QgsPointXY( 1500, 2009 ), 10 );
+  QCOMPARE( QgsMeshDatasetValue( 25 ), value );
+  value = mMdal1DLayer->datasetValue( QgsMeshDatasetIndex( 1, 0 ), QgsPointXY( 2500, 1991 ), 10 );
+  QCOMPARE( QgsMeshDatasetValue( 2.5 ), value );
+  value = mMdal1DLayer->datasetValue( QgsMeshDatasetIndex( 2, 0 ), QgsPointXY( 2500, 2500 ), 10 );
+  QCOMPARE( QgsMeshDatasetValue( 2.5, 2 ), value );
+  value = mMdal1DLayer->datasetValue( QgsMeshDatasetIndex( 3, 1 ), QgsPointXY( 2495, 2000 ), 10 );
+  QCOMPARE( QgsMeshDatasetValue( 3 ), value );
+  value = mMdal1DLayer->datasetValue( QgsMeshDatasetIndex( 4, 1 ), QgsPointXY( 2500, 2495 ), 10 );
+  QCOMPARE( QgsMeshDatasetValue( 4, 4 ), value );
+
+  //2D mesh
+  mMdalLayer->updateTriangularMesh();
+  value = mMdalLayer->datasetValue( QgsMeshDatasetIndex( 0, 0 ), QgsPointXY( 1750, 2250 ) );
+  QCOMPARE( QgsMeshDatasetValue( 32.5 ), value );
+  value = mMdalLayer->datasetValue( QgsMeshDatasetIndex( 1, 0 ), QgsPointXY( 1750, 2250 ) );
+  QCOMPARE( QgsMeshDatasetValue( 1.75 ), value );
+  value = mMdalLayer->datasetValue( QgsMeshDatasetIndex( 2, 0 ), QgsPointXY( 1750, 2250 ) );
+  QCOMPARE( QgsMeshDatasetValue( 1.75, 1.25 ), value );
+  value = mMdalLayer->datasetValue( QgsMeshDatasetIndex( 3, 1 ), QgsPointXY( 1750, 2250 ) );
+  QCOMPARE( QgsMeshDatasetValue( 2 ), value );
+  value = mMdalLayer->datasetValue( QgsMeshDatasetIndex( 4, 1 ), QgsPointXY( 1750, 2250 ) );
+  QCOMPARE( QgsMeshDatasetValue( 2, 2 ), value );
+
+}
+
 void TestQgsMeshLayer::test_temporal()
 {
   qint64 relativeTime_0 = -1000;