我试图在二维数组中找到局部最大值。我可以用固定的窗口大小来做到这一点。但是,我宁愿根据每个单元格的值使用窗口大小=x*0.1+2 之类的函数来设置可变窗口大小。我在 R 包ForestTools中找到了这种方法,该方法从Popescu 和 Wynne 2004 年的一篇论文中获得。所以基本上,如果一个单元格的值为 10,则应在窗口大小为 3(10*0.1+2=3)的窗口中查找局部最大值,如果值为 20,则窗口大小应为 4(20*0.1+ 2=4)
目前我正在尝试使用 Scikit 包中的函数 peak_local_max 在 Python 中做同样的事情,但这似乎没有任何变量窗口的选项。我一直在寻找其他 python 包/函数,例如 scipy.signal,但到目前为止没有成功。理想情况下,我想留在 Python 中,因为它非常快。作为语言,我只知道 Python 和 R。
到目前为止的代码:
chm_file = "my_chm.tif"
chm_dataset = gdal.Open(my_chm)
chm_raster = chm_dataset.GetRasterBand(1)
cols_chm = chm_dataset.RasterXSize
rows_chm = chm_dataset.RasterYSize
chm_array = chm_raster.ReadAsArray(0,0,cols_chm,rows_chm).astype(np.float)
local_maxi = peak_local_max(chm_array,indices=False, footprint=np.ones((10, 10)))
我的代码将栅格读取为数组以查找局部最大值。在这种情况下,参数footprint 给出了搜索窗口。函数 peak_local_max 有两个可以定义窗口大小的参数:size 和footprint。我尝试在那里安装一个简单的功能,例如:
def varvalue(x):
y = x*0.1+2
return y
但是我不确定如何将单元格值传递给这个函数,我想知道这是否可能。我已经研究过使用其他功能/包,但到目前为止没有任何结果。我发现的大多数 Python 相关答案都基于一维数据集,而 GIS 答案使用 QGIS 或 Arcmap。
额外信息:我这样做是为了找到树顶。较大的树木需要较大的窗户是有道理的,因为在短距离内找到两棵高大树木的机会很小。这就是为什么我的数据作为栅格读入并转换为数组的原因。玩不同的脚印会产生不同的结果。在固定窗口大小的情况下,树木的数量要么被高估,要么被低估。我使用本教程作为参考。