Skip to content

Commit

Permalink
[Processing] Fixes for zonal statistics algorithm:
Browse files Browse the repository at this point in the history
1. Mask NaN values instead of converting them to 0.
2. Handle rasters for which raster value offset and scale have not been
set.
  • Loading branch information
radosuav committed Jun 30, 2017
1 parent 6e52f1b commit 4886b36
Showing 1 changed file with 12 additions and 5 deletions.
17 changes: 12 additions & 5 deletions python/plugins/processing/algs/qgis/ZonalStatistics.py
Original file line number Diff line number Diff line change
Expand Up @@ -89,6 +89,12 @@ def processAlgorithm(self, progress):
geoTransform = rasterDS.GetGeoTransform()
rasterBand = rasterDS.GetRasterBand(bandNumber)
noData = rasterBand.GetNoDataValue()
scale = rasterBand.GetScale()
if scale is None:
scale = 1.0
offset = rasterBand.GetOffset()
if offset is None:
offset = 0.0

cellXSize = abs(geoTransform[1])
cellYSize = abs(geoTransform[5])
Expand Down Expand Up @@ -118,7 +124,7 @@ def processAlgorithm(self, progress):

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

newGeoTransform = (
geoTransform[0] + srcOffset[0] * geoTransform[1],
Expand Down Expand Up @@ -192,7 +198,7 @@ def processAlgorithm(self, progress):

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

newGeoTransform = (
geoTransform[0] + srcOffset[0] * geoTransform[1],
Expand All @@ -219,10 +225,11 @@ def processAlgorithm(self, progress):
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)))
mask=numpy.logical_or.reduce((
srcArray == noData,
numpy.logical_not(rasterizedArray),
numpy.isnan(srcArray))))

outFeat.setGeometry(geom)

Expand Down

0 comments on commit 4886b36

Please sign in to comment.