QGIS API Documentation 3.39.0-Master (67e056379ed)
Loading...
Searching...
No Matches
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 >() << static_cast< int >( Qgis::ProcessingSourceType::VectorPoint ) ) );
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" ), Qgis::ProcessingSourceType::VectorPoint ) );
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 ? QMetaType::Type::Int : QMetaType::Type::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
@ VectorPoint
Vector point layers.
DataType
Raster data types.
Definition qgis.h:288
@ 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.
Class for doing transforms between two map coordinate systems.
Defines a QGIS exception class.
Wrapper for iterator of features from vector data provider or vector layer.
bool nextFeature(QgsFeature &f)
Fetch next feature and stores in f, returns true on success.
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:58
QgsAttributes attributes
Definition qgsfeature.h:67
QgsFeatureId id
Definition qgsfeature.h:66
QgsGeometry geometry
Definition qgsfeature.h:69
bool hasGeometry() const
Returns true if the feature has an associated geometry.
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition qgsfeedback.h:53
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition qgsfeedback.h:61
Encapsulate a field in an attribute table or data source.
Definition qgsfield.h:53
Container of fields for a vector layer.
Definition qgsfields.h:46
bool append(const QgsField &field, Qgis::FieldOrigin origin=Qgis::FieldOrigin::Provider, int originIndex=-1)
Appends a field.
Definition qgsfields.cpp:69
A geometry is the spatial representation of a feature.
QgsAbstractGeometry * get()
Returns a modifiable (non-const) reference to the underlying abstract geometry primitive.
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
bool isMultipart() const
Returns true if WKB of the geometry is of WKBMulti* type.
QgsCoordinateReferenceSystem crs
Definition qgsmaplayer.h:82
A class to represent a 2D point.
Definition qgspointxy.h:60
Point geometry type, with support for z-dimension and m-values.
Definition qgspoint.h:49
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.
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).
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.
T qgsgeometry_cast(const QgsAbstractGeometry *geom)