[processing][needs-doc][FEATURE] Sample raster values to point

This commit is contained in:
matteo 2018-07-10 15:19:56 +02:00
parent 2688a9d781
commit d1cedbcf92
6 changed files with 283 additions and 1 deletions

View File

@ -423,6 +423,15 @@ qgis:rasterlayerhistogram: >
qgis:rasterlayerstatistics: >
This algorithm computes basic statistics from the values in a given band of the raster layer.
qgis:rastersampling: >
This algorithm creates a new vector layer with the same attributes of the input layer and the raster values corresponding on the point location.
Many raster layers can be chosen at the same time.
If the raster layer has more than one band, all the band values are sampled.
WARNING: raster layer(s) must be in the same CRS of the source point layer.
qgis:rasterize: >
This algorithm rasterizes map canvas content.
@ -558,4 +567,3 @@ qgis:definecurrentprojection: >
This algorithm sets an existing layer's projection to the provided CRS. Contrary to the "Assign projection" algorithm, it will not output a new layer.
For shapefile datasets, the .prj and .qpj files will be overwritten - or created if missing - to match the provided CRS.

View File

@ -113,6 +113,7 @@ from .RandomSelectionWithinSubsets import RandomSelectionWithinSubsets
from .Rasterize import RasterizeAlgorithm
from .RasterCalculator import RasterCalculator
from .RasterLayerStatistics import RasterLayerStatistics
from .RasterSampling import RasterSampling
from .RectanglesOvalsDiamondsFixed import RectanglesOvalsDiamondsFixed
from .RectanglesOvalsDiamondsVariable import RectanglesOvalsDiamondsVariable
from .RegularPoints import RegularPoints
@ -229,6 +230,7 @@ class QgisAlgorithmProvider(QgsProcessingProvider):
RasterCalculator(),
RasterizeAlgorithm(),
RasterLayerStatistics(),
RasterSampling(),
RectanglesOvalsDiamondsFixed(),
RectanglesOvalsDiamondsVariable(),
RegularPoints(),

View File

@ -0,0 +1,163 @@
# -*- coding: utf-8 -*-
"""
***************************************************************************
RasterSampling.py
-----------------------
Date : July 2018
Copyright : (C) 2018 by Matteo Ghetta
Email : matteo dot ghetta 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. *
* *
***************************************************************************
"""
__author__ = 'Matteo Ghetta'
__date__ = 'July 2018'
__copyright__ = '(C) 2018, Matteo Ghetta'
# This will get replaced with a git SHA1 when you do a git archive
__revision__ = '$Format:%H$'
import os
from qgis.PyQt.QtGui import QIcon
from qgis.PyQt.QtCore import QVariant
from qgis.core import (QgsApplication,
QgsField,
QgsFeatureSink,
QgsRaster,
QgsProcessing,
QgsProcessingParameterMultipleLayers,
QgsFields,
QgsProcessingUtils,
QgsProcessingException,
QgsProcessingParameterFeatureSource,
QgsProcessingParameterFeatureSink)
from processing.algs.qgis.QgisAlgorithm import QgisAlgorithm
class RasterSampling(QgisAlgorithm):
INPUT = 'INPUT'
RASTERCOPY = 'RASTERCOPY'
OUTPUT = 'OUTPUT'
def name(self):
return 'rastersampling'
def displayName(self):
return self.tr('Sample Raster Values')
def group(self):
return self.tr('Raster analysis')
def groupId(self):
return 'rasteranalysis'
def __init__(self):
super().__init__()
def initAlgorithm(self, config=None):
self.addParameter(
QgsProcessingParameterFeatureSource(
self.INPUT,
self.tr('Input Point Layer'),
[QgsProcessing.TypeVectorPoint]
)
)
self.addParameter(
QgsProcessingParameterMultipleLayers(
self.RASTERCOPY,
self.tr('Raster Layer to sample'),
QgsProcessing.TypeRaster
)
)
self.addParameter(
QgsProcessingParameterFeatureSink(
self.OUTPUT,
self.tr('Sampled Points')
)
)
def processAlgorithm(self, parameters, context, feedback):
source = self.parameterAsSource(
parameters,
self.INPUT,
context
)
sampled_rasters = self.parameterAsLayerList(
parameters,
self.RASTERCOPY,
context
)
if source is None:
raise QgsProcessingException(self.invalidSourceError(parameters, self.INPUT))
source_fields = source.fields()
raster_fields = QgsFields()
# append field to vector as rasterName_bandCount
for i in sampled_rasters:
for b in range(i.bandCount()):
raster_fields.append(QgsField(i.name() + str('_{}'.format(b+1)), QVariant.Double))
# combine all the vector fields
out_fields = QgsProcessingUtils.combineFields(source_fields, raster_fields)
(sink, dest_id) = self.parameterAsSink(
parameters,
self.OUTPUT,
context,
out_fields,
source.wkbType(),
source.sourceCrs()
)
if sink is None:
raise QgsProcessingException(self.invalidSinkError(parameters, self.OUTPUT))
total = 100.0 / source.featureCount() if source.featureCount() else 0
features = source.getFeatures()
for n, i in enumerate(source.getFeatures()):
for rr in sampled_rasters:
attrs = i.attributes()
if rr.bandCount() >1:
for b in range(rr.bandCount()):
attrs.append(rr.dataProvider().identify(i.geometry().asPoint(),
QgsRaster.IdentifyFormatValue).results()[b+1])
attrs.append(rr.dataProvider().identify(i.geometry().asPoint(), QgsRaster.IdentifyFormatValue).results()[1])
i.setAttributes(attrs)
sink.addFeature(i, QgsFeatureSink.FastInsert)
feedback.setProgress(int(n * total))
return {self.OUTPUT: dest_id}

View File

@ -0,0 +1,56 @@
<?xml version="1.0" encoding="utf-8" ?>
<ogr:FeatureCollection
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:schemaLocation="http://ogr.maptools.org/ raster_sampling.xsd"
xmlns:ogr="http://ogr.maptools.org/"
xmlns:gml="http://www.opengis.net/gml">
<gml:boundedBy>
<gml:Box>
<gml:coord><gml:X>18.67356158120081</gml:X><gml:Y>45.78707029401009</gml:Y></gml:coord>
<gml:coord><gml:X>18.6923582796333</gml:X><gml:Y>45.80603482010714</gml:Y></gml:coord>
</gml:Box>
</gml:boundedBy>
<gml:featureMember>
<ogr:raster_sampling fid="raster_sampling.0">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>18.6735615812008,45.8060348201071</gml:coordinates></gml:Point></ogr:geometryProperty>
<ogr:id>1</ogr:id>
<ogr:dem_1>91.9273986816406</ogr:dem_1>
</ogr:raster_sampling>
</gml:featureMember>
<gml:featureMember>
<ogr:raster_sampling fid="raster_sampling.1">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>18.6891695540064,45.8035174051385</gml:coordinates></gml:Point></ogr:geometryProperty>
<ogr:id>2</ogr:id>
<ogr:dem_1>122.21053314209</ogr:dem_1>
</ogr:raster_sampling>
</gml:featureMember>
<gml:featureMember>
<ogr:raster_sampling fid="raster_sampling.2">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>18.6923582796333,45.7875737770038</gml:coordinates></gml:Point></ogr:geometryProperty>
<ogr:id>3</ogr:id>
<ogr:dem_1>129.417877197266</ogr:dem_1>
</ogr:raster_sampling>
</gml:featureMember>
<gml:featureMember>
<ogr:raster_sampling fid="raster_sampling.3">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>18.6796033771255,45.7870702940101</gml:coordinates></gml:Point></ogr:geometryProperty>
<ogr:id>4</ogr:id>
<ogr:dem_1>109.654373168945</ogr:dem_1>
</ogr:raster_sampling>
</gml:featureMember>
<gml:featureMember>
<ogr:raster_sampling fid="raster_sampling.4">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>18.6812816537713,45.7964686432263</gml:coordinates></gml:Point></ogr:geometryProperty>
<ogr:id>5</ogr:id>
<ogr:dem_1>174.289291381836</ogr:dem_1>
</ogr:raster_sampling>
</gml:featureMember>
<gml:featureMember>
<ogr:raster_sampling fid="raster_sampling.5">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>18.6735615812008,45.7926086069411</gml:coordinates></gml:Point></ogr:geometryProperty>
<ogr:id>6</ogr:id>
<ogr:dem_1>93.1188354492188</ogr:dem_1>
</ogr:raster_sampling>
</gml:featureMember>
</ogr:FeatureCollection>

View File

@ -0,0 +1,36 @@
<?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="raster_sampling" type="ogr:raster_sampling_Type" substitutionGroup="gml:_Feature"/>
<xs:complexType name="raster_sampling_Type">
<xs:complexContent>
<xs:extension base="gml:AbstractFeatureType">
<xs:sequence>
<xs:element name="geometryProperty" type="gml:PointPropertyType" nillable="true" minOccurs="0" maxOccurs="1"/>
<xs:element name="id" nillable="true" minOccurs="0" maxOccurs="1">
<xs:simpleType>
<xs:restriction base="xs:long">
<xs:totalDigits value="10"/>
</xs:restriction>
</xs:simpleType>
</xs:element>
<xs:element name="dem_1" nillable="true" minOccurs="0" maxOccurs="1">
<xs:simpleType>
<xs:restriction base="xs:decimal">
</xs:restriction>
</xs:simpleType>
</xs:element>
</xs:sequence>
</xs:extension>
</xs:complexContent>
</xs:complexType>
</xs:schema>

View File

@ -5728,4 +5728,21 @@ tests:
name: expected/kmeans_polys.gml
type: vector
- algorithm: qgis:rastersampling
name: Raster sampling, points single band
params:
INPUT:
name: custom/points_over.shp
type: vector
RASTERCOPY:
params:
- name: dem.tif
type: raster
type: multi
results:
OUTPUT:
name: expected/raster_sampling.gml
type: vector
# See ../README.md for a description of the file format