Variance image in python using gdal and a running window approach

Question

Variance image in gdal

I want a local variance image with a 3x3 of a geospatial raster image using python. My approach so far was to read in the raster band as an array, then using matrix notation to run a moving window and write the array into a new raster image. This approach worked well for a high pass filter as described in this tutorial: http://www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_slides6.pdf

Then I tried to calculate the variance with several approaches, the last one using numpy (as np), but I just get a gray image with the same value everywhere. I am open to any kind of solution. If it gives me the average local variance in the end, that would be even better.

 rows = srcDS.RasterYSize
 #read in as array
 data = srcBand.ReadAsArray(0,0, cols, rows).astype(np.int)

 #calculate the variance for a 3x3 window
 outVariance = np.zeros((rows, cols), np.float)
 outVariance[1:rows-1,1:cols-1] = np.var([(data[0:rows-2,0:cols-2]),
   (data[0:rows-2,1:cols-1]),
   (data[0:rows-2,2:cols]  ),
   (data[1:rows-1,0:cols-2]),
   (data[1:rows-1,1:cols-1]),
   (data[1:rows-1,2:cols]  ),
   (data[2:rows,0:cols-2]  ),
   (data[2:rows,1:cols-1]  ),
   (data[2:rows,2:cols]    )])
 #output 
 outDS = driver.Create(outFN, cols, rows, 1, GDT_Float32)
 outDS.SetGeoTransform(srcDS.GetGeoTransform())
 outDS.SetProjection(srcDS.GetProjection())
 outBand = outDS.GetRasterBand(1)
 outBand.WriteArray(outVariance,0,0)
 ...

Answer 1

You could try Scipy, it has a function for running local filters on an array.

from scipy import ndimage

outVariance = ndimage.generic_filter(data, np.var, size=3)

It has a 'mode=' keyword for how the edges should be handled.

edit:

You can test it yourself, declare a 3x3 array:

a = np.random.rand(3,3)
a

[[ 0.01869967  0.14037373  0.32960675]
 [ 0.17213158  0.35287243  0.13498175]
 [ 0.29511881  0.46387688  0.89359801]]

For a 3x3 window, the variance of the center cell of the array will simply be:

print np.var(a)
0.058884734425985602

That value should be equal to the center cell of the returned array by Scipy:

print ndimage.generic_filter(a, np.var, size=3)
print ndimage.generic_filter(a, np.var, size=(3,3))
print ndimage.generic_filter(a, np.var, footprint=np.ones((3,3)))

[[ 0.01127325  0.01465338  0.00959321]
 [ 0.02001052  0.05888473  0.07897385]
 [ 0.00978547  0.06966683  0.09633447]]

Note that all other values in the array are 'edge-values' so the result depends on how Scipy handles the edges. It defaults to mode='reflect'.

See the documentation for more detailed information: http://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.filters.generic_filter.html

Answer 2

simpler solution and also faster : use uniform and a "variance trick" explained here : http://imagej.net/Integral_Image_Filters (the variance is the difference between "sum of square" and "square of sum")

import numpy as np
from scipy import ndimage 
rows, cols = 500, 500
win_rows, win_cols = 5, 5

img = np.random.rand(rows, cols)
win_mean = ndimage.uniform_filter(img,(win_rows,win_cols))
win_sqr_mean = ndimage.uniform_filter(img**2,(win_rows,win_cols))
win_var = win_sqr_mean - win_mean**2

the generic_filter is 40 times slower than the strides...