[Processing] Fixes for zonal statistics algorithm:

radosuav · radosuav · commit 4886b36bce39 · 2017-06-30T17:13:02.000+02:00
1. Mask NaN values instead of converting them to 0.
2. Handle rasters for which raster value offset and scale have not been
set.
diff --git a/python/plugins/processing/algs/qgis/ZonalStatistics.py b/python/plugins/processing/algs/qgis/ZonalStatistics.py
@@ -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])
@@ -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],
@@ -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],
@@ -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)
 
