Skip to content

Commit 2fea23f

Browse files
authoredDec 5, 2016
Merge pull request #3779 from volaya/rastercalculator
[processing] add native raster calculator
2 parents abbe281 + 5c6c18c commit 2fea23f

16 files changed

+801
-103
lines changed
 

‎python/plugins/processing/algs/gdal/GdalAlgorithmDialog.py

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -78,7 +78,8 @@ def __init__(self, parent, alg):
7878
self.parametersHaveChanged()
7979

8080
def connectParameterSignals(self):
81-
for w in list(self.widgets.values()):
81+
for wrapper in list(self.wrappers.values()):
82+
w = wrapper.widget
8283
if isinstance(w, QLineEdit):
8384
w.textChanged.connect(self.parametersHaveChanged)
8485
elif isinstance(w, QComboBox):

‎python/plugins/processing/algs/help/qgis.yaml

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -554,3 +554,18 @@ qgis:voronoipolygons: >
554554

555555
qgis:zonalstatistics:
556556

557+
qgis:rastercalculator: >
558+
This algorithm allows to perform algebraic operations using raster layers.
559+
560+
The resulting layer will have its values computed according to an expression. The expression can contain numerical values, operators and references to any of the layers in the current project. The following functions are also supported:
561+
562+
- sin(), cos(), tan(), atan2(), ln(), log10()
563+
564+
The extent and cellsize can be defined by the user. If the extent is not specified, the minimum extent that covers the input layers will be used. If the cellsize is not specified, the minimum cellsize of all input layers will be used.
565+
566+
The cellsize is assumed to be the same in both X and Y axes.
567+
568+
Layers are refered by their name as displayed in the layer list and the number of the band to use (based on 1), using the pattern 'layer_name@band number'. For instance, the first band from a layer named DEM will be referred as DEM@1.
569+
570+
When using the calculator in the batch interface or from the console, the files to use have to be specified. The corresponding layers are refered using the base name of the file (without the full path). For instance, is using a layer at path/to/my/rasterfile.tif, the first band of that layer will be refered as rasterfile.tif@1.
571+

‎python/plugins/processing/algs/qgis/FieldsMapper.py

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -86,7 +86,6 @@ def setValue(self, value):
8686
return False
8787
return False
8888

89-
9089
self.addParameter(ParameterFieldsMapping(self.FIELDS_MAPPING,
9190
self.tr('Fields mapping'),
9291
self.INPUT_LAYER))

‎python/plugins/processing/algs/qgis/QGISAlgorithmProvider.py

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -178,6 +178,7 @@
178178
from .GeometryByExpression import GeometryByExpression
179179
from .SnapGeometries import SnapGeometriesToLayer
180180
from .PoleOfInaccessibility import PoleOfInaccessibility
181+
from .RasterCalculator import RasterCalculator
181182
from .CreateAttributeIndex import CreateAttributeIndex
182183
from .DropGeometry import DropGeometry
183184
from .BasicStatistics import BasicStatisticsForField
@@ -245,7 +246,7 @@ def __init__(self):
245246
RemoveNullGeometry(), ExtractByExpression(), ExtendLines(),
246247
ExtractSpecificNodes(), GeometryByExpression(), SnapGeometriesToLayer(),
247248
PoleOfInaccessibility(), CreateAttributeIndex(), DropGeometry(),
248-
BasicStatisticsForField()
249+
BasicStatisticsForField(), RasterCalculator()
249250
]
250251

251252
if hasMatplotlib:
Lines changed: 165 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,165 @@
1+
# -*- coding: utf-8 -*-
2+
3+
"""
4+
***************************************************************************
5+
RasterLayerCalculator.py
6+
---------------------
7+
Date : November 2016
8+
Copyright : (C) 2016 by Victor Olaya
9+
Email : volayaf at gmail dot com
10+
***************************************************************************
11+
* *
12+
* This program is free software; you can redistribute it and/or modify *
13+
* it under the terms of the GNU General Public License as published by *
14+
* the Free Software Foundation; either version 2 of the License, or *
15+
* (at your option) any later version. *
16+
* *
17+
***************************************************************************
18+
"""
19+
from processing.modeler.ModelerAlgorithm import ValueFromInput, ValueFromOutput
20+
import os
21+
22+
23+
__author__ = 'Victor Olaya'
24+
__date__ = 'November 2016'
25+
__copyright__ = '(C) 2016, Victor Olaya'
26+
27+
# This will get replaced with a git SHA1 when you do a git archive
28+
29+
__revision__ = '$Format:%H$'
30+
31+
import math
32+
from processing.core.GeoAlgorithm import GeoAlgorithm
33+
from processing.core.parameters import ParameterMultipleInput, ParameterExtent, ParameterString, ParameterRaster, ParameterNumber
34+
from processing.core.outputs import OutputRaster
35+
from processing.tools import dataobjects
36+
from processing.algs.gdal.GdalUtils import GdalUtils
37+
from qgis.core import QgsRectangle
38+
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
39+
from processing.core.GeoAlgorithmExecutionException import GeoAlgorithmExecutionException
40+
from processing.algs.qgis.ui.RasterCalculatorWidgets import LayersListWidgetWrapper, ExpressionWidgetWrapper
41+
42+
43+
class RasterCalculator(GeoAlgorithm):
44+
45+
LAYERS = 'LAYERS'
46+
EXTENT = 'EXTENT'
47+
CELLSIZE = 'CELLSIZE'
48+
EXPRESSION = 'EXPRESSION'
49+
OUTPUT = 'OUTPUT'
50+
51+
def defineCharacteristics(self):
52+
self.name, self.i18n_name = self.trAlgorithm('Raster calculator')
53+
self.group, self.i18n_group = self.trAlgorithm('Raster')
54+
55+
self.addParameter(ParameterMultipleInput(self.LAYERS,
56+
self.tr('Input layers'),
57+
datatype=dataobjects.TYPE_RASTER,
58+
optional=True,
59+
metadata={'widget_wrapper': LayersListWidgetWrapper}))
60+
61+
class ParameterRasterCalculatorExpression(ParameterString):
62+
63+
def evaluateForModeler(self, value, model):
64+
for i in list(model.inputs.values()):
65+
param = i.param
66+
if isinstance(param, ParameterRaster):
67+
new = "{}@".format(os.path.basename(param.value))
68+
old = "{}@".format(param.name)
69+
value = value.replace(old, new)
70+
71+
for alg in list(model.algs.values()):
72+
for out in alg.algorithm.outputs:
73+
if isinstance(out, OutputRaster):
74+
if out.value:
75+
new = "{}@".format(os.path.basename(out.value))
76+
old = "{}:{}@".format(alg.name, out.name)
77+
value = value.replace(old, new)
78+
return value
79+
80+
self.addParameter(ParameterRasterCalculatorExpression(self.EXPRESSION, self.tr('Expression'),
81+
multiline=True,
82+
metadata={'widget_wrapper': ExpressionWidgetWrapper}))
83+
self.addParameter(ParameterNumber(self.CELLSIZE,
84+
self.tr('Cellsize (use 0 or empty to set it automatically)'),
85+
minValue=0.0, default=0.0, optional=True))
86+
self.addParameter(ParameterExtent(self.EXTENT,
87+
self.tr('Output extent'),
88+
optional=True))
89+
self.addOutput(OutputRaster(self.OUTPUT, self.tr('Output')))
90+
91+
def processAlgorithm(self, progress):
92+
expression = self.getParameterValue(self.EXPRESSION)
93+
layersValue = self.getParameterValue(self.LAYERS)
94+
layersDict = {}
95+
if layersValue:
96+
layers = [dataobjects.getObjectFromUri(f) for f in layersValue.split(";")]
97+
layersDict = {os.path.basename(lyr.source()): lyr for lyr in layers}
98+
99+
for lyr in dataobjects.getRasterLayers():
100+
name = lyr.name()
101+
if (name + "@") in expression:
102+
layersDict[name] = lyr
103+
104+
entries = []
105+
for name, lyr in layersDict.items():
106+
for n in range(lyr.bandCount()):
107+
entry = QgsRasterCalculatorEntry()
108+
entry.ref = '{:s}@{:d}'.format(name, n + 1)
109+
entry.raster = lyr
110+
entry.bandNumber = n + 1
111+
entries.append(entry)
112+
113+
output = self.getOutputValue(self.OUTPUT)
114+
extentValue = self.getParameterValue(self.EXTENT)
115+
116+
if extentValue:
117+
extent = extentValue.split(',')
118+
bbox = QgsRectangle(float(extent[0]), float(extent[2]),
119+
float(extent[1]), float(extent[3]))
120+
else:
121+
if layersDict:
122+
bbox = list(layersDict.values())[0].extent()
123+
for lyr in layersDict.values():
124+
bbox.combineExtentWith(lyr.extent())
125+
else:
126+
raise GeoAlgorithmExecutionException(self.tr("No layers selected"))
127+
128+
def _cellsize(layer):
129+
return (layer.extent().xMaximum() - layer.extent().xMinimum()) / layer.width()
130+
cellsize = self.getParameterValue(self.CELLSIZE) or min([_cellsize(lyr) for lyr in layersDict.values()])
131+
width = math.floor((bbox.xMaximum() - bbox.xMinimum()) / cellsize)
132+
height = math.floor((bbox.yMaximum() - bbox.yMinimum()) / cellsize)
133+
driverName = GdalUtils.getFormatShortNameFromFilename(output)
134+
calc = QgsRasterCalculator(expression,
135+
output,
136+
driverName,
137+
bbox,
138+
width,
139+
height,
140+
entries)
141+
142+
res = calc.processCalculation()
143+
if res == QgsRasterCalculator.ParserError:
144+
raise GeoAlgorithmExecutionException(self.tr("Error parsing formula"))
145+
146+
def processBeforeAddingToModeler(self, algorithm, model):
147+
values = []
148+
expression = algorithm.params[self.EXPRESSION]
149+
for i in list(model.inputs.values()):
150+
param = i.param
151+
if isinstance(param, ParameterRaster) and "{}@".format(param.name) in expression:
152+
values.append(ValueFromInput(param.name))
153+
154+
if algorithm.name:
155+
dependent = model.getDependentAlgorithms(algorithm.name)
156+
else:
157+
dependent = []
158+
for alg in list(model.algs.values()):
159+
if alg.name not in dependent:
160+
for out in alg.algorithm.outputs:
161+
if (isinstance(out, OutputRaster)
162+
and "{}:{}@".format(alg.name, out.name) in expression):
163+
values.append(ValueFromOutput(alg.name, out.name))
164+
165+
algorithm.params[self.LAYERS] = values

0 commit comments

Comments
 (0)
Please sign in to comment.