title: Variance image in python using gdal and a running window approach
body: <p>Variance image in gdal</p>
  
  <p>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: <a href="http://www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_slides6.pdf" rel="nofollow">http://www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_slides6.pdf</a></p>
  
  <p>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.</p>
  
  <pre><code> 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)
   ...
  </code></pre>