sharpetronics/QGIS
1
0
Fork 0
mirror of https://github.com/qgis/QGIS.git synced 2025-12-06 00:03:16 -05:00
Code Issues Packages Projects Releases Wiki Activity

Implement profile generation for point vector layers

Browse Source
This commit is contained in:
Nyall Dawson 2022-03-21 15:46:33 +10:00
parent bdc1d28bfc
commit 87c3c911c0
19 changed files with 930 additions and 49 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

58
python/core/auto_generated/elevation/qgsabstractprofilegenerator.sip.in
View File

@ -9,6 +9,45 @@
class QgsAbstractProfileResults
{
%Docstring(signature="appended")
Abstract base class for storage of elevation profiles.
.. versionadded:: 3.26
%End
%TypeHeaderCode
#include "qgsabstractprofilegenerator.h"
%End
public:
virtual ~QgsAbstractProfileResults();
virtual QString type() const = 0;
%Docstring
Returns the unique string identifier for the results type.
%End
virtual QHash< double, double > distanceToHeightMap() const = 0;
%Docstring
Returns the map of distance (chainage) to height.
%End
virtual QgsPointSequence sampledPoints() const = 0;
%Docstring
Returns a list of sampled points, with their calculated elevation
as the point z value.
%End
virtual QVector< QgsGeometry > asGeometries() const = 0;
%Docstring
Returns a list of geometries representing the calculated elevation results.
%End
};
class QgsAbstractProfileGenerator
{
%Docstring(signature="appended")
@ -31,7 +70,7 @@ The scenario will be:
additional calls to the source.
# profile job (still in GUI thread) stores the generator for later use.
# profile job (in worker thread) calls :py:func:`QgsAbstractProfileGenerator.generateProfile()`
# profile job (again in GUI thread) will check :py:func:`~errors` and report them
# profile job (again in GUI thread) will check :py:func:`~QgsAbstractProfileResults.errors` and report them
.. versionadded:: 3.26
%End
@ -51,24 +90,13 @@ Returns ``True`` if the profile was generated successfully (i.e. the generation
was not canceled early).
%End
struct Result
{
double distance;
double height;
};
QList< QgsAbstractProfileGenerator::Result > results() const;
virtual QgsAbstractProfileResults *takeResults() = 0 /TransferBack/;
%Docstring
Temporary method to return results.
%End
Takes results from the generator.
QList< QgsPoint > rawPoints() const;
%Docstring
Temporary method to return results.
Ownership is transferred to the caller.
%End
protected:
};
/************************************************************************

14
python/core/auto_generated/elevation/qgsterrainprovider.sip.in
View File

@ -78,6 +78,14 @@ Returns the unique type ID string for the provider.
Creates a clone of the provider and returns the new object.
Ownership is transferred to the caller.
%End
virtual void prepare() = 0 /Factory/;
%Docstring
Called on the main thread prior to accessing the provider from a background thread.
Subclasses must implement suitable logic in order to prepare for thread-safe calculation of
terrain heights on background threads.
%End
virtual QgsCoordinateReferenceSystem crs() const = 0;
@ -183,6 +191,8 @@ Constructor for QgsFlatTerrainProvider.
virtual QgsFlatTerrainProvider *clone() const /Factory/;
virtual void prepare();
virtual bool equals( const QgsAbstractTerrainProvider *other ) const;
};
@ -221,6 +231,8 @@ Constructor for QgsRasterDemTerrainProvider.
virtual bool equals( const QgsAbstractTerrainProvider *other ) const;
virtual void prepare();
void setLayer( QgsRasterLayer *layer );
%Docstring
@ -274,6 +286,8 @@ Constructor for QgsMeshTerrainProvider.
virtual bool equals( const QgsAbstractTerrainProvider *other ) const;
virtual void prepare();
void setLayer( QgsMeshLayer *layer );
%Docstring

4
python/core/auto_generated/vector/qgsvectorlayer.sip.in
View File

@ -18,7 +18,7 @@ typedef QList<int> QgsAttributeList;
typedef QSet<int> QgsAttributeIds;
class QgsVectorLayer : QgsMapLayer, QgsExpressionContextGenerator, QgsExpressionContextScopeGenerator, QgsFeatureSink, QgsFeatureSource
class QgsVectorLayer : QgsMapLayer, QgsExpressionContextGenerator, QgsExpressionContextScopeGenerator, QgsFeatureSink, QgsFeatureSource, QgsAbstractProfileSource
{
%Docstring(signature="appended")
Represents a vector layer which manages a vector based data sets.
@ -513,6 +513,8 @@ Uses :py:class:`QgsExpression`
virtual QgsMapLayerElevationProperties *elevationProperties();
virtual QgsAbstractProfileGenerator *createProfileGenerator( const QgsProfileRequest &request ) /Factory/;
void setProviderEncoding( const QString &encoding );
%Docstring

82
python/core/conversions.sip
View File

@ -960,6 +960,88 @@ template<double, TYPE>
};
%MappedType QHash<double, double>
{
%TypeHeaderCode
#include <QHash>
%End
%ConvertFromTypeCode
// Create the dictionary.
PyObject *d = PyDict_New();
if (!d)
return NULL;
// Set the dictionary elements.
QHash<double, double>::iterator i;
for (i = sipCpp->begin(); i != sipCpp->end(); ++i)
{
PyObject *t1obj = PyFloat_FromDouble(i.key());
PyObject *t2obj = PyFloat_FromDouble(i.value());
if (t1obj == NULL || t2obj == NULL || PyDict_SetItem(d, t1obj, t2obj) < 0)
{
Py_DECREF(d);
if (t1obj)
Py_DECREF(t1obj);
if (t2obj)
Py_DECREF(t2obj);
return NULL;
}
Py_DECREF(t1obj);
Py_DECREF(t2obj);
}
return d;
%End
%ConvertToTypeCode
PyObject *t1obj, *t2obj;
Py_ssize_t i = 0;
// Check the type if that is all that is required.
if (sipIsErr == NULL)
{
if (!PyDict_Check(sipPy))
return 0;
while (PyDict_Next(sipPy, &i, &t1obj, &t2obj))
{
if (!PyFloat_Check(t1obj))
return 0;
if (!PyFloat_Check(t2obj))
return 0;
}
return 1;
}
QHash<double, double> *qm = new QHash<double, double>;
while (PyDict_Next(sipPy, &i, &t1obj, &t2obj))
{
int state;
double k = PyFloat_AsDouble(t1obj);
double v = PyFloat_AsDouble(t2obj);
qm->insert(k, v);
}
*sipCppPtr = qm;
return sipGetState(sipTransferObj);
%End
};
template<TYPE2>
%MappedType QMap<QString, QList<TYPE2> >
{

2
src/core/CMakeLists.txt
View File

@ -837,6 +837,7 @@ set(QGIS_CORE_SRCS
vector/qgsvectorlayerexporter.cpp
vector/qgsvectorlayerjoinbuffer.cpp
vector/qgsvectorlayerjoininfo.cpp
vector/qgsvectorlayerprofilegenerator.cpp
vector/qgsvectorlayerrenderer.cpp
vector/qgsvectorlayertemporalproperties.cpp
vector/qgsvectorlayertools.cpp
@ -1793,6 +1794,7 @@ set(QGIS_CORE_HDRS
vector/qgsvectorlayerfeatureiterator.h
vector/qgsvectorlayerjoinbuffer.h
vector/qgsvectorlayerjoininfo.h
vector/qgsvectorlayerprofilegenerator.h
vector/qgsvectorlayerrenderer.h
vector/qgsvectorlayertemporalproperties.h
vector/qgsvectorlayertools.h

2
src/core/elevation/qgsabstractprofilegenerator.cpp
View File

@ -17,3 +17,5 @@
#include "qgsabstractprofilegenerator.h"
QgsAbstractProfileGenerator::~QgsAbstractProfileGenerator() = default;
QgsAbstractProfileResults::~QgsAbstractProfileResults() = default;

59
src/core/elevation/qgsabstractprofilegenerator.h
View File

@ -23,6 +23,43 @@
#include "qgspoint.h"
class QgsGeometry;
/**
* \brief Abstract base class for storage of elevation profiles.
*
* \ingroup core
* \since QGIS 3.26
*/
class CORE_EXPORT QgsAbstractProfileResults
{
public:
virtual ~QgsAbstractProfileResults();
/**
* Returns the unique string identifier for the results type.
*/
virtual QString type() const = 0;
/**
* Returns the map of distance (chainage) to height.
*/
virtual QHash< double, double > distanceToHeightMap() const = 0;
/**
* Returns a list of sampled points, with their calculated elevation
* as the point z value.
*/
virtual QgsPointSequence sampledPoints() const = 0;
/**
* Returns a list of geometries representing the calculated elevation results.
*/
virtual QVector< QgsGeometry > asGeometries() const = 0;
};
/**
* \brief Abstract base class for objects which generate elevation profiles.
*
@ -63,26 +100,12 @@ class CORE_EXPORT QgsAbstractProfileGenerator
*/
virtual bool generateProfile() = 0;
// Temporary class only!
struct CORE_EXPORT Result
{
double distance;
double height;
};
/**
* Temporary method to return results.
* Takes results from the generator.
*
* Ownership is transferred to the caller.
*/
QList< QgsAbstractProfileGenerator::Result > results() const { return mResults; }
/**
* Temporary method to return results.
*/
QList< QgsPoint > rawPoints() const { return mRawPoints; }
protected:
QList< QgsPoint > mRawPoints;
QList< Result > mResults;
virtual QgsAbstractProfileResults *takeResults() = 0 SIP_TRANSFERBACK;
};

1
src/core/elevation/qgsprofilerequest.h
View File

@ -18,7 +18,6 @@
#define QGSPROFILEREQUEST_H
#include "qgis_core.h"
#include "qgis_sip.h"
#include "qgscoordinatereferencesystem.h"
#include "qgscoordinatetransformcontext.h"

48
src/core/elevation/qgsterrainprovider.cpp
View File

@ -15,6 +15,7 @@
* *
***************************************************************************/
#include "qgsterrainprovider.h"
#include <QThread>
QgsAbstractTerrainProvider::~QgsAbstractTerrainProvider() = default;
@ -83,6 +84,12 @@ QgsFlatTerrainProvider *QgsFlatTerrainProvider::clone() const
return new QgsFlatTerrainProvider( *this );
}
void QgsFlatTerrainProvider::prepare()
{
Q_ASSERT_X( QThread::currentThread() == QCoreApplication::instance()->thread(), "QgsFlatTerrainProvider::prepare", "prepare() must be called from the main thread" );
}
bool QgsFlatTerrainProvider::equals( const QgsAbstractTerrainProvider *other ) const
{
if ( other->type() != type() )
@ -144,21 +151,27 @@ QDomElement QgsRasterDemTerrainProvider::writeXml( QDomDocument &document, const
QgsCoordinateReferenceSystem QgsRasterDemTerrainProvider::crs() const
{
return mRasterLayer ? mRasterLayer->crs() : QgsCoordinateReferenceSystem();
return mRasterProvider ? mRasterProvider->crs()
: ( mRasterLayer ? mRasterLayer->crs() : QgsCoordinateReferenceSystem() );
}
double QgsRasterDemTerrainProvider::heightAt( double x, double y ) const
{
// TODO -- may want to use a more efficient approach here, i.e. requesting whole
// blocks upfront instead of multiple sample calls
if ( mRasterLayer && mRasterLayer->isValid() )
bool ok = false;
double res = std::numeric_limits<double>::quiet_NaN();
if ( mRasterProvider )
{
bool ok = false;
const double res = mRasterLayer->dataProvider()->sample( QgsPointXY( x, y ), 1, &ok );
if ( ok )
return res * mScale + mOffset;
res = mRasterProvider->sample( QgsPointXY( x, y ), 1, &ok );
}
else if ( QThread::currentThread() == QCoreApplication::instance()->thread() && mRasterLayer && mRasterLayer->isValid() )
{
res = mRasterLayer->dataProvider()->sample( QgsPointXY( x, y ), 1, &ok );
}
if ( ok )
return res * mScale + mOffset;
return std::numeric_limits<double>::quiet_NaN();
}
@ -182,6 +195,14 @@ bool QgsRasterDemTerrainProvider::equals( const QgsAbstractTerrainProvider *othe
return true;
}
void QgsRasterDemTerrainProvider::prepare()
{
Q_ASSERT_X( QThread::currentThread() == QCoreApplication::instance()->thread(), "QgsRasterDemTerrainProvider::prepare", "prepare() must be called from the main thread" );
if ( mRasterLayer && mRasterLayer->isValid() )
mRasterProvider.reset( mRasterLayer->dataProvider()->clone() );
}
void QgsRasterDemTerrainProvider::setLayer( QgsRasterLayer *layer )
{
mRasterLayer.setLayer( layer );
@ -192,6 +213,13 @@ QgsRasterLayer *QgsRasterDemTerrainProvider::layer() const
return mRasterLayer.get();
}
QgsRasterDemTerrainProvider::QgsRasterDemTerrainProvider( const QgsRasterDemTerrainProvider &other )
: QgsAbstractTerrainProvider( other )
, mRasterLayer( other.mRasterLayer )
{
}
//
// QgsMeshTerrainProvider
@ -273,6 +301,12 @@ bool QgsMeshTerrainProvider::equals( const QgsAbstractTerrainProvider *other ) c
return true;
}
void QgsMeshTerrainProvider::prepare()
{
Q_ASSERT_X( QThread::currentThread() == QCoreApplication::instance()->thread(), "QgsMeshTerrainProvider::prepare", "prepare() must be called from the main thread" );
}
void QgsMeshTerrainProvider::setLayer( QgsMeshLayer *layer )
{
mMeshLayer.setLayer( layer );

13
src/core/elevation/qgsterrainprovider.h
View File

@ -103,6 +103,14 @@ class CORE_EXPORT QgsAbstractTerrainProvider
*/
virtual QgsAbstractTerrainProvider *clone() const = 0 SIP_FACTORY;
/**
* Called on the main thread prior to accessing the provider from a background thread.
*
* Subclasses must implement suitable logic in order to prepare for thread-safe calculation of
* terrain heights on background threads.
*/
virtual void prepare() = 0 SIP_FACTORY;
/**
* Returns the native coordinate reference system of the terrain provider.
*/
@ -208,6 +216,7 @@ class CORE_EXPORT QgsFlatTerrainProvider : public QgsAbstractTerrainProvider
QgsCoordinateReferenceSystem crs() const override;
double heightAt( double x, double y ) const override;
QgsFlatTerrainProvider *clone() const override SIP_FACTORY;
void prepare() override;
bool equals( const QgsAbstractTerrainProvider *other ) const override;
};
@ -234,6 +243,7 @@ class CORE_EXPORT QgsRasterDemTerrainProvider : public QgsAbstractTerrainProvide
double heightAt( double x, double y ) const override;
QgsRasterDemTerrainProvider *clone() const override SIP_FACTORY;
bool equals( const QgsAbstractTerrainProvider *other ) const override;
void prepare() override;
/**
* Sets the raster \a layer with elevation model to be used as the terrain source.
@ -251,8 +261,10 @@ class CORE_EXPORT QgsRasterDemTerrainProvider : public QgsAbstractTerrainProvide
private:
QgsRasterDemTerrainProvider( const QgsRasterDemTerrainProvider &other );
_LayerRef<QgsRasterLayer> mRasterLayer;
std::unique_ptr< QgsRasterDataProvider > mRasterProvider;
};
@ -280,6 +292,7 @@ class CORE_EXPORT QgsMeshTerrainProvider : public QgsAbstractTerrainProvider
double heightAt( double x, double y ) const override;
QgsMeshTerrainProvider *clone() const override SIP_FACTORY;
bool equals( const QgsAbstractTerrainProvider *other ) const override;
void prepare() override;
/**
* Sets the mesh \a layer to be used as the terrain source.

56
src/core/raster/qgsrasterlayerprofilegenerator.cpp
View File

@ -23,6 +23,45 @@
#include "qgsgeometryengine.h"
#include "qgsgeos.h"
//
// QgsRasterLayerProfileResults
//
QString QgsRasterLayerProfileResults::type() const
{
return QStringLiteral( "raster" );
}
QHash<double, double> QgsRasterLayerProfileResults::distanceToHeightMap() const
{
QHash<double, double> res;
for ( const Result &r : results )
{
res.insert( r.distance, r.height );
}
return res;
}
QgsPointSequence QgsRasterLayerProfileResults::sampledPoints() const
{
return rawPoints;
}
QVector<QgsGeometry> QgsRasterLayerProfileResults::asGeometries() const
{
QVector<QgsGeometry> res;
res.reserve( rawPoints.size() );
for ( const QgsPoint &point : rawPoints )
res.append( QgsGeometry( point.clone() ) );
return res;
}
//
// QgsRasterLayerProfileGenerator
//
QgsRasterLayerProfileGenerator::QgsRasterLayerProfileGenerator( QgsRasterLayer *layer, const QgsProfileRequest &request )
: mProfileCurve( request.profileCurve() ? request.profileCurve()->clone() : nullptr )
, mSourceCrs( layer->crs() )
@ -60,6 +99,8 @@ bool QgsRasterLayerProfileGenerator::generateProfile()
if ( !profileCurveBoundingBox.intersects( mRasterProvider->extent() ) )
return false;
mResults = std::make_unique< QgsRasterLayerProfileResults >();
std::unique_ptr< QgsGeometryEngine > curveEngine( QgsGeometry::createGeometryEngine( transformedCurve.get() ) );
curveEngine->prepareGeometry();
@ -136,7 +177,7 @@ bool QgsRasterLayerProfileGenerator::generateProfile()
{
continue;
}
mRawPoints.append( pixel );
mResults->rawPoints.append( pixel );
}
currentY -= mRasterUnitsPerPixelY;
}
@ -145,17 +186,22 @@ bool QgsRasterLayerProfileGenerator::generateProfile()
// convert x/y values back to distance/height values
QgsGeos originalCurveGeos( mProfileCurve.get() );
originalCurveGeos.prepareGeometry();
mResults.reserve( mRawPoints.size() );
mResults->results.reserve( mResults->rawPoints.size() );
QString lastError;
for ( const QgsPoint &pixel : std::as_const( mRawPoints ) )
for ( const QgsPoint &pixel : std::as_const( mResults->rawPoints ) )
{
const double distance = originalCurveGeos.lineLocatePoint( pixel, &lastError );
Result res;
QgsRasterLayerProfileResults::Result res;
res.distance = distance;
res.height = pixel.z();
mResults.push_back( res );
mResults->results.push_back( res );
}
return true;
}
QgsAbstractProfileResults *QgsRasterLayerProfileGenerator::takeResults()
{
return mResults.release();
}

30
src/core/raster/qgsrasterlayerprofilegenerator.h
View File

@ -33,6 +33,33 @@ class QgsRasterDataProvider;
#define SIP_NO_FILE
/**
* \brief Implementation of QgsAbstractProfileResults for raster layers.
*
* \note Not available in Python bindings
* \ingroup core
* \since QGIS 3.26
*/
class CORE_EXPORT QgsRasterLayerProfileResults : public QgsAbstractProfileResults
{
public:
struct Result
{
double distance;
double height;
};
QgsPointSequence rawPoints;
QList< Result > results;
QString type() const override;
QHash< double, double > distanceToHeightMap() const override;
QgsPointSequence sampledPoints() const override;
QVector< QgsGeometry > asGeometries() const override;
};
/**
* \brief Implementation of QgsAbstractProfileGenerator for raster layers.
*
@ -53,6 +80,7 @@ class CORE_EXPORT QgsRasterLayerProfileGenerator : public QgsAbstractProfileGene
~QgsRasterLayerProfileGenerator() override;
bool generateProfile() override;
QgsAbstractProfileResults *takeResults() override;
private:
@ -67,6 +95,8 @@ class CORE_EXPORT QgsRasterLayerProfileGenerator : public QgsAbstractProfileGene
std::unique_ptr< QgsRasterDataProvider > mRasterProvider;
std::unique_ptr< QgsRasterLayerProfileResults > mResults;
int mBand = 1;
double mRasterUnitsPerPixelX = 1;
double mRasterUnitsPerPixelY = 1;

6
src/core/vector/qgsvectorlayer.cpp
View File

@ -89,6 +89,7 @@
#include "qgsruntimeprofiler.h"
#include "qgsfeaturerenderergenerator.h"
#include "qgsvectorlayerutils.h"
#include "qgsvectorlayerprofilegenerator.h"
#include "diagram/qgsdiagram.h"
@ -683,6 +684,11 @@ QgsMapLayerElevationProperties *QgsVectorLayer::elevationProperties()
return mElevationProperties;
}
QgsAbstractProfileGenerator *QgsVectorLayer::createProfileGenerator( const QgsProfileRequest &request )
{
return new QgsVectorLayerProfileGenerator( this, request );
}
void QgsVectorLayer::setProviderEncoding( const QString &encoding )
{
if ( isValid() && mDataProvider && mDataProvider->encoding() != encoding )

4
src/core/vector/qgsvectorlayer.h
View File

@ -42,6 +42,7 @@
#include "qgsexpressioncontextgenerator.h"
#include "qgsexpressioncontextscopegenerator.h"
#include "qgsexpressioncontext.h"
#include "qgsabstractprofilesource.h"
class QPainter;
class QImage;
@ -389,7 +390,7 @@ typedef QSet<int> QgsAttributeIds;
*
* \see QgsVectorLayerUtils()
*/
class CORE_EXPORT QgsVectorLayer : public QgsMapLayer, public QgsExpressionContextGenerator, public QgsExpressionContextScopeGenerator, public QgsFeatureSink, public QgsFeatureSource
class CORE_EXPORT QgsVectorLayer : public QgsMapLayer, public QgsExpressionContextGenerator, public QgsExpressionContextScopeGenerator, public QgsFeatureSink, public QgsFeatureSource, public QgsAbstractProfileSource
{
Q_OBJECT
@ -627,6 +628,7 @@ class CORE_EXPORT QgsVectorLayer : public QgsMapLayer, public QgsExpressionConte
const QgsVectorDataProvider *dataProvider() const FINAL SIP_SKIP;
QgsMapLayerTemporalProperties *temporalProperties() override;
QgsMapLayerElevationProperties *elevationProperties() override;
QgsAbstractProfileGenerator *createProfileGenerator( const QgsProfileRequest &request ) override SIP_FACTORY;
/**
* Sets the text \a encoding of the data provider.

240
src/core/vector/qgsvectorlayerprofilegenerator.cpp Normal file
View File

@ -0,0 +1,240 @@
/***************************************************************************
qgsvectorlayerprofilegenerator.cpp
---------------
begin : March 2022
copyright : (C) 2022 by Nyall Dawson
email : nyall dot dawson 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 "qgsvectorlayerprofilegenerator.h"
#include "qgsprofilerequest.h"
#include "qgscurve.h"
#include "qgsvectorlayer.h"
#include "qgsvectorlayerelevationproperties.h"
#include "qgscoordinatetransform.h"
#include "qgsgeos.h"
#include "qgsvectorlayerfeatureiterator.h"
#include "qgsterrainprovider.h"
//
// QgsVectorLayerProfileResults
//
QString QgsVectorLayerProfileResults::type() const
{
return QStringLiteral( "vector" );
}
QHash<double, double> QgsVectorLayerProfileResults::distanceToHeightMap() const
{
QHash<double, double> res;
for ( const Result &r : results )
{
res.insert( r.distance, r.height );
}
return res;
}
QgsPointSequence QgsVectorLayerProfileResults::sampledPoints() const
{
return rawPoints;
}
QVector<QgsGeometry> QgsVectorLayerProfileResults::asGeometries() const
{
return geometries;
}
//
// QgsVectorLayerProfileGenerator
//
QgsVectorLayerProfileGenerator::QgsVectorLayerProfileGenerator( QgsVectorLayer *layer, const QgsProfileRequest &request )
: mProfileCurve( request.profileCurve() ? request.profileCurve()->clone() : nullptr )
, mTerrainProvider( request.terrainProvider() ? request.terrainProvider()->clone() : nullptr )
, mTolerance( request.tolerance() )
, mSourceCrs( layer->crs() )
, mTargetCrs( request.crs() )
, mTransformContext( request.transformContext() )
, mSource( std::make_unique< QgsVectorLayerFeatureSource >( layer ) )
, mOffset( layer->elevationProperties()->zOffset() )
, mScale( layer->elevationProperties()->zScale() )
, mClamping( qgis::down_cast< QgsVectorLayerElevationProperties* >( layer->elevationProperties() )->clamping() )
, mBinding( qgis::down_cast< QgsVectorLayerElevationProperties* >( layer->elevationProperties() )->binding() )
, mExtrusionEnabled( qgis::down_cast< QgsVectorLayerElevationProperties* >( layer->elevationProperties() )->extrusionEnabled() )
, mExtrusionHeight( qgis::down_cast< QgsVectorLayerElevationProperties* >( layer->elevationProperties() )->extrusionHeight() )
, mWkbType( layer->wkbType() )
{
if ( mTerrainProvider )
mTerrainProvider->prepare(); // must be done on main thread
}
QgsVectorLayerProfileGenerator::~QgsVectorLayerProfileGenerator() = default;
bool QgsVectorLayerProfileGenerator::generateProfile()
{
if ( !mProfileCurve )
return false;
// we need to transform the profile curve to the vector's CRS
mTransformedCurve.reset( mProfileCurve->clone() );
mLayerToTargetTransform = QgsCoordinateTransform( mSourceCrs, mTargetCrs, mTransformContext );
if ( mTerrainProvider )
mTargetToTerrainProviderTransform = QgsCoordinateTransform( mTargetCrs, mTerrainProvider->crs(), mTransformContext );
try
{
mTransformedCurve->transform( mLayerToTargetTransform, Qgis::TransformDirection::Reverse );
}
catch ( QgsCsException & )
{
QgsDebugMsg( QStringLiteral( "Error transforming profile line to vector CRS" ) );
return false;
}
mResults = std::make_unique< QgsVectorLayerProfileResults >();
switch ( QgsWkbTypes::geometryType( mWkbType ) )
{
case QgsWkbTypes::PointGeometry:
if ( !generateProfileForPoints() )
return false;
break;
}
const QgsRectangle profileCurveBoundingBox = mTransformedCurve->boundingBox();
//if ( !profileCurveBoundingBox.intersects( mRasterProvider->extent() ) )
// return false;
std::unique_ptr< QgsGeometryEngine > curveEngine( QgsGeometry::createGeometryEngine( mTransformedCurve.get() ) );
curveEngine->prepareGeometry();
// convert x/y values back to distance/height values
QgsGeos originalCurveGeos( mProfileCurve.get() );
originalCurveGeos.prepareGeometry();
mResults->results.reserve( mResults->rawPoints.size() );
QString lastError;
for ( const QgsPoint &pixel : std::as_const( mResults->rawPoints ) )
{
const double distance = originalCurveGeos.lineLocatePoint( pixel, &lastError );
QgsVectorLayerProfileResults::Result res;
res.distance = distance;
res.height = pixel.z();
mResults->results.push_back( res );
}
return true;
}
QgsAbstractProfileResults *QgsVectorLayerProfileGenerator::takeResults()
{
return mResults.release();
}
bool QgsVectorLayerProfileGenerator::generateProfileForPoints()
{
// get features from layer
QgsFeatureRequest request;
request.setDestinationCrs( mTargetCrs, mTransformContext );
request.setDistanceWithin( QgsGeometry( mProfileCurve->clone() ), mTolerance );
auto processPoint = [this]( const QgsPoint * point )
{
const double height = featureZToHeight( point->x(), point->y(), point->z() );
mResults->rawPoints.append( QgsPoint( point->x(), point->y(), height ) );
if ( mExtrusionEnabled )
{
mResults->geometries.append( QgsGeometry( new QgsLineString( QgsPoint( point->x(), point->y(), height ),
QgsPoint( point->x(), point->y(), height + mExtrusionHeight ) ) ) );
}
else
{
mResults->geometries.append( QgsGeometry( new QgsPoint( point->x(), point->y(), height ) ) );
}
};
QgsFeature feature;
QgsFeatureIterator it = mSource->getFeatures( request );
while ( it.nextFeature( feature ) )
{
const QgsGeometry g = feature.geometry();
if ( g.isMultipart() )
{
for ( auto it = g.const_parts_begin(); it != g.const_parts_end(); ++it )
{
processPoint( qgsgeometry_cast< const QgsPoint * >( *it ) );
}
}
else
{
processPoint( qgsgeometry_cast< const QgsPoint * >( g.constGet() ) );
}
}
return true;
}
double QgsVectorLayerProfileGenerator::featureZToHeight( double x, double y, double z )
{
switch ( mClamping )
{
case Qgis::AltitudeClamping::Absolute:
break;
case Qgis::AltitudeClamping::Relative:
case Qgis::AltitudeClamping::Terrain:
{
if ( !mTerrainProvider )
break;
// transform feature point to terrain provider crs
try
{
double dummyZ = 0;
mTargetToTerrainProviderTransform.transformInPlace( x, y, dummyZ );
}
catch ( QgsCsException & )
{
break;
}
const double terrainZ = mTerrainProvider->heightAt( x, y );
if ( !std::isnan( terrainZ ) )
{
switch ( mClamping )
{
case Qgis::AltitudeClamping::Relative:
if ( std::isnan( z ) )
z = terrainZ;
else
z += terrainZ;
break;
case Qgis::AltitudeClamping::Terrain:
z = terrainZ;
break;
case Qgis::AltitudeClamping::Absolute:
break;
}
}
break;
}
}
return z * mScale + mOffset;
}

125
src/core/vector/qgsvectorlayerprofilegenerator.h Normal file
View File

@ -0,0 +1,125 @@
/***************************************************************************
qgsvectorlayerprofilegenerator.h
---------------
begin : March 2022
copyright : (C) 2022 by Nyall Dawson
email : nyall dot dawson 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 QGSVECTORLAYERPROFILEGENERATOR_H
#define QGSVECTORLAYERPROFILEGENERATOR_H
#include "qgis_core.h"
#include "qgis_sip.h"
#include "qgsabstractprofilegenerator.h"
#include "qgscoordinatereferencesystem.h"
#include "qgscoordinatetransformcontext.h"
#include "qgscoordinatetransform.h"
#include <memory>
class QgsProfileRequest;
class QgsCurve;
class QgsVectorLayer;
class QgsVectorLayerFeatureSource;
class QgsAbstractTerrainProvider;
#define SIP_NO_FILE
/**
* \brief Implementation of QgsAbstractProfileResults for vector layers.
*
* \note Not available in Python bindings
* \ingroup core
* \since QGIS 3.26
*/
class CORE_EXPORT QgsVectorLayerProfileResults : public QgsAbstractProfileResults
{
public:
// Temporary class only!
struct Result
{
double distance;
double height;
};
QgsPointSequence rawPoints;
QList< Result > results;
QVector< QgsGeometry > geometries;
QString type() const override;
QHash< double, double > distanceToHeightMap() const override;
QgsPointSequence sampledPoints() const override;
QVector< QgsGeometry > asGeometries() const override;
};
/**
* \brief Implementation of QgsAbstractProfileGenerator for vector layers.
*
* \note Not available in Python bindings
* \ingroup core
* \since QGIS 3.26
*/
class CORE_EXPORT QgsVectorLayerProfileGenerator : public QgsAbstractProfileGenerator
{
public:
/**
* Constructor for QgsVectorLayerProfileGenerator.
*/
QgsVectorLayerProfileGenerator( QgsVectorLayer *layer, const QgsProfileRequest &request );
~QgsVectorLayerProfileGenerator() override;
bool generateProfile() override;
QgsAbstractProfileResults *takeResults() override;
private:
bool generateProfileForPoints();
double featureZToHeight( double x, double y, double z );
std::unique_ptr< QgsCurve > mProfileCurve;
std::unique_ptr< QgsAbstractTerrainProvider > mTerrainProvider;
std::unique_ptr< QgsCurve > mTransformedCurve;
double mTolerance = 0;
QgsCoordinateReferenceSystem mSourceCrs;
QgsCoordinateReferenceSystem mTargetCrs;
QgsCoordinateTransformContext mTransformContext;
std::unique_ptr< QgsVectorLayerFeatureSource > mSource;
double mOffset = 0;
double mScale = 1;
Qgis::AltitudeClamping mClamping = Qgis::AltitudeClamping::Terrain;
Qgis::AltitudeBinding mBinding = Qgis::AltitudeBinding::Centroid;
bool mExtrusionEnabled = false;
double mExtrusionHeight = 0;
QgsWkbTypes::Type mWkbType = QgsWkbTypes::Unknown;
QgsCoordinateTransform mLayerToTargetTransform;
QgsCoordinateTransform mTargetToTerrainProviderTransform;
std::unique_ptr< QgsVectorLayerProfileResults > mResults;
};
#endif // QGSVECTORLAYERPROFILEGENERATOR_H

1
tests/src/python/CMakeLists.txt
View File

@ -362,6 +362,7 @@ ADD_PYTHON_TEST(PyQgsVectorLayerCache test_qgsvectorlayercache.py)
ADD_PYTHON_TEST(PyQgsVectorLayerEditBuffer test_qgsvectorlayereditbuffer.py)
ADD_PYTHON_TEST(PyQgsVectorLayerEditBufferGroup test_qgsvectorlayereditbuffergroup.py)
ADD_PYTHON_TEST(PyQgsVectorLayerNamedStyle test_qgsvectorlayer_namedstyle.py)
ADD_PYTHON_TEST(PyQgsVectorLayerProfileGenerator test_qgsvectorlayerprofilegenerator.py)
ADD_PYTHON_TEST(PyQgsVectorLayerRenderer test_qgsvectorlayerrenderer.py)
ADD_PYTHON_TEST(PyQgsVectorLayerSelectedFeatureSource test_qgsvectorlayerselectedfeaturesource.py)
ADD_PYTHON_TEST(PyQgsVectorLayerShapefile test_qgsvectorlayershapefile.py)

3
tests/src/python/test_qgsrasterlayerprofilegenerator.py
View File

@ -60,7 +60,8 @@ class TestQgsRasterLayerProfileGenerator(unittest.TestCase):
generator = rl.createProfileGenerator(req)
self.assertTrue(generator.generateProfile())
results = {r.distance: r.height for r in generator.results()}
r = generator.takeResults()
results = r.distanceToHeightMap()
first_point = min(results.keys())
last_point = max(results.keys())
self.assertEqual(results[first_point], 154)

231
tests/src/python/test_qgsvectorlayerprofilegenerator.py Normal file
View File

@ -0,0 +1,231 @@
# -*- coding: utf-8 -*-
"""QGIS Unit tests for QgsVectorLayer profile generation
.. note:: 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.
"""
__author__ = 'Nyall Dawson'
__date__ = '18/03/2022'
__copyright__ = 'Copyright 2022, The QGIS Project'
import os
import qgis # NOQA
from qgis.core import (
QgsRasterLayer,
QgsLineString,
QgsProfileRequest,
QgsCoordinateReferenceSystem,
QgsVectorLayer,
Qgis,
QgsRasterDemTerrainProvider,
QgsFeature,
QgsGeometry
)
from qgis.testing import start_app, unittest
from utilities import unitTestDataPath
start_app()
class TestQgsVectorLayerProfileGenerator(unittest.TestCase):
def testPointGenerationAbsolute(self):
vl = QgsVectorLayer(os.path.join(unitTestDataPath(), '3d', 'points_with_z.shp'), 'trees')
self.assertTrue(vl.isValid())
vl.elevationProperties().setClamping(Qgis.AltitudeClamping.Absolute)
vl.elevationProperties().setZScale(2.5)
vl.elevationProperties().setZOffset(10)
curve = QgsLineString()
curve.fromWkt(
'LineString (-347692.88994020794052631 6632796.97473032586276531, -346715.98066368483705446 6632458.84484416991472244, -346546.9152928851544857 6632444.29659010842442513, -346522.40419847227167338 6632307.79493184108287096)')
req = QgsProfileRequest(curve)
generator = vl.createProfileGenerator(req)
self.assertIsNotNone(generator)
# the request did not have the crs of the linestring set, so the whole linestring falls outside the vector extent
# TODO
# self.assertFalse(generator.generateProfile())
# set correct crs for linestring and re-try
req.setCrs(QgsCoordinateReferenceSystem('EPSG:3857'))
generator = vl.createProfileGenerator(req)
self.assertTrue(generator.generateProfile())
# no tolerance => no hits
results = generator.takeResults()
self.assertFalse(results.distanceToHeightMap())
req.setTolerance(20)
generator = vl.createProfileGenerator(req)
self.assertTrue(generator.generateProfile())
results = generator.takeResults()
self.assertEqual(results.distanceToHeightMap(),
{167.60402625236236: 274.0, 22.897466329720825: 274.0, 1159.1281061593656: 232.75000000000003,
1171.6873910420782: 235.5, 1197.5553515228498: 241.0, 1245.704759018164: 227.25,
1291.0771889584548: 221.75})
# lower tolerance
req.setTolerance(4.5)
generator = vl.createProfileGenerator(req)
self.assertTrue(generator.generateProfile())
results = generator.takeResults()
self.assertEqual(results.distanceToHeightMap(),
{167.60402625236236: 274.0, 1171.6873910420782: 235.5, 1197.5553515228498: 241.0})
def testPointGenerationTerrain(self):
"""
Points layer with terrain clamping
"""
vl = QgsVectorLayer(os.path.join(unitTestDataPath(), '3d', 'points_with_z.shp'), 'trees')
self.assertTrue(vl.isValid())
vl.elevationProperties().setClamping(Qgis.AltitudeClamping.Terrain)
vl.elevationProperties().setZScale(2.5)
vl.elevationProperties().setZOffset(10)
rl = QgsRasterLayer(os.path.join(unitTestDataPath(), '3d', 'dtm.tif'), 'DTM')
self.assertTrue(rl.isValid())
curve = QgsLineString()
curve.fromWkt(
'LineString (-347692.88994020794052631 6632796.97473032586276531, -346715.98066368483705446 6632458.84484416991472244, -346546.9152928851544857 6632444.29659010842442513, -346522.40419847227167338 6632307.79493184108287096)')
req = QgsProfileRequest(curve)
terrain_provider = QgsRasterDemTerrainProvider()
terrain_provider.setLayer(rl)
terrain_provider.setScale(0.3)
terrain_provider.setOffset(-5)
req.setTerrainProvider(terrain_provider)
req.setCrs(QgsCoordinateReferenceSystem('EPSG:3857'))
generator = vl.createProfileGenerator(req)
self.assertTrue(generator.generateProfile())
# no tolerance => no hits
results = generator.takeResults()
self.assertFalse(results.distanceToHeightMap())
req.setTolerance(4.5)
generator = vl.createProfileGenerator(req)
self.assertTrue(generator.generateProfile())
results = generator.takeResults()
self.assertEqual(results.distanceToHeightMap(),
{167.60402625236236: 69.5, 1171.6873910420782: 58.99999999999999, 1197.5553515228498: 60.5})
def testPointGenerationRelative(self):
"""
Points layer with relative clamping
"""
vl = QgsVectorLayer(os.path.join(unitTestDataPath(), '3d', 'points_with_z.shp'), 'trees')
self.assertTrue(vl.isValid())
vl.elevationProperties().setClamping(Qgis.AltitudeClamping.Relative)
vl.elevationProperties().setZScale(2.5)
vl.elevationProperties().setZOffset(10)
rl = QgsRasterLayer(os.path.join(unitTestDataPath(), '3d', 'dtm.tif'), 'DTM')
self.assertTrue(rl.isValid())
curve = QgsLineString()
curve.fromWkt(
'LineString (-347692.88994020794052631 6632796.97473032586276531, -346715.98066368483705446 6632458.84484416991472244, -346546.9152928851544857 6632444.29659010842442513, -346522.40419847227167338 6632307.79493184108287096)')
req = QgsProfileRequest(curve)
terrain_provider = QgsRasterDemTerrainProvider()
terrain_provider.setLayer(rl)
terrain_provider.setScale(0.3)
terrain_provider.setOffset(-5)
req.setTerrainProvider(terrain_provider)
req.setCrs(QgsCoordinateReferenceSystem('EPSG:3857'))
generator = vl.createProfileGenerator(req)
self.assertTrue(generator.generateProfile())
# no tolerance => no hits
results = generator.takeResults()
self.assertFalse(results.distanceToHeightMap())
req.setTolerance(4.5)
generator = vl.createProfileGenerator(req)
self.assertTrue(generator.generateProfile())
results = generator.takeResults()
self.assertEqual(results.distanceToHeightMap(),
{167.60402625236236: 333.5, 1171.6873910420782: 284.5, 1197.5553515228498: 291.5})
def testPointGenerationRelativeExtrusion(self):
"""
Points layer with relative clamping and extrusion
"""
vl = QgsVectorLayer(os.path.join(unitTestDataPath(), '3d', 'points_with_z.shp'), 'trees')
self.assertTrue(vl.isValid())
vl.elevationProperties().setClamping(Qgis.AltitudeClamping.Relative)
vl.elevationProperties().setZScale(2.5)
vl.elevationProperties().setZOffset(10)
vl.elevationProperties().setExtrusionEnabled(True)
vl.elevationProperties().setExtrusionHeight(7)
rl = QgsRasterLayer(os.path.join(unitTestDataPath(), '3d', 'dtm.tif'), 'DTM')
self.assertTrue(rl.isValid())
curve = QgsLineString()
curve.fromWkt(
'LineString (-347692.88994020794052631 6632796.97473032586276531, -346715.98066368483705446 6632458.84484416991472244, -346546.9152928851544857 6632444.29659010842442513, -346522.40419847227167338 6632307.79493184108287096)')
req = QgsProfileRequest(curve)
terrain_provider = QgsRasterDemTerrainProvider()
terrain_provider.setLayer(rl)
terrain_provider.setScale(0.3)
terrain_provider.setOffset(-5)
req.setTerrainProvider(terrain_provider)
req.setCrs(QgsCoordinateReferenceSystem('EPSG:3857'))
generator = vl.createProfileGenerator(req)
self.assertTrue(generator.generateProfile())
req.setTolerance(4.5)
generator = vl.createProfileGenerator(req)
self.assertTrue(generator.generateProfile())
results = generator.takeResults()
self.assertEqual(results.distanceToHeightMap(),
{167.60402625236236: 333.5, 1171.6873910420782: 284.5, 1197.5553515228498: 291.5})
self.assertCountEqual([g.asWkt(1) for g in results.asGeometries()],
['LineStringZ (-347535.1 6632740.6 333.5, -347535.1 6632740.6 340.5)',
'LineStringZ (-346578.2 6632450.9 284.5, -346578.2 6632450.9 291.5)',
'LineStringZ (-346552.5 6632448.8 291.5, -346552.5 6632448.8 298.5)'])
def testPointGenerationMultiPoint(self):
vl = QgsVectorLayer('MultipointZ?crs=EPSG:27700', 'trees', 'memory')
self.assertTrue(vl.isValid())
f = QgsFeature()
f.setGeometry(QgsGeometry.fromWkt('MultiPointZ(322069 129893 89.1, 322077 129889 90.2, 322093 129888 92.4)'))
self.assertTrue(vl.dataProvider().addFeature(f))
vl.elevationProperties().setClamping(Qgis.AltitudeClamping.Absolute)
vl.elevationProperties().setZScale(2.5)
vl.elevationProperties().setZOffset(10)
curve = QgsLineString()
curve.fromWkt(
'LineString (-347692.88994020794052631 6632796.97473032586276531, -346715.98066368483705446 6632458.84484416991472244, -346546.9152928851544857 6632444.29659010842442513, -346522.40419847227167338 6632307.79493184108287096)')
req = QgsProfileRequest(curve)
req.setCrs(QgsCoordinateReferenceSystem('EPSG:3857'))
req.setTolerance(20)
generator = vl.createProfileGenerator(req)
self.assertTrue(generator.generateProfile())
results = generator.takeResults()
self.assertEqual(results.distanceToHeightMap(),
{1158.036781085277: 232.75, 1171.3215208424097: 235.5, 1196.7677444714948: 241.0})
if __name__ == '__main__':
unittest.main()