# How to calculate peaks in a DEM

How could I extract all peaks from a DEM (large mountainous area)? Peak means something like: a pixel which is sourrounded by other pixels with smaler elevation numbers in a area arround the "peak-pixel" of about 100m.

• "Extract" to what format? – whuber Dec 7 '13 at 20:22
• any format. Vector, Raster, coordinates – MAP Dec 7 '13 at 20:24
• and what defines an "area"? Otherwise you might end up with ether a single or hundreds of peaks in the landscape – Curlew Dec 7 '13 at 20:24
• This question is clear--and it answers itself: it asks to compare the DEM to a 100m focal maximum of the same DEM. The result will be an indicator (true/false) grid showing the peaks. – whuber Dec 7 '13 at 20:27

Use the grass module r.param.scale with the "param=feature" option. THis creates an output map with each pixel categorized as peak, ridge, channel, plane, etc. Peaks are given category 6 (ridges=5, etc)

• Where as this does return a landoform classification it is not under the criteria that was defined in the post. – Jeffrey Evans Dec 8 '13 at 16:55
• ...but I guess I can work with it! Still testing... – MAP Dec 11 '13 at 8:38
• Finaly I can get peaks from a DEM but I have to do a lot of work after using r.param.scale. The WOODS-Algorythm generates a peak area, and not a discrete peak. I have to vectorize to polygone and find its centroid. But the centroid is not always the pixel with the highes elevation. Im working on a python-script doing the work... – MAP Dec 16 '13 at 20:00

A pure pixel based approach could be done with "r.mapcalc" ([..]) but will likely not lead to exciting results. Hence r.param.scale as suggested by Micha or this Addon: r.prominence which calculates the average difference between a central cell and its neighbors. It approximated the terrain 'ruggedness' by looking at average differences in elevation within a given neighborhood. See http://grasswiki.osgeo.org/wiki/AddOns/GRASS_6#r.prominence

I took the reverse approach to finding peaks, I found the peak pixels first, and am now trying to select those exceeding a certain prominence (as described on peakbagger.com, they appear to manually inspect topos to get their peaks).

A solution for finding the single pixel peaks is to use r.terraflow, followed by r.mapcalc looking for the minimum: Flow Accumulation == 1. As terraflow uses multiple flow directions this method eliminates ridge lines and gives you do with a single flow direction model like D8. This is slow (2 hours for a 15,000 x 15,000 cell DEM), and I am sure there is a more elegant way to do it, but it works.

• I believe you can do better, as you hint at the end. Since you can find peaks, you likely can also find the "cols." Construct all the contours at all the col elevations, assign its elevation to each contour, and perform a Euclidean allocation to fill in the regions enclosed by the contours. The elevation of each peak, minus the value of this grid, tells you the relative height of its prominence: perform your selection accordingly. For your DEM these calculations should take seconds, not hours. – whuber Dec 24 '14 at 21:54

I took a completely different approach. I exported as ASCII using How to produce a CSV file from .tif file with elevation data? and then found the peaks outside QGIS on a 5*5 template using my favourite computer language. In QGIS 2.18.2 this is select layer then Raster->Conversion->Translate and then set the output to ASCII XYZ from the drop-down menu when you specify the output file. I would have preferred Arc/Info ASCII grid (*.asc) format, but the fix for the gotcha described in the above link can't work the same way in 2.18.2.