4

enter image description here

I have a DEM available and would like to identify all the peaks that are likely to be bare rock. I want to be able to isolate these areas and create a 5m elevation buffer around them to represent the non-productive forest. I have used tools like fill but find the output to be to large/general and would like to get more specific..


In summary, I want to find and extract each peak in the DEM to turn it into a polygon out to the closed contour that encompasses the peak. Ideally like the picture attached.

  • Perhaps you could use a focal statistics function with a Con, if the value of a cell is greater than the FocalMean, or even better FocalMax with a kernel that excludes the subject cell. Can you draw a picture of what you have and what you want to get so we can better understand what you're after? As it stands there is a good chance that this post will be closed for being unclear what you're asking, it would be good if you clarified your requirements to avoid that. – Michael Stimson Oct 19 '17 at 21:24
  • @MichaelStimson Hi Micheal, the first step I am trying to complete is extracting the peaks in the DEM. Do you know of a way to do this? I have attached an image of what I am trying to get. – user107580 Oct 19 '17 at 22:07
  • Have a go with Con and Focal Stats, that is assuming you have the spatial analyst extension. If you don't have spatial analyst I would say that you can't do this with ArcGIS and need to find a different software package or obtain this extension from your Esri representative. Your data looks somewhat simplified, are you trying to adapt this from sample data or from real world elevations? – Michael Stimson Oct 19 '17 at 22:11
  • @MichaelStimson What is Con? I do have spatial analyst so I will attempt this. That photo is just taken from an article, I have real world data but am attempting to replicate something similar to the attached photo. Thanks for your suggestions! – user107580 Oct 19 '17 at 22:26
  • Con resources.arcgis.com/en/help/main/10.2/index.html#//… it's short for Conditional Raster (or something like that) It would be best to try this with arcpy or Raster Calculator resources.arcgis.com/en/help/main/10.2/index.html#//… to get the concepts clear. Also Focal Statistics resources.arcgis.com/en/help/main/10.2/index.html#//… and NbrAnnulus resources.arcgis.com/en/help/main/10.2/index.html#//… to get the stats of the outside without the focus pixel. – Michael Stimson Oct 19 '17 at 22:30
2

WORKFLOW:

import arcpy
from arcpy.sa import *
from arcpy import env
env.overwriteOutput = True
height=0.5
DEM=arcpy.Raster("dem")
invert=8848.1-DEM
shallow = Fill(invert,0.95*height)
deep = Fill(shallow)
big=Con(deep>shallow,1)
outRgnGrp = RegionGroup (big, "EIGHT")
outZmin = ZonalStatistics(outRgnGrp, "VALUE", invert,"MINIMUM")
AOI=Con((outZmin+height)>invert,outRgnGrp)

pGons=r"..\SCRATCH.gdb\pgons"
peaks=r"..\peaks.shp"
arcpy.RasterToPolygon_conversion(AOI, pGons)
arcpy.Sort_management(pGons, peaks, sort_field="Shape_Area DESCENDING")
arcpy.DeleteIdentical_management(peaks, fields="gridcode")

OUTPUT:

enter image description here

NOTE: I presented workflow as a script, just to make it look shorter. In reality it is just a sequence of existing tools placed in order. You'll need to replace 3 lines, height assignment and output feature classes to make it work.

  • @FelixlP Thanks so much for posting this script! I am a little confused about the output locations for "pgons" and "peaks" (still new to arcpy). Are these shapes supposed to be made beforehand? I have attempted to put a suitable output location but keep getting an error! – user107580 Oct 20 '17 at 20:44
  • No, not at all. Script will create them. Important thing is location of pgons. It has to be created in file gdb. That way shape area will be created and populated automatically. This info used in next line, sorting. Location of peaks is up to you. Don't forget to set processing extent, snap raster and cell size pointing to Dem in environment settings – FelixIP Oct 20 '17 at 20:59
  • @FelixlP what does the height (0.5 in your example) represent? A series of peaks were left out when I ran the tools, specifically ones that are surrounded by a cluster of other peaks. Do you know of a way to prevent this? – user107580 Oct 20 '17 at 21:44
  • 0.5 is equivalent of yours 5 m. I didn't have less flat DEM handy. – FelixIP Oct 20 '17 at 22:25
  • @FelixlP it's all good. Your script has narrowed down my project area immensely. Thank you for sharing – user107580 Oct 20 '17 at 22:37
0

If you're willing/able to step away from Arc, the GRASS addon module r.geomorphon will extract the summits for you - in fact, it's the major use case example in the documentation: https://grass.osgeo.org/grass72/manuals/addons/r.geomorphon.html

  • Hi there, that looks like exactly what I am trying to do! However, I have never used Grass GIS before. Would it take a while to become familiar with the software? Do you have any suggestions for tutorials? Thanks! – user107580 Oct 20 '17 at 16:51
  • Sure, there's a bunch of good stuff here grass.osgeo.org/documentation/first-time-users . FelixLP's script is very good, though. – obrl_soil Oct 21 '17 at 1:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy