QGIS API Documentation 3.28.0-Firenze (ed3ad0430f)
qgsalgorithmrastersampling.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsalgorithmrastersampling.cpp
3 --------------------------
4 begin : August 2020
5 copyright : (C) 2020 by Mathieu Pellerin
6 email : nirvn dot asia at gmail dot com
7 ***************************************************************************/
8
9/***************************************************************************
10 * *
11 * This program is free software; you can redistribute it and/or modify *
12 * it under the terms of the GNU General Public License as published by *
13 * the Free Software Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 ***************************************************************************/
17
19#include "qgsmultipoint.h"
20
22
23QString QgsRasterSamplingAlgorithm::name() const
24{
25 return QStringLiteral( "rastersampling" );
26}
27
28QString QgsRasterSamplingAlgorithm::displayName() const
29{
30 return QObject::tr( "Sample raster values" );
31}
32
33QStringList QgsRasterSamplingAlgorithm::tags() const
34{
35 return QObject::tr( "extract,point,pixel,value" ).split( ',' );
36}
37
38QString QgsRasterSamplingAlgorithm::group() const
39{
40 return QObject::tr( "Raster analysis" );
41}
42
43QString QgsRasterSamplingAlgorithm::groupId() const
44{
45 return QStringLiteral( "rasteranalysis" );
46}
47
48QString QgsRasterSamplingAlgorithm::shortDescription() const
49{
50 return QObject::tr( "Samples raster values under a set of points." );
51}
52
53QString QgsRasterSamplingAlgorithm::shortHelpString() const
54{
55 return QObject::tr( "This algorithm creates a new vector layer with the same attributes of the input layer and the raster values corresponding on the point location.\n\n"
56 "If the raster layer has more than one band, all the band values are sampled." );
57}
58
59QgsRasterSamplingAlgorithm *QgsRasterSamplingAlgorithm::createInstance() const
60{
61 return new QgsRasterSamplingAlgorithm();
62}
63
64void QgsRasterSamplingAlgorithm::initAlgorithm( const QVariantMap & )
65{
66 addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "INPUT" ), QObject::tr( "Input layer" ), QList< int >() << QgsProcessing::TypeVectorPoint ) );
67 addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "RASTERCOPY" ), QObject::tr( "Raster layer" ) ) );
68
69 addParameter( new QgsProcessingParameterString( QStringLiteral( "COLUMN_PREFIX" ), QObject::tr( "Output column prefix" ), QStringLiteral( "SAMPLE_" ), false, true ) );
70 addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT" ), QObject::tr( "Sampled" ), QgsProcessing::TypeVectorPoint ) );
71}
72
73bool QgsRasterSamplingAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback * )
74{
75 QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral( "RASTERCOPY" ), context );
76 if ( !layer )
77 throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "RASTERCOPY" ) ) );
78
79 mBandCount = layer->bandCount();
80 mCrs = layer->crs();
81 mDataProvider.reset( static_cast< QgsRasterDataProvider *>( layer->dataProvider()->clone() ) );
82
83 return true;
84}
85
86QVariantMap QgsRasterSamplingAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
87{
88 std::unique_ptr< QgsFeatureSource > source( parameterAsSource( parameters, QStringLiteral( "INPUT" ), context ) );
89 if ( !source )
90 throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "INPUT" ) ) );
91
92 const QString fieldPrefix = parameterAsString( parameters, QStringLiteral( "COLUMN_PREFIX" ), context );
93 QgsFields newFields;
94 QgsAttributes emptySampleAttributes;
95 for ( int band = 1; band <= mBandCount; band++ )
96 {
97 const Qgis::DataType dataType = mDataProvider->dataType( band );
98 const bool intSafe = ( dataType == Qgis::DataType::Byte || dataType == Qgis::DataType::UInt16 || dataType == Qgis::DataType::Int16 || dataType == Qgis::DataType::UInt32 ||
99 dataType == Qgis::DataType::Int32 || dataType == Qgis::DataType::CInt16 || dataType == Qgis::DataType::CInt32 );
100
101 newFields.append( QgsField( QStringLiteral( "%1%2" ).arg( fieldPrefix, QString::number( band ) ), intSafe ? QVariant::Int : QVariant::Double ) );
102 emptySampleAttributes += QVariant();
103 }
104 const QgsFields fields = QgsProcessingUtils::combineFields( source->fields(), newFields );
105
106 QString dest;
107 std::unique_ptr< QgsFeatureSink > sink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, dest, fields,
108 source->wkbType(), source->sourceCrs() ) );
109 if ( !sink )
110 throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT" ) ) );
111
112 const long count = source->featureCount();
113 const double step = count > 0 ? 100.0 / count : 1;
114 long current = 0;
115
116 const QgsCoordinateTransform ct( source->sourceCrs(), mCrs, context.transformContext() );
117 QgsFeatureIterator it = source->getFeatures( QgsFeatureRequest() );
118 QgsFeature feature;
119 while ( it.nextFeature( feature ) )
120 {
121 if ( feedback->isCanceled() )
122 {
123 break;
124 }
125 feedback->setProgress( current * step );
126 current++;
127
128 QgsAttributes attributes = feature.attributes();
129 QgsFeature outputFeature( feature );
130 if ( !feature.hasGeometry() )
131 {
132 attributes += emptySampleAttributes;
133 outputFeature.setAttributes( attributes );
134 sink->addFeature( outputFeature, QgsFeatureSink::FastInsert );
135 feedback->reportError( QObject::tr( "No geometry attached to feature %1." ).arg( feature.id() ) );
136 continue;
137 }
138
139 QgsGeometry geometry = feature.geometry();
140 if ( geometry.isMultipart() && geometry.get()->partCount() != 1 )
141 {
142 attributes += emptySampleAttributes;
143 outputFeature.setAttributes( attributes );
144 if ( !sink->addFeature( outputFeature, QgsFeatureSink::FastInsert ) )
145 throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QStringLiteral( "OUTPUT" ) ) );
146 feedback->reportError( QObject::tr( "Impossible to sample data of multipart feature %1." ).arg( feature.id() ) );
147 continue;
148 }
149 QgsPointXY point( *( geometry.isMultipart() ? qgsgeometry_cast< const QgsPoint * >( qgsgeometry_cast< const QgsMultiPoint * >( geometry.constGet() )->geometryN( 0 ) ) :
150 qgsgeometry_cast< const QgsPoint * >( geometry.constGet() ) ) );
151 try
152 {
153 point = ct.transform( point );
154 }
155 catch ( const QgsException & )
156 {
157 attributes += emptySampleAttributes;
158 outputFeature.setAttributes( attributes );
159 if ( !sink->addFeature( outputFeature, QgsFeatureSink::FastInsert ) )
160 throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QStringLiteral( "OUTPUT" ) ) );
161 feedback->reportError( QObject::tr( "Could not reproject feature %1 to raster CRS." ).arg( feature.id() ) );
162 continue;
163 }
164
165 for ( int band = 1; band <= mBandCount; band ++ )
166 {
167 bool ok = false;
168 const double value = mDataProvider->sample( point, band, &ok );
169 attributes += ok ? value : QVariant();
170 }
171 outputFeature.setAttributes( attributes );
172 if ( !sink->addFeature( outputFeature, QgsFeatureSink::FastInsert ) )
173 throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QStringLiteral( "OUTPUT" ) ) );
174 }
175
176 QVariantMap outputs;
177 outputs.insert( QStringLiteral( "OUTPUT" ), dest );
178 return outputs;
179}
180
DataType
Raster data types.
Definition: qgis.h:129
@ CInt32
Complex Int32.
@ Int16
Sixteen bit signed integer (qint16)
@ UInt16
Sixteen bit unsigned integer (quint16)
@ Byte
Eight bit unsigned integer (quint8)
@ Int32
Thirty two bit signed integer (qint32)
@ CInt16
Complex Int16.
@ UInt32
Thirty two bit unsigned integer (quint32)
virtual int partCount() const =0
Returns count of parts contained in the geometry.
A vector of attributes.
Definition: qgsattributes.h:59
Class for doing transforms between two map coordinate systems.
Defines a QGIS exception class.
Definition: qgsexception.h:35
Wrapper for iterator of features from vector data provider or vector layer.
bool nextFeature(QgsFeature &f)
This class wraps a request for features to a vector layer (or directly its vector data provider).
@ FastInsert
Use faster inserts, at the cost of updating the passed features to reflect changes made at the provid...
The feature class encapsulates a single feature including its unique ID, geometry and a list of field...
Definition: qgsfeature.h:56
QgsAttributes attributes
Definition: qgsfeature.h:65
QgsGeometry geometry
Definition: qgsfeature.h:67
bool hasGeometry() const
Returns true if the feature has an associated geometry.
Definition: qgsfeature.cpp:233
Q_GADGET QgsFeatureId id
Definition: qgsfeature.h:64
bool isCanceled() const SIP_HOLDGIL
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:54
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:63
Encapsulate a field in an attribute table or data source.
Definition: qgsfield.h:51
Container of fields for a vector layer.
Definition: qgsfields.h:45
bool append(const QgsField &field, FieldOrigin origin=OriginProvider, int originIndex=-1)
Appends a field. The field must have unique name, otherwise it is rejected (returns false)
Definition: qgsfields.cpp:59
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:164
const QgsAbstractGeometry * constGet() const SIP_HOLDGIL
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
QgsAbstractGeometry * get()
Returns a modifiable (non-const) reference to the underlying abstract geometry primitive.
bool isMultipart() const SIP_HOLDGIL
Returns true if WKB of the geometry is of WKBMulti* type.
QgsCoordinateReferenceSystem crs
Definition: qgsmaplayer.h:79
A class to represent a 2D point.
Definition: qgspointxy.h:59
Contains information about the context in which a processing algorithm is executed.
QgsCoordinateTransformContext transformContext() const
Returns the coordinate transform context.
Custom exception class for processing related exceptions.
Definition: qgsexception.h:83
Base class for providing feedback from a processing algorithm.
virtual void reportError(const QString &error, bool fatalError=false)
Reports that the algorithm encountered an error while executing.
A feature sink output for processing algorithms.
An input feature source (such as vector layers) parameter for processing algorithms.
A raster layer parameter for processing algorithms.
A string parameter for processing algorithms.
static QgsFields combineFields(const QgsFields &fieldsA, const QgsFields &fieldsB, const QString &fieldsBPrefix=QString())
Combines two field lists, avoiding duplicate field names (in a case-insensitive manner).
@ TypeVectorPoint
Vector point layers.
Definition: qgsprocessing.h:49
Base class for raster data providers.
QgsRasterDataProvider * clone() const override=0
Clone itself, create deep copy.
Represents a raster layer.
int bandCount() const
Returns the number of bands in this layer.
QgsRasterDataProvider * dataProvider() override
Returns the source data provider.