I am trying to find local maxima in a DEM (canopy height model) that represents tree tops. I can do this with a fixed window size. However I would rather have a variable window size based on the height of each cell with a function such as windowSize=x*0.1+1. I found this approach in a in R package ForestTools, based got the method from a paper by Popescu and Wynne from 2004.
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 but so far without success. Ideally I would like to only use Python. The Forest tools has some issues with large files, I could also move to R if it has a better alternative.
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. I have tried fitting a simple function there such as:
def varvalue(x): y = x*0.1+1 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.
Does anyone have an idea or hint how to approach this problem?