Branch data Line data Source code
1 : : /***************************************************************************
2 : : qgsalgorithmlinedensity.cpp
3 : : ---------------------
4 : : begin : December 2019
5 : : copyright : (C) 2019 by Clemens Raffler
6 : : email : clemens dot raffler 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 "qgsalgorithmlinedensity.h"
19 : : #include "qgscircle.h"
20 : : #include "qgsgeometryengine.h"
21 : : #include "qgsrasterfilewriter.h"
22 : :
23 : : ///@cond PRIVATE
24 : :
25 : 0 : QString QgsLineDensityAlgorithm::name() const
26 : : {
27 : 0 : return QStringLiteral( "linedensity" );
28 : : }
29 : :
30 : 0 : QString QgsLineDensityAlgorithm::displayName() const
31 : : {
32 : 0 : return QObject::tr( "Line density" );
33 : : }
34 : :
35 : 0 : QStringList QgsLineDensityAlgorithm::tags() const
36 : : {
37 : 0 : return QObject::tr( "density,kernel,line,line density,interpolation,weight" ).split( ',' );
38 : 0 : }
39 : :
40 : 0 : QString QgsLineDensityAlgorithm::group() const
41 : : {
42 : 0 : return QObject::tr( "Interpolation" );
43 : : }
44 : :
45 : 0 : QString QgsLineDensityAlgorithm::groupId() const
46 : : {
47 : 0 : return QStringLiteral( "interpolation" );
48 : : }
49 : :
50 : 0 : void QgsLineDensityAlgorithm::initAlgorithm( const QVariantMap & )
51 : : {
52 : 0 : addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "INPUT" ), QObject::tr( "Input line layer" ), QList<int>() << QgsProcessing::TypeVectorLine ) );
53 : 0 : addParameter( new QgsProcessingParameterField( QStringLiteral( "WEIGHT" ), QObject::tr( "Weight field " ), QVariant(), QStringLiteral( "INPUT" ), QgsProcessingParameterField::Numeric, false, true ) );
54 : 0 : addParameter( new QgsProcessingParameterDistance( QStringLiteral( "RADIUS" ), QObject::tr( "Search radius" ), 10, QStringLiteral( "INPUT" ), false, 0 ) );
55 : 0 : addParameter( new QgsProcessingParameterDistance( QStringLiteral( "PIXEL_SIZE" ), QObject::tr( "Pixel size" ), 10, QStringLiteral( "INPUT" ), false ) );
56 : :
57 : 0 : addParameter( new QgsProcessingParameterRasterDestination( QStringLiteral( "OUTPUT" ), QObject::tr( "Line density raster" ) ) );
58 : 0 : }
59 : :
60 : 0 : QString QgsLineDensityAlgorithm::shortHelpString() const
61 : : {
62 : 0 : return QObject::tr( "The line density interpolation algorithm calculates a density measure of linear features "
63 : : "which is obtained in a circular neighborhood within each raster cell. "
64 : : "First, the length of the segment of each line that is intersected by the circular neighborhood "
65 : : "is multiplied with the lines weight factor. In a second step, all length values are summed and "
66 : : "divided by the area of the circular neighborhood. This process is repeated for all raster cells."
67 : : );
68 : : }
69 : :
70 : 0 : QgsLineDensityAlgorithm *QgsLineDensityAlgorithm::createInstance() const
71 : : {
72 : 0 : return new QgsLineDensityAlgorithm();
73 : 0 : }
74 : :
75 : 0 : bool QgsLineDensityAlgorithm::prepareAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
76 : : {
77 : : Q_UNUSED( feedback );
78 : 0 : mSource.reset( parameterAsSource( parameters, QStringLiteral( "INPUT" ), context ) );
79 : 0 : if ( !mSource )
80 : 0 : throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "INPUT" ) ) );
81 : :
82 : 0 : mWeightField = parameterAsString( parameters, QStringLiteral( "WEIGHT" ), context );
83 : :
84 : 0 : mPixelSize = parameterAsDouble( parameters, QStringLiteral( "PIXEL_SIZE" ), context );
85 : :
86 : 0 : mSearchRadius = parameterAsDouble( parameters, QStringLiteral( "RADIUS" ), context );
87 : 0 : if ( mSearchRadius < 0.5 * mPixelSize * std::sqrt( 2 ) )
88 : 0 : throw QgsProcessingException( QObject::tr( "Raster cells must be fully contained by the search circle. Therefore, "
89 : : "the search radius must not be smaller than half of the pixel diagonal." ) );
90 : :
91 : 0 : mExtent = mSource->sourceExtent();
92 : 0 : mCrs = mSource->sourceCrs();
93 : 0 : mDa = QgsDistanceArea();
94 : 0 : mDa.setEllipsoid( context.ellipsoid() );
95 : 0 : mDa.setSourceCrs( mCrs, context.transformContext() );
96 : :
97 : : //get cell midpoint from top left cell
98 : 0 : QgsPoint firstCellMidpoint = QgsPoint( mExtent.xMinimum() + ( mPixelSize / 2 ), mExtent.yMaximum() - ( mPixelSize / 2 ) );
99 : 0 : QgsCircle searchCircle = QgsCircle( firstCellMidpoint, mSearchRadius );
100 : 0 : mSearchGeometry = QgsGeometry( searchCircle.toPolygon() );
101 : :
102 : : return true;
103 : 0 : }
104 : :
105 : 0 : QVariantMap QgsLineDensityAlgorithm::processAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
106 : : {
107 : 0 : mIndex = QgsSpatialIndex( QgsSpatialIndex::FlagStoreFeatureGeometries );
108 : :
109 : 0 : QStringList weightName = QStringList( mWeightField );
110 : 0 : QgsFields attrFields = mSource->fields();
111 : :
112 : 0 : QgsFeatureRequest r = QgsFeatureRequest();
113 : 0 : r.setSubsetOfAttributes( weightName, attrFields );
114 : 0 : QgsFeatureIterator fit = mSource->getFeatures( r );
115 : 0 : QgsFeature f;
116 : :
117 : 0 : while ( fit.nextFeature( f ) )
118 : : {
119 : 0 : mIndex.addFeature( f, QgsFeatureSink::FastInsert );
120 : :
121 : : //only populate hash if weight field is given
122 : 0 : if ( !mWeightField.isEmpty() )
123 : : {
124 : 0 : double analysisWeight = f.attribute( mWeightField ).toDouble();
125 : 0 : mFeatureWeights.insert( f.id(), analysisWeight );
126 : 0 : }
127 : : }
128 : :
129 : 0 : const QString outputFile = parameterAsOutputLayer( parameters, QStringLiteral( "OUTPUT" ), context );
130 : 0 : QFileInfo fi( outputFile );
131 : 0 : const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() );
132 : :
133 : 0 : int rows = std::max( std::ceil( mExtent.height() / mPixelSize ), 1.0 );
134 : 0 : int cols = std::max( std::ceil( mExtent.width() / mPixelSize ), 1.0 );
135 : :
136 : : //build new raster extent based on number of columns and cellsize
137 : : //this prevents output cellsize being calculated too small
138 : 0 : QgsRectangle rasterExtent = QgsRectangle( mExtent.xMinimum(), mExtent.yMaximum() - ( rows * mPixelSize ), mExtent.xMinimum() + ( cols * mPixelSize ), mExtent.yMaximum() );
139 : :
140 : 0 : QgsRasterFileWriter writer = QgsRasterFileWriter( outputFile );
141 : 0 : writer.setOutputProviderKey( QStringLiteral( "gdal" ) );
142 : 0 : writer.setOutputFormat( outputFormat );
143 : 0 : std::unique_ptr<QgsRasterDataProvider > provider( writer.createOneBandRaster( Qgis::Float32, cols, rows, rasterExtent, mCrs ) );
144 : 0 : if ( !provider )
145 : 0 : throw QgsProcessingException( QObject::tr( "Could not create raster output: %1" ).arg( outputFile ) );
146 : 0 : if ( !provider->isValid() )
147 : 0 : throw QgsProcessingException( QObject::tr( "Could not create raster output %1: %2" ).arg( outputFile, provider->error().message( QgsErrorMessage::Text ) ) );
148 : :
149 : 0 : provider->setNoDataValue( 1, -9999 );
150 : :
151 : 0 : qgssize totalCellcnt = static_cast<qgssize>( rows ) * cols;
152 : 0 : int cellcnt = 0;
153 : :
154 : 0 : std::unique_ptr< QgsRasterBlock > rasterDataLine = std::make_unique< QgsRasterBlock >( Qgis::Float32, cols, 1 );
155 : :
156 : 0 : for ( int row = 0; row < rows; row++ )
157 : : {
158 : 0 : for ( int col = 0; col < cols; col++ )
159 : : {
160 : 0 : if ( feedback->isCanceled() )
161 : : {
162 : 0 : break;
163 : : }
164 : :
165 : 0 : if ( col > 0 )
166 : 0 : mSearchGeometry.translate( mPixelSize, 0 );
167 : :
168 : 0 : const QList<QgsFeatureId> fids = mIndex.intersects( mSearchGeometry.boundingBox() );
169 : :
170 : 0 : if ( !fids.isEmpty() )
171 : : {
172 : 0 : std::unique_ptr< QgsGeometryEngine > engine( QgsGeometry::createGeometryEngine( mSearchGeometry.constGet() ) );
173 : 0 : engine->prepareGeometry();
174 : :
175 : 0 : double absDensity = 0;
176 : 0 : for ( QgsFeatureId id : fids )
177 : : {
178 : 0 : QgsGeometry lineGeom = mIndex.geometry( id );
179 : :
180 : 0 : if ( engine->intersects( lineGeom.constGet() ) )
181 : : {
182 : 0 : double analysisLineLength = mDa.measureLength( QgsGeometry( engine->intersection( mIndex.geometry( id ).constGet() ) ) );
183 : 0 : double weight = 1;
184 : :
185 : 0 : if ( !mWeightField.isEmpty() )
186 : : {
187 : 0 : weight = mFeatureWeights.value( id );
188 : 0 : }
189 : :
190 : 0 : absDensity += ( analysisLineLength * weight );
191 : 0 : }
192 : 0 : }
193 : :
194 : 0 : double lineDensity = 0;
195 : 0 : if ( absDensity > 0 )
196 : : {
197 : : //only calculate ellipsoidal area if abs density is greater 0
198 : 0 : double analysisSearchGeometryArea = mDa.measureArea( mSearchGeometry );
199 : 0 : lineDensity = absDensity / analysisSearchGeometryArea;
200 : 0 : }
201 : 0 : rasterDataLine->setValue( 0, col, lineDensity );
202 : 0 : }
203 : : else
204 : : {
205 : : //no lines found in search radius
206 : 0 : rasterDataLine->setValue( 0, col, 0.0 );
207 : : }
208 : :
209 : 0 : feedback->setProgress( static_cast<double>( cellcnt ) / static_cast<double>( totalCellcnt ) * 100 );
210 : 0 : cellcnt++;
211 : 0 : }
212 : 0 : provider->writeBlock( rasterDataLine.get(), 1, 0, row );
213 : :
214 : : //'carriage return and newline' for search geometry
215 : 0 : mSearchGeometry.translate( ( cols - 1 ) * -mPixelSize, -mPixelSize );
216 : 0 : }
217 : :
218 : 0 : QVariantMap outputs;
219 : 0 : outputs.insert( QStringLiteral( "OUTPUT" ), outputFile );
220 : 0 : return outputs;
221 : 0 : }
222 : :
223 : :
224 : : ///@endcond
|
