[mesh] [feature] function to get value for the point on map #7582

PeterPetrik · 2018-08-09T14:51:11Z

To be able to create some nice plots, we need a way how to get a (scalar) value from a dataset for particular point. The face (triangle) of the that contains point is identified by spatial index. Value is interpolated similarly to raster (contours) rendering of the mesh layer.

Probably the plots for mesh layers will not make it to QGIS 3.4, but see the example of the usage of the new function from python (Crayfish library for QGIS 3.4, not released yet = WIP)

m-kuhn · 2018-08-09T15:30:53Z

src/core/mesh/qgstriangularmesh.h

@@ -77,6 +80,12 @@ class CORE_EXPORT QgsTriangularMesh
    //! Returns mapping between triangles and original faces
    const QVector<int> &trianglesToNativeFaces() const ;

+    /**
+     * Returns triangle index that contains the given point, -1 if no such triangle exists
+     * It uses spatial indexing


m-kuhn · 2018-08-09T15:32:06Z

src/core/mesh/qgstriangularmesh.cpp

@@ -143,6 +146,15 @@ 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 )


Not sure what kind of datatype mTrisangularMesh.faces is, but if it's a container, better use an iterator

I also need i few lines lower to add geometry to spatial index

If it's performance critical code, better iterate and use a separate counter. If not, no worries.

m-kuhn · 2018-08-09T15:33:42Z

src/core/mesh/qgsmeshlayerinterpolator.h

@@ -59,6 +59,12 @@ class QgsMeshLayerInterpolator : public QgsRasterInterface
    int bandCount() const override;
    QgsRasterBlock *block( int, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback = nullptr ) override;

+
+    static double interpolateFromVerticesData( const QgsPointXY &p1,


I think this deserves some API docs too

m-kuhn · 2018-08-09T15:34:41Z

src/core/mesh/qgsmeshlayer.h

+     *
+     * \since QGIS 3.4
+     */
+    void activeScalarDatasetChanged( QgsMeshDatasetIndex index );


Is that a candidate for a const reference or a super trivial structure?

it is 2 integers

In other places it's passed as const ref (e.g. QgsMeshDatasetValue QgsMeshLayer::datasetValue( const QgsMeshDatasetIndex &index, const QgsPointXY &point ) const), @nyalldawson I think you once investigated when it's appropriate to copy vs. ref?

nyalldawson · 2018-08-09T23:34:12Z

src/core/mesh/qgsmeshlayer.cpp

+      bool isOnFaces = dataProvider()->datasetGroupMetadata( index ).isOnFaces();
+      if ( isOnFaces )
+      {
+        int native_face_index = mTriangularMesh->trianglesToNativeFaces().at( face_index );


native_face_index -> camel case

nyalldawson · 2018-08-09T23:34:52Z

src/core/mesh/qgsmeshlayer.h

+      * Returns NaN for NaN values, for values outside range or when
+      * triangular mesh is not yet initialized on given point
+      */
+    QgsMeshDatasetValue datasetValue( const QgsMeshDatasetIndex &index, const QgsPointXY &point ) const;


\since QGIS 3.4

nyalldawson · 2018-08-09T23:35:43Z

src/core/mesh/qgstriangularmesh.cpp

+  {
+    const QgsMeshFace &face = mTriangularMesh.faces.at( i ) ;
+    QgsGeometry geom = QgsMeshUtils::toGeometry( face, mTriangularMesh.vertices );
+    bool success = mSpatialIndex.insertFeature( i, geom.boundingBox() );


instead of this just go

(void)mSpatialIndex.insertFeature( i, geom.boundingBox() );

nyalldawson · 2018-08-09T23:35:58Z

src/core/mesh/qgstriangularmesh.cpp

+  const QList<QgsFeatureId> face_indexes = mSpatialIndex.intersects( QgsRectangle( point, point ) );
+  for ( QgsFeatureId fid : face_indexes )
+  {
+    int face_index = static_cast<int>( fid );


nyalldawson · 2018-08-09T23:38:56Z

src/core/mesh/qgstriangularmesh.cpp

+
+QgsGeometry QgsMeshUtils::toGeometry( const QgsMeshFace &face, const QVector<QgsMeshVertex> &vertices )
+{
+  QVector<QgsPointXY> ring;


Would be more efficient to avoid the QgsPointXY methods here. Instead I'd suggest

QVector<QgsPoint> ring; ... std::unique_ptr< QgsPolygon > polygon = qgis::make_unique< QgsPolygon >(); polygon.setExteriorRing( new QgsLineString( ring ) ); return QgsGeometry( std::move( polygon ) );

nyalldawson · 2018-08-09T23:39:06Z

src/core/mesh/qgstriangularmesh.cpp

+  QVector<QgsPointXY> ring;
+  for ( int j = 0; j < face.size(); ++j )
+  {
+    int vertex_id = face[j];


nirvn · 2018-08-10T06:52:10Z

@PeterPetrik , are you planning to add mesh layer support to the identify map tool? Users will probably expect it to work.

PeterPetrik · 2018-08-10T07:02:15Z

@nirvn yes, but we would like to do it as a separate PR

wonder-sk · 2018-08-10T07:49:29Z

src/core/mesh/qgsmeshlayer.h

+      *
+      * \since QGIS 3.4
+      */
+    QgsMeshDatasetValue datasetValue( const QgsMeshDatasetIndex &index, const QgsPointXY &point ) const;


The "point" argument is in map coordinates rather than layer coordinates, right? That's worth noting in documentation (the assumption is that we pass layer coordinates)

Also, worth mentioning that it uses cached previously used triangular mesh and so if the layer has not been rendered previously (e.g. when used in a script) it would return invalid value

wonder-sk · 2018-08-10T07:51:40Z

src/core/mesh/qgstriangularmesh.cpp

@@ -165,3 +178,31 @@ const QVector<int> &QgsTriangularMesh::trianglesToNativeFaces() const
  return mTrianglesToNativeFaces;
 }

+int QgsTriangularMesh::faceIndexForPoint( const QgsPointXY &point ) const
+{
+  const QList<QgsFeatureId> face_indexes = mSpatialIndex.intersects( QgsRectangle( point, point ) );


wonder-sk · 2018-08-10T07:51:53Z

src/core/mesh/qgsmeshlayer.cpp

+
+  if ( mTriangularMesh && dataProvider() && dataProvider()->isValid() && index.isValid() )
+  {
+    int face_index = mTriangularMesh->faceIndexForPoint( point ) ;


wonder-sk · 2018-08-10T07:53:07Z

src/core/mesh/qgsmeshdataprovider.h

+    /**
+     * \brief Returns whether dataset group data is defined on faces
+     */
+    bool isOnFaces() const;


Maybe it would be nicer to replace isOnVertices() / isOnFaces() with a "type" enum. It would be easier to understand that "on vertices" and "on faces" are mutually exclusive

wonder-sk · 2018-08-10T13:45:33Z

Tests fail only due to recently broken dbmanager test on master in f4d08eb

[mesh] [feature] add function to identify value on the point

c79e1d0

PeterPetrik mentioned this pull request Aug 9, 2018

qgis3 port lutraconsulting/qgis-crayfish-plugin#345

Merged

m-kuhn reviewed Aug 9, 2018

View reviewed changes

nyalldawson reviewed Aug 9, 2018

View reviewed changes

nyalldawson added the Feature label Aug 10, 2018

nyalldawson added this to the 3.4 milestone Aug 10, 2018

format code and add documentation

9ea0bc2

wonder-sk reviewed Aug 10, 2018

View reviewed changes

introduce DataType for mesh dataset data definition and small fixes

271bab1

wonder-sk approved these changes Aug 10, 2018

View reviewed changes

wonder-sk merged commit ad4ddb1 into qgis:master Aug 10, 2018

qgib mentioned this pull request May 25, 2019

Identify tool does not work on mesh layer #27951

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[mesh] [feature] function to get value for the point on map #7582

[mesh] [feature] function to get value for the point on map #7582

PeterPetrik commented Aug 9, 2018

m-kuhn Aug 9, 2018

PeterPetrik Aug 10, 2018

m-kuhn Aug 9, 2018

PeterPetrik Aug 10, 2018

m-kuhn Aug 10, 2018

m-kuhn Aug 9, 2018

m-kuhn Aug 9, 2018

PeterPetrik Aug 10, 2018

m-kuhn Aug 10, 2018

nyalldawson Aug 9, 2018

nyalldawson Aug 9, 2018

nyalldawson Aug 9, 2018

nyalldawson Aug 9, 2018

nyalldawson Aug 9, 2018

nyalldawson Aug 9, 2018

nirvn commented Aug 10, 2018

PeterPetrik commented Aug 10, 2018

wonder-sk Aug 10, 2018

wonder-sk Aug 10, 2018

wonder-sk Aug 10, 2018

wonder-sk Aug 10, 2018

wonder-sk commented Aug 10, 2018

[mesh] [feature] function to get value for the point on map #7582

[mesh] [feature] function to get value for the point on map #7582

Conversation

PeterPetrik commented Aug 9, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

nirvn commented Aug 10, 2018

PeterPetrik commented Aug 10, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

wonder-sk commented Aug 10, 2018