2

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)
 ...
4

2 回答 2

5

您可以尝试 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

于 2013-04-22T07:11:13.587 回答
5

更简单的解决方案也更快：使用统一和此处解释的“方差技巧”：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 比 strides 慢 40 倍......

于 2015-11-03T15:02:19.070 回答