[feature][processing] Add Line denisty algorithm

qgis · Jan 10, 2020 · ce342f7 · ce342f7
1 parent 5765cab
commit ce342f7
Show file tree

Hide file tree

Showing 8 changed files with 337 additions and 0 deletions.
diff --git a/images/images.qrc b/images/images.qrc
@@ -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>

diff --git a/images/themes/default/algorithms/mAlgorithmLineDensity.svg b/images/themes/default/algorithms/mAlgorithmLineDensity.svg
diff --git a/python/plugins/processing/tests/testdata/linedensity.gml b/python/plugins/processing/tests/testdata/linedensity.gml
@@ -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>
diff --git a/python/plugins/processing/tests/testdata/qgis_algorithm_tests2.yaml b/python/plugins/processing/tests/testdata/qgis_algorithm_tests2.yaml
@@ -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:

diff --git a/src/analysis/CMakeLists.txt b/src/analysis/CMakeLists.txt
@@ -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

diff --git a/src/analysis/processing/qgsalgorithmlinedensity.cpp b/src/analysis/processing/qgsalgorithmlinedensity.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
diff --git a/src/analysis/processing/qgsalgorithmlinedensity.h b/src/analysis/processing/qgsalgorithmlinedensity.h
@@ -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
diff --git a/src/analysis/processing/qgsnativealgorithms.cpp b/src/analysis/processing/qgsnativealgorithms.cpp
@@ -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() );