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,它有一个在数组上运行本地过滤器的函数。在
它有一个“mode=”关键字来表示边缘应该如何处理。在
编辑:
您可以自己测试,声明一个3x3数组:
^{pr2}$对于3x3窗口,数组中心单元格的方差将简单地为:
该值应等于Scipy返回的数组的中心单元格:
请注意,数组中的所有其他值都是“边值”,因此结果取决于Scipy如何处理边。它默认为
mode='reflect'
。在有关详细信息,请参阅文档: http://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.filters.generic_filter.html
更简单、更快的解决方案:使用统一的 这里解释了一个“方差技巧”:http://imagej.net/Integral_Image_Filters(方差是“平方和”和“和的平方之差”)
通用的过滤器比步伐慢40倍。。。在
相关问题 更多 >
编程相关推荐