Skip to content
Permalink
Browse files

[processing] keep only one zonal statistics algorithm

  • Loading branch information
nirvn committed Jun 29, 2017
1 parent 58f6f93 commit 26d9c74bc40e02b5fae3215b0314913a4bd1ec9b
@@ -72,7 +72,7 @@
from .SpatialiteExecuteSQL import SpatialiteExecuteSQL
from .SymmetricalDifference import SymmetricalDifference
from .VectorSplit import VectorSplit
from .ZonalStatisticsQgis import ZonalStatisticsQgis
from .ZonalStatistics import ZonalStatistics

# from .ExtractByLocation import ExtractByLocation
# from .PointsInPolygon import PointsInPolygon
@@ -119,7 +119,6 @@
# from .JoinAttributes import JoinAttributes
# from .CreateConstantRaster import CreateConstantRaster
# from .PointsDisplacement import PointsDisplacement
# from .ZonalStatistics import ZonalStatistics
# from .PointsFromPolygons import PointsFromPolygons
# from .PointsFromLines import PointsFromLines
# from .RandomPointsExtent import RandomPointsExtent
@@ -210,7 +209,7 @@ def getAlgs(self):
# EquivalentNumField(),
# StatisticsByCategories(), ConcaveHull(),
# RasterLayerStatistics(), PointsDisplacement(),
# ZonalStatistics(), PointsFromPolygons(),
# PointsFromPolygons(),
# PointsFromLines(), RandomPointsExtent(),
# RandomPointsLayer(), RandomPointsPolygonsFixed(),
# RandomPointsPolygonsVariable(),
@@ -273,7 +272,7 @@ def getAlgs(self):
SpatialiteExecuteSQL(),
SymmetricalDifference(),
VectorSplit(),
ZonalStatisticsQgis()
ZonalStatistics()
]

if hasPlotly:
@@ -4,8 +4,8 @@
***************************************************************************
ZonalStatistics.py
---------------------
Date : August 2013
Copyright : (C) 2013 by Alexander Bruy
Date : September 2016
Copyright : (C) 2016 by Alexander Bruy
Email : alexander dot bruy at gmail dot com
***************************************************************************
* *
@@ -16,41 +16,33 @@
* *
***************************************************************************
"""
from builtins import str

__author__ = 'Alexander Bruy'
__date__ = 'August 2013'
__copyright__ = '(C) 2013, Alexander Bruy'
__date__ = 'September 2016'
__copyright__ = '(C) 2016, Alexander Bruy'

# This will get replaced with a git SHA1 when you do a git archive

__revision__ = '$Format:%H$'

import numpy
import os

try:
from scipy.stats.mstats import mode
hasSciPy = True
except:
hasSciPy = False
from qgis.PyQt.QtGui import QIcon

from osgeo import gdal, ogr, osr
from qgis.core import (QgsApplication,
QgsFeatureSink,
QgsRectangle,
QgsGeometry,
QgsFeature,
QgsProcessingUtils)
from qgis.analysis import QgsZonalStatistics
from qgis.core import (QgsFeatureSink,
QgsProcessingUtils,
QgsProcessingParameterDefinition,
QgsProcessingParameterVectorLayer,
QgsProcessingParameterRasterLayer,
QgsProcessingParameterString,
QgsProcessingParameterNumber,
QgsProcessingParameterEnum,
QgsProcessingOutputVectorLayer)

from processing.algs.qgis.QgisAlgorithm import QgisAlgorithm
from processing.core.parameters import ParameterVector
from processing.core.parameters import ParameterRaster
from processing.core.parameters import ParameterString
from processing.core.parameters import ParameterNumber
from processing.core.parameters import ParameterBoolean
from processing.core.outputs import OutputVector
from processing.tools.raster import mapToPixel
from processing.tools import dataobjects, vector

pluginPath = os.path.split(os.path.split(os.path.dirname(__file__))[0])[0]


class ZonalStatistics(QgisAlgorithm):
@@ -59,26 +51,48 @@ class ZonalStatistics(QgisAlgorithm):
RASTER_BAND = 'RASTER_BAND'
INPUT_VECTOR = 'INPUT_VECTOR'
COLUMN_PREFIX = 'COLUMN_PREFIX'
GLOBAL_EXTENT = 'GLOBAL_EXTENT'
OUTPUT_LAYER = 'OUTPUT_LAYER'
STATISTICS = 'STATS'

def icon(self):
return QIcon(os.path.join(pluginPath, 'images', 'zonalstats.png'))

def group(self):
return self.tr('Raster tools')

def __init__(self):
super().__init__()
self.addParameter(ParameterRaster(self.INPUT_RASTER,
self.tr('Raster layer')))
self.addParameter(ParameterNumber(self.RASTER_BAND,
self.tr('Raster band'), 1, 999, 1))
self.addParameter(ParameterVector(self.INPUT_VECTOR,
self.tr('Vector layer containing zones'),
[dataobjects.TYPE_VECTOR_POLYGON]))
self.addParameter(ParameterString(self.COLUMN_PREFIX,
self.tr('Output column prefix'), '_'))
self.addParameter(ParameterBoolean(self.GLOBAL_EXTENT,
self.tr('Load whole raster in memory')))
self.addOutput(OutputVector(self.OUTPUT_LAYER, self.tr('Zonal statistics'), datatype=[dataobjects.TYPE_VECTOR_POLYGON]))
self.STATS = {self.tr('Count'): QgsZonalStatistics.Count,
self.tr('Sum'): QgsZonalStatistics.Sum,
self.tr('Mean'): QgsZonalStatistics.Mean,
self.tr('Median'): QgsZonalStatistics.Median,
self.tr('Std. dev.'): QgsZonalStatistics.StDev,
self.tr('Min'): QgsZonalStatistics.Min,
self.tr('Max'): QgsZonalStatistics.Max,
self.tr('Range'): QgsZonalStatistics.Range,
self.tr('Minority'): QgsZonalStatistics.Minority,
self.tr('Majority (mode)'): QgsZonalStatistics.Majority,
self.tr('Variety'): QgsZonalStatistics.Variety,
self.tr('Variance'): QgsZonalStatistics.Variance,
self.tr('All'): QgsZonalStatistics.All
}

self.addParameter(QgsProcessingParameterRasterLayer(self.INPUT_RASTER,
self.tr('Raster layer')))
self.addParameter(QgsProcessingParameterNumber(self.RASTER_BAND,
self.tr('Raster band'),
minValue=1, maxValue=999, defaultValue=1))
self.addParameter(QgsProcessingParameterVectorLayer(self.INPUT_VECTOR,
self.tr('Vector layer containing zones'),
[QgsProcessingParameterDefinition.TypeVectorPolygon]))
self.addParameter(QgsProcessingParameterString(self.COLUMN_PREFIX,
self.tr('Output column prefix'), '_'))
self.addParameter(QgsProcessingParameterEnum(self.STATISTICS,
self.tr('Statistics to calculate'),
list(self.STATS.keys()),
allowMultiple=True))
self.addOutput(QgsProcessingOutputVectorLayer(self.INPUT_VECTOR,
self.tr('Zonal statistics'),
QgsProcessingParameterDefinition.TypeVectorPolygon))

def name(self):
return 'zonalstatistics'
@@ -87,191 +101,23 @@ def displayName(self):
return self.tr('Zonal Statistics')

def processAlgorithm(self, parameters, context, feedback):
""" Based on code by Matthew Perry
https://gist.github.com/perrygeo/5667173
:param parameters:
:param context:
"""

layer = QgsProcessingUtils.mapLayerFromString(self.getParameterValue(self.INPUT_VECTOR), context)

rasterPath = str(self.getParameterValue(self.INPUT_RASTER))
bandNumber = self.getParameterValue(self.RASTER_BAND)
columnPrefix = self.getParameterValue(self.COLUMN_PREFIX)
useGlobalExtent = self.getParameterValue(self.GLOBAL_EXTENT)

rasterDS = gdal.Open(rasterPath, gdal.GA_ReadOnly)
geoTransform = rasterDS.GetGeoTransform()
rasterBand = rasterDS.GetRasterBand(bandNumber)
noData = rasterBand.GetNoDataValue()

cellXSize = abs(geoTransform[1])
cellYSize = abs(geoTransform[5])
rasterXSize = rasterDS.RasterXSize
rasterYSize = rasterDS.RasterYSize

rasterBBox = QgsRectangle(geoTransform[0],
geoTransform[3] - cellYSize * rasterYSize,
geoTransform[0] + cellXSize * rasterXSize,
geoTransform[3])

rasterGeom = QgsGeometry.fromRect(rasterBBox)

crs = osr.SpatialReference()
crs.ImportFromProj4(str(layer.crs().toProj4()))

if useGlobalExtent:
xMin = rasterBBox.xMinimum()
xMax = rasterBBox.xMaximum()
yMin = rasterBBox.yMinimum()
yMax = rasterBBox.yMaximum()

(startColumn, startRow) = mapToPixel(xMin, yMax, geoTransform)
(endColumn, endRow) = mapToPixel(xMax, yMin, geoTransform)

width = endColumn - startColumn
height = endRow - startRow

srcOffset = (startColumn, startRow, width, height)
srcArray = rasterBand.ReadAsArray(*srcOffset)
srcArray = srcArray * rasterBand.GetScale() + rasterBand.GetOffset()

newGeoTransform = (
geoTransform[0] + srcOffset[0] * geoTransform[1],
geoTransform[1],
0.0,
geoTransform[3] + srcOffset[1] * geoTransform[5],
0.0,
geoTransform[5],
)

memVectorDriver = ogr.GetDriverByName('Memory')
memRasterDriver = gdal.GetDriverByName('MEM')

fields = layer.fields()
(idxMin, fields) = vector.findOrCreateField(layer, fields,
columnPrefix + 'min', 21, 6)
(idxMax, fields) = vector.findOrCreateField(layer, fields,
columnPrefix + 'max', 21, 6)
(idxSum, fields) = vector.findOrCreateField(layer, fields,
columnPrefix + 'sum', 21, 6)
(idxCount, fields) = vector.findOrCreateField(layer, fields,
columnPrefix + 'count', 21, 6)
(idxMean, fields) = vector.findOrCreateField(layer, fields,
columnPrefix + 'mean', 21, 6)
(idxStd, fields) = vector.findOrCreateField(layer, fields,
columnPrefix + 'std', 21, 6)
(idxUnique, fields) = vector.findOrCreateField(layer, fields,
columnPrefix + 'unique', 21, 6)
(idxRange, fields) = vector.findOrCreateField(layer, fields,
columnPrefix + 'range', 21, 6)
(idxVar, fields) = vector.findOrCreateField(layer, fields,
columnPrefix + 'var', 21, 6)
(idxMedian, fields) = vector.findOrCreateField(layer, fields,
columnPrefix + 'median', 21, 6)
if hasSciPy:
(idxMode, fields) = vector.findOrCreateField(layer, fields,
columnPrefix + 'mode', 21, 6)

writer = self.getOutputFromName(self.OUTPUT_LAYER).getVectorWriter(fields, layer.wkbType(),
layer.crs(), context)

outFeat = QgsFeature()

outFeat.initAttributes(len(fields))
outFeat.setFields(fields)

features = QgsProcessingUtils.getFeatures(layer, context)
total = 100.0 / layer.featureCount() if layer.featureCount() else 0
for current, f in enumerate(features):
geom = f.geometry()

intersectedGeom = rasterGeom.intersection(geom)
ogrGeom = ogr.CreateGeometryFromWkt(intersectedGeom.exportToWkt())

if not useGlobalExtent:
bbox = intersectedGeom.boundingBox()

xMin = bbox.xMinimum()
xMax = bbox.xMaximum()
yMin = bbox.yMinimum()
yMax = bbox.yMaximum()

(startColumn, startRow) = mapToPixel(xMin, yMax, geoTransform)
(endColumn, endRow) = mapToPixel(xMax, yMin, geoTransform)

width = endColumn - startColumn
height = endRow - startRow

if width == 0 or height == 0:
continue

srcOffset = (startColumn, startRow, width, height)
srcArray = rasterBand.ReadAsArray(*srcOffset)
srcArray = srcArray * rasterBand.GetScale() + rasterBand.GetOffset()

newGeoTransform = (
geoTransform[0] + srcOffset[0] * geoTransform[1],
geoTransform[1],
0.0,
geoTransform[3] + srcOffset[1] * geoTransform[5],
0.0,
geoTransform[5],
)

# Create a temporary vector layer in memory
memVDS = memVectorDriver.CreateDataSource('out')
memLayer = memVDS.CreateLayer('poly', crs, ogr.wkbPolygon)

ft = ogr.Feature(memLayer.GetLayerDefn())
ft.SetGeometry(ogrGeom)
memLayer.CreateFeature(ft)
ft.Destroy()

# Rasterize it
rasterizedDS = memRasterDriver.Create('', srcOffset[2],
srcOffset[3], 1, gdal.GDT_Byte)
rasterizedDS.SetGeoTransform(newGeoTransform)
gdal.RasterizeLayer(rasterizedDS, [1], memLayer, burn_values=[1])
rasterizedArray = rasterizedDS.ReadAsArray()

srcArray = numpy.nan_to_num(srcArray)
masked = numpy.ma.MaskedArray(srcArray,
mask=numpy.logical_or(srcArray == noData,
numpy.logical_not(rasterizedArray)))

outFeat.setGeometry(geom)

attrs = f.attributes()
v = float(masked.min())
attrs.insert(idxMin, None if numpy.isnan(v) else v)
v = float(masked.max())
attrs.insert(idxMax, None if numpy.isnan(v) else v)
v = float(masked.sum())
attrs.insert(idxSum, None if numpy.isnan(v) else v)
attrs.insert(idxCount, int(masked.count()))
v = float(masked.mean())
attrs.insert(idxMean, None if numpy.isnan(v) else v)
v = float(masked.std())
attrs.insert(idxStd, None if numpy.isnan(v) else v)
attrs.insert(idxUnique, numpy.unique(masked.compressed()).size)
v = float(masked.max()) - float(masked.min())
attrs.insert(idxRange, None if numpy.isnan(v) else v)
v = float(masked.var())
attrs.insert(idxVar, None if numpy.isnan(v) else v)
v = float(numpy.ma.median(masked))
attrs.insert(idxMedian, None if numpy.isnan(v) else v)
if hasSciPy:
attrs.insert(idxMode, float(mode(masked, axis=None)[0][0]))

outFeat.setAttributes(attrs)
writer.addFeature(outFeat, QgsFeatureSink.FastInsert)

memVDS = None
rasterizedDS = None

feedback.setProgress(int(current * total))

rasterDS = None

del writer
bandNumber = self.parameterAsInt(parameters, self.RASTER_BAND, context)
columnPrefix = self.parameterAsString(parameters, self.COLUMN_PREFIX, context)
st = self.parameterAsEnums(parameters, self.STATISTICS, context)

vectorLayer = self.parameterAsVectorLayer(parameters, self.INPUT_VECTOR, context)
rasterLayer = self.parameterAsRasterLayer(parameters, self.INPUT_RASTER, context)

keys = list(self.STATS.keys())
selectedStats = 0
for i in st:
selectedStats |= self.STATS[keys[i]]

zs = QgsZonalStatistics(vectorLayer,
rasterLayer,
columnPrefix,
bandNumber,
selectedStats)
zs.calculateStatistics(feedback)

return {self.INPUT_VECTOR: vectorLayer}

0 comments on commit 26d9c74

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