Skip to content
Permalink
Browse files

[processing] add zonal statistics algorithm

  • Loading branch information
alexbruy committed Sep 13, 2013
1 parent 72536f6 commit 8481ebcc9181d56778263a3ffdc6101a9d273fd8
@@ -87,6 +87,7 @@
from processing.algs.PointsLayerFromTable import PointsLayerFromTable

from processing.algs.PointsDisplacement import PointsDisplacement
from sextante.algs.ZonalStatistics import ZonalStatistics

#from processing.algs.VectorLayerHistogram import VectorLayerHistogram
#from processing.algs.VectorLayerScatterplot import VectorLayerScatterplot
@@ -141,7 +142,8 @@ def __init__(self):
#VectorLayerHistogram(), VectorLayerScatterplot(), RasterLayerHistogram(),
#MeanAndStdDevPlot(), BarPlot(), PolarPlot()
# ------ vector ------
PointsDisplacement()
PointsDisplacement(),
ZonalStatistics()
]

def initializeSettings(self):
@@ -0,0 +1,98 @@
# -*- coding: utf-8 -*-

"""
***************************************************************************
QGISUtils.py
---------------------
Date : August 2013
Copyright : (C) 2013 by Alexander Bruy
Email : alexander dot bruy 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__ = 'Alexander Bruy'
__date__ = 'August 2013'
__copyright__ = '(C) 2013, Alexander Bruy'
# This will get replaced with a git SHA1 when you do a git archive
__revision__ = '$Format:%H$'

import math

from PyQt4.QtCore import *

from qgis.core import *


def mapToPixel(mX, mY, geoTransform):
'''Convert map coordinates to pixel coordinates.
@param mX Input map X coordinate (double)
@param mY Input map Y coordinate (double)
@param geoTransform Input geotransform (six doubles)
@return pX, pY Output coordinates (two doubles)
'''
if geoTransform[2] + geoTransform[4] == 0:
pX = (mX - geoTransform[0]) / geoTransform[1]
pY = (mY - geoTransform[3]) / geoTransform[5]
else:
pX, pY = applyGeoTransform(mX, mY, invertGeoTransform(geoTransform))
return int(pX), int(pY)

def pixelToMap(pX, pY, geoTransform):
'''Convert pixel coordinates to map coordinates.
@param pX Input pixel X coordinate (double)
@param pY Input pixel Y coordinate (double)
@param geoTransform Input geotransform (six doubles)
@return mX, mY Output coordinates (two doubles)
'''
mX, mY = applyGeoTransform(pX + 0.5, pY + 0.5, geoTransform)
return mX, mY

def applyGeoTransform(inX, inY, geoTransform):
'''Apply a geotransform to coordinates.
@param inX Input coordinate (double)
@param inY Input coordinate (double)
@param geoTransform Input geotransform (six doubles)
@return outX, outY Output coordinates (two doubles)
'''
outX = geoTransform[0] + inX * geoTransform[1] + inY * geoTransform[2]
outY = geoTransform[3] + inX * geoTransform[4] + inY * geoTransform[5]
return outX, outY

def invertGeoTransform(geoTransform):
'''Invert standard 3x2 set of geotransform coefficients.
@param geoTransform Input GeoTransform (six doubles - unaltered)
@return outGeoTransform Output GeoTransform (six doubles - updated)
on success, None if the equation is uninvertable
'''
# we assume a 3rd row that is [1 0 0]
# compute determinate
det = geoTransform[1] * geoTransform[5] - geoTransform[2] * geoTransform[4]

if abs(det) < 0.000000000000001:
return

invDet = 1.0 / det

# compute adjoint and divide by determinate
outGeoTransform = [0, 0, 0, 0, 0, 0]
outGeoTransform[1] = geoTransform[5] * invDet
outGeoTransform[4] = -geoTransform[4] * invDet

outGeoTransform[2] = -geoTransform[2] * invDet
outGeoTransfrom[5] = geoTransform[1] * invDet

outGeoTransform[0] = (geoTransform[2] * geoTransform[3] - geoTransform[0] * geoTransform[5]) * invDet
outGeoTransform[3] = (-geoTransform[1] * geoTransform[3] + geoTransform[0] * geoTransform[4]) * invDet

return outGeoTransform
@@ -0,0 +1,226 @@
# -*- coding: utf-8 -*-

"""
***************************************************************************
ZonalStatistics.py
---------------------
Date : August 2013
Copyright : (C) 2013 by Alexander Bruy
Email : alexander dot bruy 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__ = 'Alexander Bruy'
__date__ = 'August 2013'
__copyright__ = '(C) 2013, Alexander Bruy'
# This will get replaced with a git SHA1 when you do a git archive
__revision__ = '$Format:%H$'

from PyQt4.QtCore import *

import numpy
from osgeo import gdal, ogr, osr

from qgis.core import *

from sextante.core.GeoAlgorithm import GeoAlgorithm
from sextante.core.QGisLayers import QGisLayers

from sextante.parameters.ParameterVector import ParameterVector
from sextante.parameters.ParameterRaster import ParameterRaster
from sextante.parameters.ParameterString import ParameterString
from sextante.parameters.ParameterNumber import ParameterNumber
from sextante.parameters.ParameterBoolean import ParameterBoolean

from sextante.outputs.OutputVector import OutputVector

from sextante.algs.ftools import FToolsUtils as ftools_utils
from sextante.algs import QGISUtils as utils

class ZonalStatistics(GeoAlgorithm):

INPUT_RASTER = "INPUT_RASTER"
RASTER_BAND = "RASTER_BAND"
INPUT_VECTOR = "INPUT_VECTOR"
COLUMN_PREFIX = "COLUMN_PREFIX"
GLOBAL_EXTENT = "GLOBAL_EXTENT"
OUTPUT_LAYER = "OUTPUT_LAYER"

def defineCharacteristics(self):
self.name = "Zonal Statistics"
self.group = "Raster tools"

self.addParameter(ParameterRaster(self.INPUT_RASTER, "Raster layer"))
self.addParameter(ParameterNumber(self.RASTER_BAND, "Raster band", 1, 999, 1))
self.addParameter(ParameterVector(self.INPUT_VECTOR, "Vector layer containing zones", [ParameterVector.VECTOR_TYPE_POLYGON]))
self.addParameter(ParameterString(self.COLUMN_PREFIX, "Output column prefix", "_"))
self.addParameter(ParameterBoolean(self.GLOBAL_EXTENT, "Load whole raster in memory"))
self.addOutput(OutputVector(self.OUTPUT_LAYER, "Output layer"))

def processAlgorithm(self, progress):
layer = QGisLayers.getObjectFromUri(self.getParameterValue(self.INPUT_VECTOR))

rasterPath = unicode(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 = utils.mapToPixel(xMin, yMax, geoTransform)
endColumn, endRow = utils.mapToPixel(xMax, yMin, geoTransform)

width = endColumn - startColumn
height = endRow - startRow

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

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.pendingFields()
idxMin, fields = ftools_utils.findOrCreateField(layer, fields, columnPrefix + "min", 21, 6)
idxMax, fields = ftools_utils.findOrCreateField(layer, fields, columnPrefix + "max", 21, 6)
idxSum, fields = ftools_utils.findOrCreateField(layer, fields, columnPrefix + "sum", 21, 6)
idxCount, fields = ftools_utils.findOrCreateField(layer, fields, columnPrefix + "count", 21, 6)
idxMean, fields = ftools_utils.findOrCreateField(layer, fields, columnPrefix + "mean", 21, 6)
idxStd, fields = ftools_utils.findOrCreateField(layer, fields, columnPrefix + "std", 21, 6)
idxUnique, fields = ftools_utils.findOrCreateField(layer, fields, columnPrefix + "unique", 21, 6)
idxRange, fields = ftools_utils.findOrCreateField(layer, fields, columnPrefix + "range", 21, 6)
idxCV, fields = ftools_utils.findOrCreateField(layer, fields, columnPrefix + "cv", 21, 6)
#idxMedian, fields = ftools_utils.findOrCreateField(layer, fields, columnPrefix + "median", 21, 6)

writer = self.getOutputFromName(self.OUTPUT_LAYER).getVectorWriter(fields.toList(),
layer.dataProvider().geometryType(), layer.crs())

outFeat = QgsFeature()

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

current = 0
features = QGisLayers.features(layer)
total = 100.0 / len(features)
for f in 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 = utils.mapToPixel(xMin, yMax, geoTransform)
endColumn, endRow = utils.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)

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()
attrs.insert(idxMin, float(masked.min()))
attrs.insert(idxMax, float(masked.max()))
attrs.insert(idxSum, float(masked.sum()))
attrs.insert(idxCount, int(masked.count()))
attrs.insert(idxMean, float(masked.mean()))
attrs.insert(idxStd, float(masked.std()))
attrs.insert(idxUnique, numpy.unique(masked.compressed()).size)
attrs.insert(idxRange, float(masked.max()) - float(masked.min()))
attrs.insert(idxCV, float(masked.var()))
#attrs.insert(idxMedian, float(masked.median()))

outFeat.setAttributes(attrs)
writer.addFeature(outFeat)

memVDS = None
rasterizedDS = None

current += 1
progress.setPercentage(int(current * total))

rasterDS = None

del writer

0 comments on commit 8481ebc

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