[FEATURE][processing] New 'Raster pixels to polygons' algorithm

Converts a raster layer into a vector layer, with a polygon feature corresponding to each pixel from the raster and a single field containing the band value from the raster. Sponsored by SMEC/SJ
qgis · Jun 22, 2018 · 5b7eefa · 5b7eefa
1 parent 39f0922
commit 5b7eefa
Show file tree

Hide file tree

Showing 7 changed files with 1,019 additions and 1 deletion.
diff --git a/python/plugins/processing/tests/testdata/expected/vectorize.gml b/python/plugins/processing/tests/testdata/expected/vectorize.gml
diff --git a/python/plugins/processing/tests/testdata/expected/vectorize.xsd b/python/plugins/processing/tests/testdata/expected/vectorize.xsd
@@ -0,0 +1,31 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<xs:schema targetNamespace="http://ogr.maptools.org/" xmlns:ogr="http://ogr.maptools.org/" xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:gml="http://www.opengis.net/gml" elementFormDefault="qualified" version="1.0">
+<xs:import namespace="http://www.opengis.net/gml" schemaLocation="http://schemas.opengis.net/gml/2.1.2/feature.xsd"/>
+<xs:element name="FeatureCollection" type="ogr:FeatureCollectionType" substitutionGroup="gml:_FeatureCollection"/>
+<xs:complexType name="FeatureCollectionType">
+  <xs:complexContent>
+    <xs:extension base="gml:AbstractFeatureCollectionType">
+      <xs:attribute name="lockId" type="xs:string" use="optional"/>
+      <xs:attribute name="scope" type="xs:string" use="optional"/>
+    </xs:extension>
+  </xs:complexContent>
+</xs:complexType>
+<xs:element name="vectorize" type="ogr:vectorize_Type" substitutionGroup="gml:_Feature"/>
+<xs:complexType name="vectorize_Type">
+  <xs:complexContent>
+    <xs:extension base="gml:AbstractFeatureType">
+      <xs:sequence>
+        <xs:element name="geometryProperty" type="gml:PolygonPropertyType" nillable="true" minOccurs="0" maxOccurs="1"/>
+        <xs:element name="pix_val" nillable="true" minOccurs="0" maxOccurs="1">
+          <xs:simpleType>
+            <xs:restriction base="xs:decimal">
+              <xs:totalDigits value="21"/>
+              <xs:fractionDigits value="8"/>
+            </xs:restriction>
+          </xs:simpleType>
+        </xs:element>
+      </xs:sequence>
+    </xs:extension>
+  </xs:complexContent>
+</xs:complexType>
+</xs:schema>
diff --git a/python/plugins/processing/tests/testdata/qgis_algorithm_tests.yaml b/python/plugins/processing/tests/testdata/qgis_algorithm_tests.yaml
@@ -5647,6 +5647,17 @@ tests:
         hash: c29d14f71e8686f7445d53be646fce84702644f159fd0164ac38e861
         type: rasterhash
 
-
+  - algorithm: native:pixelstopolygons
+    name: Pixels to polygons
+    params:
+      FIELD_NAME: pix_val
+      INPUT_RASTER:
+        name: raster.tif
+        type: raster
+      RASTER_BAND: 1
+    results:
+      OUTPUT:
+        name: expected/vectorize.gml
+        type: vector
 
 # See ../README.md for a description of the file format
diff --git a/src/analysis/CMakeLists.txt b/src/analysis/CMakeLists.txt
@@ -82,6 +82,7 @@ SET(QGIS_ANALYSIS_SRCS
   processing/qgsalgorithmtranslate.cpp
   processing/qgsalgorithmunion.cpp
   processing/qgsalgorithmuniquevalueindex.cpp
+  processing/qgsalgorithmvectorize.cpp
   processing/qgsalgorithmwedgebuffers.cpp
   processing/qgsalgorithmzonalhistogram.cpp
 

diff --git a/src/analysis/processing/qgsalgorithmvectorize.cpp b/src/analysis/processing/qgsalgorithmvectorize.cpp
@@ -0,0 +1,167 @@
+/***************************************************************************
+                         qgsalgorithmvectorize.cpp
+                         ---------------------
+    begin                : June, 2018
+    copyright            : (C) 2018 by Nyall Dawson
+    email                : nyall dot dawson 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 "qgsalgorithmvectorize.h"
+#include "qgis.h"
+#include "qgsprocessing.h"
+
+///@cond PRIVATE
+
+QString QgsVectorizeAlgorithm::name() const
+{
+  return QStringLiteral( "pixelstopolygons" );
+}
+
+QString QgsVectorizeAlgorithm::displayName() const
+{
+  return QObject::tr( "Raster pixels to polygons" );
+}
+
+QStringList QgsVectorizeAlgorithm::tags() const
+{
+  return QObject::tr( "vectorize,polygonize,raster,convert,pixels" ).split( ',' );
+}
+
+QString QgsVectorizeAlgorithm::shortHelpString() const
+{
+  return QObject::tr( "This algorithm converts a raster layer to a vector layer, by creating polygon features for each individual pixel in the raster layer." );
+}
+
+QgsVectorizeAlgorithm *QgsVectorizeAlgorithm::createInstance() const
+{
+  return new QgsVectorizeAlgorithm();
+}
+
+QString QgsVectorizeAlgorithm::group() const
+{
+  return QObject::tr( "Vector creation" );
+}
+
+QString QgsVectorizeAlgorithm::groupId() const
+{
+  return QStringLiteral( "vectorcreation" );
+}
+
+void QgsVectorizeAlgorithm::initAlgorithm( const QVariantMap & )
+{
+  addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "INPUT_RASTER" ),
+                QObject::tr( "Raster layer" ) ) );
+  addParameter( new QgsProcessingParameterBand( QStringLiteral( "RASTER_BAND" ),
+                QObject::tr( "Band number" ), 1, QStringLiteral( "INPUT_RASTER" ) ) );
+  addParameter( new QgsProcessingParameterString( QStringLiteral( "FIELD_NAME" ),
+                QObject::tr( "Field name" ), QStringLiteral( "VALUE" ) ) );
+
+  addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT" ), QObject::tr( "Vectorized layer" ), QgsProcessing::TypeVectorPolygon ) );
+}
+
+bool QgsVectorizeAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback * )
+{
+  QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral( "INPUT_RASTER" ), context );
+
+  if ( !layer )
+    throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "INPUT_RASTER" ) ) );
+
+  mBand = parameterAsInt( parameters, QStringLiteral( "RASTER_BAND" ), context );
+  if ( mBand < 1 || mBand > layer->bandCount() )
+    throw QgsProcessingException( QObject::tr( "Invalid band number for RASTER_BAND (%1): Valid values for input raster are 1 to %2" ).arg( mBand )
+                                  .arg( layer->bandCount() ) );
+
+  mInterface.reset( layer->dataProvider()->clone() );
+  mExtent = layer->extent();
+  mCrs = layer->crs();
+  mRasterUnitsPerPixelX = std::abs( layer->rasterUnitsPerPixelX() );
+  mRasterUnitsPerPixelY = std::abs( layer->rasterUnitsPerPixelY() );
+  mNbCellsXProvider = mInterface->xSize();
+  mNbCellsYProvider = mInterface->ySize();
+  return true;
+}
+
+QVariantMap QgsVectorizeAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
+{
+  const QString fieldName = parameterAsString( parameters, QStringLiteral( "FIELD_NAME" ), context );
+  QgsFields fields;
+  fields.append( QgsField( fieldName, QVariant::Double, QString(), 20, 8 ) );
+
+  QString dest;
+  std::unique_ptr< QgsFeatureSink > sink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, dest, fields, QgsWkbTypes::Polygon, mCrs ) );
+  if ( !sink )
+    throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT" ) ) );
+
+
+  int maxWidth = QgsRasterIterator::DEFAULT_MAXIMUM_TILE_WIDTH;
+  int maxHeight = QgsRasterIterator::DEFAULT_MAXIMUM_TILE_HEIGHT;
+
+  QgsRasterIterator iter( mInterface.get() );
+  iter.startRasterRead( mBand, mNbCellsXProvider, mNbCellsYProvider, mExtent );
+
+  int nbBlocksWidth = static_cast< int >( std::ceil( 1.0 * mNbCellsXProvider / maxWidth ) );
+  int nbBlocksHeight = static_cast< int >( std::ceil( 1.0 * mNbCellsYProvider / maxHeight ) );
+  int nbBlocks = nbBlocksWidth * nbBlocksHeight;
+
+  double hCellSizeX = mRasterUnitsPerPixelX / 2.0;
+  double hCellSizeY = mRasterUnitsPerPixelY / 2.0;
+
+  int iterLeft = 0;
+  int iterTop = 0;
+  int iterCols = 0;
+  int iterRows = 0;
+  std::unique_ptr< QgsRasterBlock > rasterBlock;
+  QgsRectangle blockExtent;
+  while ( iter.readNextRasterPart( mBand, iterCols, iterRows, rasterBlock, iterLeft, iterTop, &blockExtent ) )
+  {
+    if ( feedback )
+      feedback->setProgress( 100 * ( ( iterTop / maxHeight * nbBlocksWidth ) + iterLeft / maxWidth ) / nbBlocks );
+    if ( feedback && feedback->isCanceled() )
+      break;
+
+    double currentY = blockExtent.yMaximum() - 0.5 * mRasterUnitsPerPixelY;
+
+    for ( int row = 0; row < iterRows; row++ )
+    {
+      if ( feedback && feedback->isCanceled() )
+        break;
+
+      double currentX = blockExtent.xMinimum() + 0.5 * mRasterUnitsPerPixelX;
+
+      for ( int column = 0; column < iterCols; column++ )
+      {
+        if ( !rasterBlock->isNoData( row, column ) )
+        {
+
+          QgsGeometry pixelRectGeometry = QgsGeometry::fromRect( QgsRectangle( currentX - hCellSizeX, currentY - hCellSizeY, currentX + hCellSizeX, currentY + hCellSizeY ) );
+
+          double value = rasterBlock->value( row, column );
+
+          QgsFeature f;
+          f.setGeometry( pixelRectGeometry );
+          f.setAttributes( QgsAttributes() << value );
+          sink->addFeature( f, QgsFeatureSink::FastInsert );
+        }
+        currentX += mRasterUnitsPerPixelX;
+      }
+      currentY -= mRasterUnitsPerPixelY;
+    }
+  }
+
+  QVariantMap outputs;
+  outputs.insert( QStringLiteral( "OUTPUT" ), dest );
+  return outputs;
+}
+
+///@endcond
+
+
diff --git a/src/analysis/processing/qgsalgorithmvectorize.h b/src/analysis/processing/qgsalgorithmvectorize.h
@@ -0,0 +1,72 @@
+/***************************************************************************
+                         qgsalgorithmvectorize.h
+                         ---------------------
+    begin                : June, 2018
+    copyright            : (C) 2018 by Nyall Dawson
+    email                : nyall dot dawson 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 QGSALGORITHMVECTORIZE_H
+#define QGSALGORITHMVECTORIZE_H
+
+#define SIP_NO_FILE
+
+#include "qgis.h"
+#include "qgsprocessingalgorithm.h"
+#include "qgsreclassifyutils.h"
+
+///@cond PRIVATE
+
+/**
+ * Native vectorize algorithm
+ */
+class QgsVectorizeAlgorithm : public QgsProcessingAlgorithm
+{
+  public:
+
+    QgsVectorizeAlgorithm() = default;
+    QString name() const override;
+    QString displayName() const override;
+    QStringList tags() const override;
+    QString shortHelpString() const override;
+    QgsVectorizeAlgorithm *createInstance() const override SIP_FACTORY;
+    QString group() const override final;
+    QString groupId() const override final;
+    void initAlgorithm( const QVariantMap &configuration = QVariantMap() ) override;
+
+  protected:
+
+    bool prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override;
+
+    QVariantMap processAlgorithm( const QVariantMap &parameters,
+                                  QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override;
+
+    std::unique_ptr< QgsRasterInterface > mInterface;
+
+    Qgis::DataType mDataType = Qgis::Float32;
+    double mNoDataValue = -9999;
+    int mBand = 1;
+    QgsRectangle mExtent;
+    QgsCoordinateReferenceSystem mCrs;
+    double mRasterUnitsPerPixelX = 0;
+    double mRasterUnitsPerPixelY = 0;
+    int mNbCellsXProvider = 0;
+    int mNbCellsYProvider = 0;
+    QgsReclassifyUtils::RasterClass::BoundsType mBoundsType = QgsReclassifyUtils::RasterClass::IncludeMax;
+    bool mUseNoDataForMissingValues = false;
+};
+
+///@endcond PRIVATE
+
+#endif // QGSALGORITHMVECTORIZE_H
+
+
diff --git a/src/analysis/processing/qgsnativealgorithms.cpp b/src/analysis/processing/qgsnativealgorithms.cpp
@@ -79,6 +79,7 @@
 #include "qgsalgorithmtranslate.h"
 #include "qgsalgorithmunion.h"
 #include "qgsalgorithmuniquevalueindex.h"
+#include "qgsalgorithmvectorize.h"
 #include "qgsalgorithmwedgebuffers.h"
 #include "qgsalgorithmzonalhistogram.h"
 
@@ -189,6 +190,7 @@ void QgsNativeAlgorithms::loadAlgorithms()
   addAlgorithm( new QgsTranslateAlgorithm() );
   addAlgorithm( new QgsUnionAlgorithm() );
   addAlgorithm( new QgsVariableWidthBufferByMAlgorithm() );
+  addAlgorithm( new QgsVectorizeAlgorithm() );
   addAlgorithm( new QgsWedgeBuffersAlgorithm() );
   addAlgorithm( new QgsZonalHistogramAlgorithm() );
 }