Skip to content
Permalink
Browse files

[feature][processing] Add Line denisty algorithm

  • Loading branch information
root676 authored and nyalldawson committed Dec 26, 2019
1 parent 5765cab commit ce342f7ac4ebec772c5e8142e1549634de34ed5f
@@ -102,6 +102,7 @@
<file>themes/default/algorithms/mAlgorithmExtractLayerExtent.svg</file>
<file>themes/default/algorithms/mAlgorithmIntersect.svg</file>
<file>themes/default/algorithms/mAlgorithmConcaveHull.svg</file>
<file>themes/default/algorithms/mAlgorithmLineDensity.svg</file>
<file>themes/default/algorithms/mAlgorithmLineIntersections.svg</file>
<file>themes/default/algorithms/mAlgorithmLineToPolygon.svg</file>
<file>themes/default/algorithms/mAlgorithmMeanCoordinates.svg</file>
@@ -0,0 +1 @@
<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 6.35 6.35" height="24" width="24"><path d="M.194 2.682L5.796.537M.968 5.814s.461-.596.825-.893c.364-.296.75-.699 1.097-.872.17-.085.51-.137.958.044.449.182 1.02-.558 1.02-.558l.282-.35M3.102 2.274C4.88 1.019 2.153.93 1.642.219" fill="none" stroke="#517083" stroke-width=".529" stroke-linecap="round"/><path d="M3.896 4.098s-.499.708-.338 1.021c.16.314.627 1.021.627 1.021" fill="none" stroke="#517083" stroke-width=".529" stroke-linecap="round"/><path fill="none" stroke="#1a1a1a" stroke-width=".244" stroke-linecap="round" d="M2.863 2.725h.823v.831h-.823z"/><path d="M3.896 4.098s-.499.708-.338 1.021c.08.146.153.254.153.254M3.102 2.274C4.534 1.26 3.108 1.04 2.886.896" fill="none" stroke="#92a9c2" stroke-width=".529" stroke-linecap="round"/><path d="M1.807 4.921c.347-.289.736-.699 1.083-.872.17-.085.51-.137.958.044.449.182 1.02-.558 1.02-.558l.282-.35M1.108 2.337l3.144-1.212" fill="none" stroke="#92a9c2" stroke-width=".529" stroke-linecap="round"/><circle cx="3.206" cy="3.145" r="2.485" fill="none" stroke="#1a1a1a" stroke-width=".212" stroke-linecap="round" stroke-linejoin="round"/></svg>
@@ -0,0 +1,50 @@
<?xml version="1.0" encoding="utf-8" ?>
<ogr:FeatureCollection
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:schemaLocation="http://ogr.maptools.org/ linedensity.xsd"
xmlns:ogr="http://ogr.maptools.org/"
xmlns:gml="http://www.opengis.net/gml">
<gml:boundedBy>
<gml:Box>
<gml:coord><gml:X>0.3605633802816897</gml:X><gml:Y>-3.816901408450707</gml:Y></gml:coord>
<gml:coord><gml:X>9.930985915492959</gml:X><gml:Y>5.512676056338028</gml:Y></gml:coord>
</gml:Box>
</gml:boundedBy>

<gml:featureMember>
<ogr:linedensity fid="1">
<ogr:geometryProperty><gml:LineString srsName="EPSG:4326"><gml:coordinates>0.753521126760563,2.11549295774648 1.65352112676056,4.04225352112676 3.1112676056338,3.96619718309859 5.19014084507042,5.37323943661972 6.66056338028169,4.86619718309859 7.7887323943662,4.94225352112676 8.70140845070423,5.51267605633803 9.5887323943662,5.03098591549296</gml:coordinates></gml:LineString></ogr:geometryProperty>
<ogr:weight>1.000</ogr:weight>
</ogr:linedensity>
</gml:featureMember>
<gml:featureMember>
<ogr:linedensity fid="2">
<ogr:geometryProperty><gml:LineString srsName="EPSG:4326"><gml:coordinates>3.66901408450704,4.16901408450704 4.55633802816901,3.31971830985915 3.6056338028169,2.45774647887324 4.78450704225352,1.78591549295775</gml:coordinates></gml:LineString></ogr:geometryProperty>
<ogr:weight>1.000</ogr:weight>
</ogr:linedensity>
</gml:featureMember>
<gml:featureMember>
<ogr:linedensity fid="3">
<ogr:geometryProperty><gml:LineString srsName="EPSG:4326"><gml:coordinates>3.45352112676056,-1.01549295774648 5.67183098591549,0.239436619718309 6.05211267605634,1.59577464788732 6.58450704225352,3.49718309859155 7.62394366197183,3.40845070422535 8.67605633802817,4.71408450704225 9.22112676056338,4.37183098591549 9.93098591549296,4.15633802816901</gml:coordinates></gml:LineString></ogr:geometryProperty>
<ogr:weight>2.000</ogr:weight>
</ogr:linedensity>
</gml:featureMember>
<gml:featureMember>
<ogr:linedensity fid="4">
<ogr:geometryProperty><gml:LineString srsName="EPSG:4326"><gml:coordinates>1.20985915492958,1.03802816901408 3.45352112676056,0.797183098591549</gml:coordinates></gml:LineString></ogr:geometryProperty>
<ogr:weight>1.000</ogr:weight>
</ogr:linedensity>
</gml:featureMember>
<gml:featureMember>
<ogr:linedensity fid="5">
<ogr:geometryProperty><gml:LineString srsName="EPSG:4326"><gml:coordinates>0.36056338028169,4.38450704225352 1.60281690140845,2.36901408450704 2.05915492957746,2.49577464788732 2.05915492957746,3.59859154929577 2.26197183098592,3.94084507042253 2.83239436619718,2.72394366197183 2.62957746478873,4.57464788732394 2.8830985915493,5.29718309859155 3.5169014084507,4.65070422535211 3.75774647887324,3.91549295774648 3.37746478873239,3.58591549295775 3.04788732394366,3.09154929577465 3.16197183098592,2.48309859154929</gml:coordinates></gml:LineString></ogr:geometryProperty>
<ogr:weight>3.000</ogr:weight>
</ogr:linedensity>
</gml:featureMember>
<gml:featureMember>
<ogr:linedensity fid="6">
<ogr:geometryProperty><gml:LineString srsName="EPSG:4326"><gml:coordinates>1.07042253521127,-2.52394366197183 3.33943661971831,-3.81690140845071 6.34366197183099,-2.09295774647887 8.29577464788733,-2.14366197183099 9.32253521126761,-0.736619718309861</gml:coordinates></gml:LineString></ogr:geometryProperty>
<ogr:weight>1.000</ogr:weight>
</ogr:linedensity>
</gml:featureMember>
</ogr:FeatureCollection>
@@ -982,6 +982,20 @@ tests:
bottom:
precision: 7

- algorithm: native:linedensity
name: Line Density
params:
INPUT:
name: linedensity.gml|layername=linedensity
type: vector
PIXEL_SIZE: 1.0
RADIUS: 2.0
WEIGHT: weight
results:
OUTPUT:
hash: 37b153a98004a6656a27a8c5dcdb8dbefcd4af4938a512ec0cf5cd43
type: rasterhash

- algorithm: native:deleteholes
name: Delete holes (no min)
params:
@@ -79,6 +79,7 @@ SET(QGIS_ANALYSIS_SRCS
processing/qgsalgorithmjoinbynearest.cpp
processing/qgsalgorithmjoinwithlines.cpp
processing/qgsalgorithmkmeansclustering.cpp
processing/qgsalgorithmlinedensity.cpp
processing/qgsalgorithmlineintersection.cpp
processing/qgsalgorithmlinesubstring.cpp
processing/qgsalgorithmloadlayer.cpp
@@ -0,0 +1,193 @@
/***************************************************************************
qgsalgorithmlinedensity.cpp
---------------------
begin : December 2019
copyright : (C) 2019 by Clemens Raffler
email : clemens dot raffler 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 "qgsalgorithmlinedensity.h"
#include "qgsgeometryengine.h"
#include "qgsrasterfilewriter.h"

///@cond PRIVATE

QString QgsLineDensityAlgorithm::name() const
{
return QStringLiteral( "linedensity" );
}

QString QgsLineDensityAlgorithm::displayName() const
{
return QObject::tr( "Line density" );
}

QStringList QgsLineDensityAlgorithm::tags() const
{
return QObject::tr( "density,kernel,line,line density,interpolation,weight" ).split( ',' );
}

QString QgsLineDensityAlgorithm::group() const
{
return QObject::tr( "Interpolation" );
}

QString QgsLineDensityAlgorithm::groupId() const
{
return QStringLiteral( "interpolation" );
}

void QgsLineDensityAlgorithm::initAlgorithm( const QVariantMap & )
{
addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "INPUT" ), QObject::tr( "Input line layer" ), QList<int>() << QgsProcessing::TypeVectorLine ) );
addParameter( new QgsProcessingParameterField( QStringLiteral( "WEIGHT" ), QObject::tr( "Weight field "), QVariant(), QStringLiteral( "INPUT" ), QgsProcessingParameterField::Numeric ) );
addParameter( new QgsProcessingParameterDistance( QStringLiteral( "RADIUS" ), QObject::tr( "Search radius" ), 10, QStringLiteral( "INPUT" ), false, 0) );
addParameter( new QgsProcessingParameterDistance( QStringLiteral( "PIXEL_SIZE" ), QObject::tr( "Pixel size" ), 10, QStringLiteral( "INPUT" ), false) );

addParameter( new QgsProcessingParameterRasterDestination( QStringLiteral( "OUTPUT" ), QObject::tr( "Line density raster" ) ) );
}

QString QgsLineDensityAlgorithm::shortHelpString() const
{
return QObject::tr( "The line density interpolation algorithm calculates a density measure of linear features "
"which is obtained in a circular neighborhood within each raster cell. "
"First, the length of the segment of each line that is intersected by the circular neighborhood "
"is multiplied with the lines weight factor. In a second step, all length values are summed and "
"divided by the area of the circular neighborhood. This process is repeated for all raster cells."
);
}

QgsLineDensityAlgorithm *QgsLineDensityAlgorithm::createInstance() const
{
return new QgsLineDensityAlgorithm();
}

bool QgsLineDensityAlgorithm::prepareAlgorithm(const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback)
{
Q_UNUSED(feedback);
mSource.reset( parameterAsSource( parameters, QStringLiteral( "INPUT" ), context ) );
if ( !mSource )
throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "INPUT" ) ) );
mWeight = parameterAsString( parameters, QStringLiteral( "WEIGHT" ), context );
mSearchRadius = parameterAsDouble(parameters, QStringLiteral("RADIUS"), context);
mPixelSize = parameterAsDouble(parameters, QStringLiteral("PIXEL_SIZE"), context);

mExtent = mSource->sourceExtent();
mCrs = mSource->sourceCrs();

if ( context.project() )
{
mDa.setEllipsoid( context.project()->ellipsoid() );
}
mDa.setSourceCrs( mCrs, context.transformContext() );


//get cell midpoint from top left cell
QgsGeometry firstCellMidpoint = QgsGeometry(new QgsPoint(mExtent.xMinimum() + (mPixelSize/2), mExtent.yMaximum() - (mPixelSize/2)));
mSearchGeometry = firstCellMidpoint.buffer(mSearchRadius, 5);
mSearchGeometryArea = mDa.measureArea(mSearchGeometry);

mIndex = QgsSpatialIndex( *mSource, nullptr, QgsSpatialIndex::FlagStoreFeatureGeometries );


return true;
}

QVariantMap QgsLineDensityAlgorithm::processAlgorithm(const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback)
{
const QString outputFile = parameterAsOutputLayer( parameters, QStringLiteral( "OUTPUT" ), context );
QFileInfo fi( outputFile );
const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() );

int rows = std::max( std::ceil( mExtent.height() / mPixelSize ), 1.0 );
int cols = std::max( std::ceil( mExtent.width() / mPixelSize ), 1.0 );

//build new raster extent based on number of columns and cellsize
//this prevents output cellsize being calculated too small
QgsRectangle rasterExtent = QgsRectangle(mExtent.xMinimum(), mExtent.yMaximum() - (rows * mPixelSize), mExtent.xMinimum() + (cols * mPixelSize), mExtent.yMaximum());

std::unique_ptr< QgsRasterFileWriter > writer = qgis::make_unique< QgsRasterFileWriter >( outputFile );
writer->setOutputProviderKey( QStringLiteral( "gdal" ) );
writer->setOutputFormat( outputFormat );
std::unique_ptr<QgsRasterDataProvider > provider( writer->createOneBandRaster( Qgis::Float32, cols, rows, rasterExtent, mCrs ) );
if ( !provider )
throw QgsProcessingException( QObject::tr( "Could not create raster output: %1" ).arg( outputFile ) );
if ( !provider->isValid() )
throw QgsProcessingException( QObject::tr( "Could not create raster output %1: %2" ).arg( outputFile, provider->error().message( QgsErrorMessage::Text ) ) );

provider->setNoDataValue( 1, -9999 );

int totalCellcnt = rows * cols;
int cellcnt = 0;

for (int row = 0; row < rows; row++)
{
std::unique_ptr< QgsRasterBlock > rasterDataLine = qgis::make_unique< QgsRasterBlock >( Qgis::Float32, cols, 1 );
for (int col = 0; col < cols; col++)
{
if ( feedback->isCanceled() )
{
break;
}

if(col > 0)
mSearchGeometry.translate( mPixelSize, 0 );

QList<QgsFeatureId> fids = mIndex.intersects(mSearchGeometry.boundingBox());

std::unique_ptr< QgsGeometryEngine > engine( QgsGeometry::createGeometryEngine( mSearchGeometry.constGet() ) );
engine->prepareGeometry();

//remove non-intersecting fids
int i = 0;
while( i < fids.size() )
{
QgsFeatureId fid = fids.at(i);
QgsGeometry lineGeom = mIndex.geometry(fid);

if(engine->intersects(lineGeom.constGet()) == false)
{
fids.removeAt(i);
}
i++;
}

QgsFeatureIds isectingFids = QSet< QgsFeatureId >::fromList( fids );
QgsFeatureIterator fit = mSource->getFeatures(QgsFeatureRequest().setFilterFids(isectingFids));
QgsFeature f;
double abs_density = 0;
while (fit.nextFeature( f ))
{
double analysisLineLength = mDa.measureLength(QgsGeometry(engine->intersection(f.geometry().constGet())));
double analysisWeight = f.attribute( mWeight ).toDouble();
abs_density += ( analysisLineLength * analysisWeight );
}

double lineDensity = abs_density/mSearchGeometryArea;
rasterDataLine->setValue(0, col, lineDensity);

feedback->setProgress(static_cast<double>(cellcnt)/static_cast<double>(totalCellcnt) * 100);
cellcnt++;
}
provider->writeBlock( rasterDataLine.get(), 1, 0, row);

//'carriage return and newline' for search geometry
mSearchGeometry.translate((cols - 1) * -mPixelSize, -mPixelSize);
}

QVariantMap outputs;
outputs.insert( QStringLiteral( "OUTPUT" ), outputFile);
return outputs;
}


///@endcond
@@ -0,0 +1,75 @@
/***************************************************************************
qgsalgorithmlinedensity.h
---------------------
begin : December 2019
copyright : (C) 2019 by Clemens Raffler
email : clemens dot raffler 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 QGSALGORITHMLINEDENSITY_H
#define QGSALGORITHMLINEDENSITY_H

#define SIP_NO_FILE

#include "qgis_sip.h"
#include "qgsprocessingalgorithm.h"
#include "qgsapplication.h"

///@cond PRIVATE

/**
* Line Density Algorithm as implemented in ESRI ArcGIS Spatial Analyst
*
* Literature:
* Silverman, B.w. Density Estimation for Statistics and Data Analsis.
* New York: Chapman and Hall, 1986
*
*/
class QgsLineDensityAlgorithm : public QgsProcessingAlgorithm
{
public:

QgsLineDensityAlgorithm() = default;
void initAlgorithm( const QVariantMap &configuration = QVariantMap() ) override;
QIcon icon() const override { return QgsApplication::getThemeIcon( QStringLiteral( "/algorithms/mAlgorithmLineDensity.svg" ) ); }
QString svgIconPath() const override { return QgsApplication::iconPath( QStringLiteral( "/algorithms/mAlgorithmLineDensity.svg" ) ); }
QString name() const override;
QString displayName() const override;
QStringList tags() const override;
QString group() const override;
QString groupId() const override;
QString shortHelpString() const override;
QgsLineDensityAlgorithm *createInstance() const override SIP_FACTORY;

protected:
bool prepareAlgorithm(const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback);
QVariantMap processAlgorithm( const QVariantMap &parameters,
QgsProcessingContext &context,
QgsProcessingFeedback *feedback ) override;

private:
std::unique_ptr< QgsFeatureSource > mSource;
QString mWeight;
double mSearchRadius;
double mSearchGeometryArea;
double mPixelSize;
QgsGeometry mSearchGeometry;
QgsRectangle mExtent;
QgsCoordinateReferenceSystem mCrs;
QgsDistanceArea mDa;
QgsSpatialIndex mIndex;

};

///@endcond PRIVATE

#endif // QGSALGORITHMLINEDENSITY_H
@@ -74,6 +74,7 @@
#include "qgsalgorithminterpolatepoint.h"
#include "qgsalgorithmintersection.h"
#include "qgsalgorithmkmeansclustering.h"
#include "qgsalgorithmlinedensity.h"
#include "qgsalgorithmlineintersection.h"
#include "qgsalgorithmlinesubstring.h"
#include "qgsalgorithmloadlayer.h"
@@ -264,6 +265,7 @@ void QgsNativeAlgorithms::loadAlgorithms()
addAlgorithm( new QgsKMeansClusteringAlgorithm() );
addAlgorithm( new QgsLayerToBookmarksAlgorithm() );
addAlgorithm( new QgsLayoutMapExtentToLayerAlgorithm() );
addAlgorithm( new QgsLineDensityAlgorithm() );
addAlgorithm( new QgsLineIntersectionAlgorithm() );
addAlgorithm( new QgsLineSubstringAlgorithm() );
addAlgorithm( new QgsLoadLayerAlgorithm() );

0 comments on commit ce342f7

Please sign in to comment.
You can’t perform that action at this time.