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

Question

gdal中的方差图像

我想要一个局部方差图像与3x3的地理空间光栅图像使用python.到目前为止,我的方法是将光栅带作为数组读入,然后使用矩阵表示法运行移动窗口并将数组写入新的光栅图像.这种方法适用于本教程中描述的高通滤波器: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 倍......