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)))
