Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
Merge pull request #7393 from ghtmtt/processing_sampling
[processing][FEATURE] Sample raster values to point
  • Loading branch information
nyalldawson committed Jul 16, 2018
2 parents f9fc708 + 14c74d8 commit 0c62abe
Show file tree
Hide file tree
Showing 8 changed files with 418 additions and 1 deletion.
6 changes: 5 additions & 1 deletion python/plugins/processing/algs/help/qgis.yaml
Expand Up @@ -423,6 +423,11 @@ 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.

If the raster layer has more than one band, all the band values are sampled.

qgis:rasterize: >
This algorithm rasterizes map canvas content.

Expand Down Expand Up @@ -558,4 +563,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.

2 changes: 2 additions & 0 deletions python/plugins/processing/algs/qgis/QgisAlgorithmProvider.py
Expand Up @@ -113,6 +113,7 @@
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
Expand Down Expand Up @@ -229,6 +230,7 @@ def getAlgs(self):
RasterCalculator(),
RasterizeAlgorithm(),
RasterLayerStatistics(),
RasterSampling(),
RectanglesOvalsDiamondsFixed(),
RectanglesOvalsDiamondsVariable(),
RegularPoints(),
Expand Down
202 changes: 202 additions & 0 deletions python/plugins/processing/algs/qgis/RasterSampling.py
@@ -0,0 +1,202 @@
# -*- 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,
QgsProcessingParameterRasterLayer,
QgsProcessingParameterString,
QgsProcessingParameterDefinition,
QgsCoordinateTransform,
QgsFields,
QgsProcessingUtils,
QgsProcessingException,
QgsProcessingParameterFeatureSource,
QgsProcessingParameterFeatureSink)

from processing.algs.qgis.QgisAlgorithm import QgisAlgorithm


class RasterSampling(QgisAlgorithm):

INPUT = 'INPUT'
RASTERCOPY = 'RASTERCOPY'
COLUMN_PREFIX = 'COLUMN_PREFIX'
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(
QgsProcessingParameterRasterLayer(
self.RASTERCOPY,
self.tr('Raster Layer to sample'),
)
)

columnPrefix = QgsProcessingParameterString(
self.COLUMN_PREFIX,
self.tr('Output column prefix'), 'rvalue'
)
columnPrefix.setFlags(columnPrefix.flags() | QgsProcessingParameterDefinition.FlagAdvanced)
self.addParameter(columnPrefix)

self.addParameter(
QgsProcessingParameterFeatureSink(
self.OUTPUT,
self.tr('Sampled Points')
)
)

def processAlgorithm(self, parameters, context, feedback):

source = self.parameterAsSource(
parameters,
self.INPUT,
context
)

sampled_raster = self.parameterAsRasterLayer(
parameters,
self.RASTERCOPY,
context
)

columnPrefix = self.parameterAsString(
parameters,
self.COLUMN_PREFIX,
context
)

if source is None:
raise QgsProcessingException(self.invalidSourceError(parameters, self.INPUT))

source_fields = source.fields()
raster_fields = QgsFields()

# append field to vector as columnPrefix_bandCount
for b in range(sampled_raster.bandCount()):
raster_fields.append(QgsField(
columnPrefix + 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()

# create the coordinates transformation context
ct = QgsCoordinateTransform(source.sourceCrs(), sampled_raster.crs(), context.transformContext())

for n, i in enumerate(source.getFeatures()):

attrs = i.attributes()

if i.geometry().isMultipart():
raise QgsProcessingException(self.tr('''Impossible to sample data
of a Multipart layer. Please use the Multipart to single part
algorithm to transform the layer.'''))

# get the feature geometry as point
point = i.geometry().asPoint()

# reproject to raster crs
try:
point = ct.transform(point)
except QgsCsException:
for b in range(sampled_raster.bandCount()):
attrs.append(None)
i.setAttributes(attrs)
sink.addFeature(i, QgsFeatureSink.FastInsert)
feedback.setProgress(int(n * total))
feedback.reportError(self.tr('Could not reproject feature {} to raster CRS').format(i.id()))
continue

if sampled_raster.bandCount() > 1:

for b in range(sampled_raster.bandCount()):
attrs.append(
sampled_raster.dataProvider().identify(
point,
QgsRaster.IdentifyFormatValue).results()[b + 1]
)

attrs.append(
sampled_raster.dataProvider().identify(
point,
QgsRaster.IdentifyFormatValue).results()[1]
)

i.setAttributes(attrs)

sink.addFeature(i, QgsFeatureSink.FastInsert)
feedback.setProgress(int(n * total))

return {self.OUTPUT: dest_id}
@@ -0,0 +1,64 @@
<?xml version="1.0" encoding="utf-8" ?>
<ogr:FeatureCollection
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:schemaLocation="http://ogr.maptools.org/ sampling_points.xsd"
xmlns:ogr="http://ogr.maptools.org/"
xmlns:gml="http://www.opengis.net/gml">
<gml:boundedBy>
<gml:Box>
<gml:coord><gml:X>18.66907337319408</gml:X><gml:Y>45.78164855298462</gml:Y></gml:coord>
<gml:coord><gml:X>18.69811881698106</gml:X><gml:Y>45.80559559440474</gml:Y></gml:coord>
</gml:Box>
</gml:boundedBy>

<gml:featureMember>
<ogr:sampling_points fid="sampling_points.0">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>18.6716225743775,45.8055955944047</gml:coordinates></gml:Point></ogr:geometryProperty>
</ogr:sampling_points>
</gml:featureMember>
<gml:featureMember>
<ogr:sampling_points fid="sampling_points.1">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>18.6825918643183,45.8040506239905</gml:coordinates></gml:Point></ogr:geometryProperty>
</ogr:sampling_points>
</gml:featureMember>
<gml:featureMember>
<ogr:sampling_points fid="sampling_points.2">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>18.6905484619515,45.7973300026888</gml:coordinates></gml:Point></ogr:geometryProperty>
</ogr:sampling_points>
</gml:featureMember>
<gml:featureMember>
<ogr:sampling_points fid="sampling_points.3">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>18.6882310063302,45.7870559494343</gml:coordinates></gml:Point></ogr:geometryProperty>
</ogr:sampling_points>
</gml:featureMember>
<gml:featureMember>
<ogr:sampling_points fid="sampling_points.4">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>18.6778024560343,45.7857427245822</gml:coordinates></gml:Point></ogr:geometryProperty>
</ogr:sampling_points>
</gml:featureMember>
<gml:featureMember>
<ogr:sampling_points fid="sampling_points.5">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>18.6833643495254,45.7918453577183</gml:coordinates></gml:Point></ogr:geometryProperty>
</ogr:sampling_points>
</gml:featureMember>
<gml:featureMember>
<ogr:sampling_points fid="sampling_points.6">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>18.6974235802947,45.8035871328663</gml:coordinates></gml:Point></ogr:geometryProperty>
</ogr:sampling_points>
</gml:featureMember>
<gml:featureMember>
<ogr:sampling_points fid="sampling_points.7">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>18.6981188169811,45.7852019849373</gml:coordinates></gml:Point></ogr:geometryProperty>
</ogr:sampling_points>
</gml:featureMember>
<gml:featureMember>
<ogr:sampling_points fid="sampling_points.8">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>18.6690733731941,45.7816485529846</gml:coordinates></gml:Point></ogr:geometryProperty>
</ogr:sampling_points>
</gml:featureMember>
<gml:featureMember>
<ogr:sampling_points fid="sampling_points.9">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>18.6707728406497,45.7915363636355</gml:coordinates></gml:Point></ogr:geometryProperty>
</ogr:sampling_points>
</gml:featureMember>
</ogr:FeatureCollection>
@@ -0,0 +1,23 @@
<?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="sampling_points" type="ogr:sampling_points_Type" substitutionGroup="gml:_Feature"/>
<xs:complexType name="sampling_points_Type">
<xs:complexContent>
<xs:extension base="gml:AbstractFeatureType">
<xs:sequence>
<xs:element name="geometryProperty" type="gml:PointPropertyType" nillable="true" minOccurs="0" maxOccurs="1"/>
</xs:sequence>
</xs:extension>
</xs:complexContent>
</xs:complexType>
</xs:schema>
@@ -0,0 +1,74 @@
<?xml version="1.0" encoding="utf-8" ?>
<ogr:FeatureCollection
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:schemaLocation="http://ogr.maptools.org/ sampling_raster.xsd"
xmlns:ogr="http://ogr.maptools.org/"
xmlns:gml="http://www.opengis.net/gml">
<gml:boundedBy>
<gml:Box>
<gml:coord><gml:X>18.6690733731941</gml:X><gml:Y>45.7816485529846</gml:Y></gml:coord>
<gml:coord><gml:X>18.6981188169811</gml:X><gml:Y>45.8055955944047</gml:Y></gml:coord>
</gml:Box>
</gml:boundedBy>

<gml:featureMember>
<ogr:sampling_raster fid="sampling_points.0">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>18.6716225743775,45.8055955944047</gml:coordinates></gml:Point></ogr:geometryProperty>
<ogr:rvalue_1>91.9159927368164</ogr:rvalue_1>
</ogr:sampling_raster>
</gml:featureMember>
<gml:featureMember>
<ogr:sampling_raster fid="sampling_points.1">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>18.6825918643183,45.8040506239905</gml:coordinates></gml:Point></ogr:geometryProperty>
<ogr:rvalue_1>88.300178527832</ogr:rvalue_1>
</ogr:sampling_raster>
</gml:featureMember>
<gml:featureMember>
<ogr:sampling_raster fid="sampling_points.2">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>18.6905484619515,45.7973300026888</gml:coordinates></gml:Point></ogr:geometryProperty>
<ogr:rvalue_1>182.22380065918</ogr:rvalue_1>
</ogr:sampling_raster>
</gml:featureMember>
<gml:featureMember>
<ogr:sampling_raster fid="sampling_points.3">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>18.6882310063302,45.7870559494343</gml:coordinates></gml:Point></ogr:geometryProperty>
<ogr:rvalue_1>126.56761932373</ogr:rvalue_1>
</ogr:sampling_raster>
</gml:featureMember>
<gml:featureMember>
<ogr:sampling_raster fid="sampling_points.4">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>18.6778024560343,45.7857427245822</gml:coordinates></gml:Point></ogr:geometryProperty>
<ogr:rvalue_1>131.698287963867</ogr:rvalue_1>
</ogr:sampling_raster>
</gml:featureMember>
<gml:featureMember>
<ogr:sampling_raster fid="sampling_points.5">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>18.6833643495254,45.7918453577183</gml:coordinates></gml:Point></ogr:geometryProperty>
<ogr:rvalue_1>150</ogr:rvalue_1>
</ogr:sampling_raster>
</gml:featureMember>
<gml:featureMember>
<ogr:sampling_raster fid="sampling_points.6">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>18.6974235802947,45.8035871328663</gml:coordinates></gml:Point></ogr:geometryProperty>
<ogr:rvalue_1>123.953346252441</ogr:rvalue_1>
</ogr:sampling_raster>
</gml:featureMember>
<gml:featureMember>
<ogr:sampling_raster fid="sampling_points.7">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>18.6981188169811,45.7852019849373</gml:coordinates></gml:Point></ogr:geometryProperty>
<ogr:rvalue_1>188.614822387695</ogr:rvalue_1>
</ogr:sampling_raster>
</gml:featureMember>
<gml:featureMember>
<ogr:sampling_raster fid="sampling_points.8">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>18.6690733731941,45.7816485529846</gml:coordinates></gml:Point></ogr:geometryProperty>
<ogr:rvalue_1>168.492599487305</ogr:rvalue_1>
</ogr:sampling_raster>
</gml:featureMember>
<gml:featureMember>
<ogr:sampling_raster fid="sampling_points.9">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>18.6707728406497,45.7915363636355</gml:coordinates></gml:Point></ogr:geometryProperty>
<ogr:rvalue_1>96.8800735473633</ogr:rvalue_1>
</ogr:sampling_raster>
</gml:featureMember>
</ogr:FeatureCollection>

0 comments on commit 0c62abe

Please sign in to comment.