Skip to content
Permalink
Browse files

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

  • Loading branch information
ghtmtt committed Jul 10, 2018
1 parent 2688a9d commit d1cedbcf92f2679fdcfbb3d3f63bfd5cf638ae45
@@ -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.

@@ -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
@@ -229,6 +230,7 @@ def getAlgs(self):
RasterCalculator(),
RasterizeAlgorithm(),
RasterLayerStatistics(),
RasterSampling(),
RectanglesOvalsDiamondsFixed(),
RectanglesOvalsDiamondsVariable(),
RegularPoints(),
@@ -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}
@@ -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>
@@ -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>
@@ -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

0 comments on commit d1cedbc

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