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

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数组：

^{pr2}$

对于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