gpt4 book ai didi

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

转载 作者:行者123 更新时间:2023-11-28 20:25:07 27 4
gpt4 key购买 nike

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

最佳答案

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

关于python - 使用 gdal 和运行窗口方法的 python 中的方差图像,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/16107671/

27 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com