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:
See here for a detailed description of the answer: