Skip to content
Permalink
Browse files

identify tool for 1d mesh

  • Loading branch information
vcloarec authored and nyalldawson committed Apr 27, 2020
1 parent f7b001f commit a6ef994cd28f5565a1a0f7268107ff732176fe13
@@ -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 )
@@ -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
@@ -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 )
{
@@ -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
@@ -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();

@@ -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;

0 comments on commit a6ef994

Please sign in to comment.
You can’t perform that action at this time.