I am trying to find local maxima in a 2d array. I can do this with a fixed window size. However I would rather have a variable window size based on the value of each cell with a function such as window size=x*0.1+2. I found this approach in a in R package ForestTools that got the method from a paper by Popescu and Wynne from 2004. So basically, if a cell has a value of 10 the local maxima should be looked for in a window size of 3 (10*0.1+2=3) and if the value is 20 the window size should be 4 (20*0.1+2=4)
At the moment I am trying to do the same in Python using the function peak_local_max from the Scikit package but this does not seem to have any options for variable windows. I have looked for other python packages/functions, such as scipy.signal, but so far without success. Ideally I would like to stay within Python as it is quite fast. As languages I only know Python and R.
Code so far:
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)))
My code reads a raster as an array to find the local maxima. The parameter footprint gives the search window in this case. The function peak_local_max has two parameters that can define the window size called: size and footprint. I have tried fitting a simple function there such as:
def varvalue(x):
y = x*0.1+2
return y
However I am not sure how to pass the cell value to this function and I wonder if this is even possible. I have looked at using other functions/packages but so far without any result. Most Python related answers I find are based on one dimensional data sets and GIS answers use QGIS or Arcmap.
Extra info: I am doing this to find the tops of trees. It makes sense that larger trees need a bigger window because the chances of finding two tall trees within a short distance is small. That is why my data is read in as a raster and converted to an array. Playing with the different footprints gives different results. With a fixed window size there will either be an over- or underestimation of the amount of trees. I used this tutorial as reference.