Spend a good portion of yesterday exploring code for finding peaks in a digital elevation model. I tried a few different methods, one searching for highest points within the x and y axis, looking for inclines and declines, and picking the points in the middle, but that found significantly more points then desired as it seems there is a lot of small up and downs in the elevation, even if they aren’t the actual summit.

What I ended up settling on was this code, that first finds the highest point in the DEM, then deletes all points within 100 meters in all directions from a copy of the DEM then recurvisely calls the function. This works fairly well, although it is memory intensive after about 10 calls, as you are making a lot of copies of DEM in memory. But it seems to work well, assuming you set a reasonable size for the buffer based on the size of the mountain.

library(tidyverse)
library(terra)
library(units)
library(sf)
rm(list=ls())

dem <- rast('/tmp/merge_5925072950_2_meter.tif')
demdf <- as.data.frame(dem, xy=T)  %>%
  mutate(Layer_1 = Layer_1 %>% set_units('m') %>% set_units('ft')) %>% drop_units()

dems <- aggregate(dem, fact=4)
demsdf <- as.data.frame(dems, xy=T) %>%
  mutate(Layer_1 = Layer_1 %>% set_units('m') %>% set_units('ft')) %>% drop_units()


findPeaks <- function(df, size=100, n=5, peakdf=NA) {
  dfp <- df
  newPeak <- dfp %>% slice_max(Layer_1)
  
  ifelse(is.na(peakdf),  peakdf <- newPeak, peakdf <- rbind(peakdf, newPeak))
  
  dfp <- 
    dfp %>% 
    filter( !(between(x, newPeak[1,1]-size, newPeak[1,1]+size) &
                between(y, newPeak[1,2]-size, newPeak[1,2]+size)))
            
  if (nrow(peakdf) >= n) {
    return(peakdf)                          
  }
  else {
    findPeaks(dfp, size, n, peakdf)
  }
}


findPeaks(demdf, 300, 10) %>% 
  st_as_sf(coords=c('x', 'y'), crs=26918,  agr = "constant") %>% 
  ggplot() +
  geom_raster(data=demsdf, aes(x=x, y=y, fill=Layer_1)) +
  geom_sf(color='white') +
  ggsflabel::geom_sf_label_repel(aes(label=floor(Layer_1)))