7

How could I extract all peaks from a DEM (large mountainous area) using GRASS?

Peak means something like: a pixel which is sourrounded by other pixels with smaller elevation numbers in a area arround the "peak-pixel" of about 100m.

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
MartinMap
  • 8,262
  • 7
  • 50
  • 114

6 Answers6

7

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)

Micha
  • 15,555
  • 23
  • 29
  • 3
    Where as this does return a landoform classification it is not under the criteria that was defined in the post. – Jeffrey Evans Dec 08 '13 at 16:55
  • ...but I guess I can work with it! Still testing... – MartinMap Dec 11 '13 at 08: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... – MartinMap Dec 16 '13 at 20:00
  • I have some python code that can help you but am busy rest of today, if you are still still working on this tomorrow I will contribute. Long story short, use moving window operation to search for local maxima. – Kartograaf Aug 25 '22 at 23:53
2

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

markusN
  • 12,976
  • 31
  • 48
1

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.

thayer
  • 141
  • 7
  • 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
1

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.

Community
  • 1
mollwollfumble
  • 51
  • 2
0

I haven't found a prebuilt algorithm, but "prominence" is a defined thing. Basically Everest has a prominence equal to its elevation. So does the highest point on any contiguous land mass. To find the prominence of every peak you step through the contours and find the highest point in each "island" surrounded by the contour "ring" or polygon. The prominence is the elevation minus the lowest contour in which that spot height is still the highest point.

anwjones
  • 1
0

Another approach is as follows (this was my solution when faced with the same problem):

==

In pseudocode:

  • Turn the DEM into contours (polygons) at 1m intervals
  • Make these polygons 'single-part' - i.e. disconnected contour polygons at the same elevation are treated as separate polygons, not a single multipart polygon
  • starting with an empty summits table

Then:

  • FOR elevation = highest_contour_elevation DOWN TO sea_level IN 1m decrements; do the following:
    • For each distinct contour polygon at this elevation; do the following:
      • Look for any existing summits of unknown prominence within the contour polygon
        • if none found (this is a new contour polygon with no existing summits, and no higher contours within it - so must be a new summit), add an new summit at the centre of the contour, summit elevation = contour elevation and mark it as 'unknown prominence'
        • if more than one existing 'unknown prominence' summits are found within the same contour polygon (then we've reached a saddle between 2 summits that we already know about but don't have the prominence for, and this contour polygon contains the defining saddle for the lower summit(s))
          • Mark the prominence of all but the highest contained summit as it's summit elevation minus the current contour elevation, and mark the summit as 'complete'
          • If the summit prominence is less than our minimum (default=30m) then delete the summit as a spurious bump
  • If we've reached sea level, mark any remaining 'unknown prominence' summits as prominence = their elevation, 'complete'.

If, like mine, your DEM does not cover the entire landmass, and additional check is required for each contour polygon:

  • If contour polygon intersects with DEM extent boundary then mark any contained summits of 'unknown prominence' as 'unconfirmed prominence' and set their prominence to their elevation minus the current contour. i.e. these peaks have at least this prominence, but we can't confirm the exact value as we've hit the edge of the DEM.

This took a couple of hours to implement in Ruby + PostGIS. And a couple of days to run for the 87GB Nelson/Tasman 1m x 1m LiDAR dataset (on a budget single-core cloud-based linux server).

==

Cons:

  • slow
  • not scalable
  • runs to only 1m elevation resolution

Pros:

  • 100% reproducable (no fuzzy AI!)
  • Highlights any summits for which it fails to compute prominence
  • Works reliably on ~80GB 1m x 1m LiDAR DEMs for NZ regional council areas
  • Suited my requirements.
madpom
  • 163
  • 6