[GIS] Local variance image in python using gdal and a running window approach

filtergdalpython

I want a local variance image with a 3×3 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 Thanks, Mace

  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)
  ...

Best Answer

This trick from another user on Stackoverflow solved the problem:

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

See here for a detailed description of the answer: