使用gdal和运行窗口方法的python中的方差图像

Question

gdal中的方差图像

我想要一个使用python的3x3地理空间光栅图像的局部方差图像。到目前为止，我的方法是将光栅带作为一个数组读取，然后使用矩阵表示法运行移动窗口，并将该数组写入新的光栅图像。此方法适用于本教程中描述的高通滤波器：http://www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_slides6.pdf

然后，我尝试用几种方法计算方差，最后一种方法使用numpy (作为np)，但我只得到了一幅到处都是相同值的灰度图像。我对任何类型的解决方案都持开放态度。如果它最终给出了平均局部方差，那就更好了。

 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

你可以试试Scipy，它有一个在数组上运行本地过滤器的函数。

from scipy import ndimage

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

它有一个'mode=‘关键字，表示应该如何处理边缘。

编辑：

您可以自己测试它，声明一个3x3数组：

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]]

对于3x3窗口，阵列中心单元的方差将简单地为：

print np.var(a)
0.058884734425985602

该值应等于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]]

请注意，数组中的所有其他值都是“边值”，因此结果取决于Scipy如何处理边。默认为mode='reflect'。

有关更详细的信息，请参阅文档：http://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.filters.generic_filter.html

Answer 2

更简单的解决方案，也更快:使用统一的和这里解释的“方差技巧”：http://imagej.net/Integral_Image_Filters (方差是“平方和”和“和的平方”之间的差)

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

generic_filter比步幅慢40倍...