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
18 #include "qgsalgorithmrastersampling.h"
19 #include "qgsmultipoint.h"
20
21 ///@cond PRIVATE
22
name() const23 QString QgsRasterSamplingAlgorithm::name() const
24 {
25 return QStringLiteral( "rastersampling" );
26 }
27
displayName() const28 QString QgsRasterSamplingAlgorithm::displayName() const
29 {
30 return QObject::tr( "Sample raster values" );
31 }
32
tags() const33 QStringList QgsRasterSamplingAlgorithm::tags() const
34 {
35 return QObject::tr( "extract,point,pixel,value" ).split( ',' );
36 }
37
group() const38 QString QgsRasterSamplingAlgorithm::group() const
39 {
40 return QObject::tr( "Raster analysis" );
41 }
42
groupId() const43 QString QgsRasterSamplingAlgorithm::groupId() const
44 {
45 return QStringLiteral( "rasteranalysis" );
46 }
47
shortDescription() const48 QString QgsRasterSamplingAlgorithm::shortDescription() const
49 {
50 return QObject::tr( "Samples raster values under a set of points." );
51 }
52
shortHelpString() const53 QString 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
createInstance() const59 QgsRasterSamplingAlgorithm *QgsRasterSamplingAlgorithm::createInstance() const
60 {
61 return new QgsRasterSamplingAlgorithm();
62 }
63
initAlgorithm(const QVariantMap &)64 void 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
prepareAlgorithm(const QVariantMap & parameters,QgsProcessingContext & context,QgsProcessingFeedback *)73 bool QgsRasterSamplingAlgorithm::prepareAlgorithm( const QVariantMap ¶meters, 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
processAlgorithm(const QVariantMap & parameters,QgsProcessingContext & context,QgsProcessingFeedback * feedback)86 QVariantMap QgsRasterSamplingAlgorithm::processAlgorithm( const QVariantMap ¶meters, 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
181 ///@endcond
182