Skip to content
Permalink
Browse files

[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
  • Loading branch information
nyalldawson committed Jun 14, 2018
1 parent 39f0922 commit 5b7eefa6aed53c0dd42086e88dde28faee0f4ae0

Large diffs are not rendered by default.

@@ -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>
@@ -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
@@ -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

@@ -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


@@ -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


@@ -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() );
}

0 comments on commit 5b7eefa

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