I have a DEM with small pits (one pit is approximately 5-10 cm wide, here shown as Hillshade).

enter image description here

I want to extract the ridges of these pits to have polygons out of each pit (like this where I did it manually)

enter image description here

The main problem is, that the for one pit the ridge on the right side f.e. is much higher than the one at the left side, because the overall of the wall is quite curved and seen from the front then the ridge on the left could have the same value than the "deepest" point of the pit. I tried to visualize what I mean - that is the wall drawn if you see it from the top.

enter image description here

So, on the one hand I tried to calculate the bigger mean of the whole area and then to nivellate the whole DEM to get it like in the picture above on the right side. However the problem is that my whole area is not only going convex but also concave and in all possible directions so that it is not possible to nivellate it properly.

The solution I reached up to now is, that I calculate the curvature and then thinning the areas. However, the thinning resolved on lines being not properly on the ridges and some weird polygons that do not match at all - see below. I also tried it with Topographic Position Index, however, the ridges are not marked well enough in many cases to close the pits into single polygons.

enter image description here

Does someone have another idea how to be able to extract those ridges as good as possible?

I am open to all kinds of software and solutions: ArcGIS, QGIS, R, Python,...


Not sure this can be fully automated, but workflow below will give you very good estimation of ridges:

arcpy.Curvature_3d(in_raster="DEM_clipped.tif", out_curvature_raster="C:/SCRATCH/CURVA", z_factor="1", out_profile_curve_raster="", out_plan_curve_raster="")
arcpy.gp.RasterCalculator_sa('("CURVA" - "CURVA".minimum)/1000000', "C:/SCRATCH/FLATTENED")
arcpy.gp.Fill_sa("FLATTENED", "C:/SCRATCH/FILLED", "")
arcpy.gp.RasterCalculator_sa('Con("FILLED" > "FLATTENED","FILLED" - "FLATTENED")', "C:/SCRATCH/DEPTH")
arcpy.gp.FocalStatistics_sa("DEPTH", "C:/SCRATCH/FS", "Rectangle 15 15 CELL", "MAXIMUM", "DATA")
arcpy.gp.RasterCalculator_sa('Con("DEPTH" >= "FS",1)', "C:/SCRATCH/TWO")
arcpy.gp.RasterCalculator_sa('Con(IsNull("TWO"),"FLATTENED")', "C:/SCRATCH/HOLES")
arcpy.gp.Fill_sa("HOLES", "C:/SCRATCH/F_HOLES", "")
arcpy.gp.FlowDirection_sa("F_HOLES", "C:/SCRATCH/H_FDIR", "NORMAL", "")
arcpy.gp.RegionGroup_sa("TWO", "C:/SCRATCH/RGP", "FOUR", "WITHIN", "NO_LINK", "")
arcpy.gp.Watershed_sa("H_FDIR", "RGP", "C:/SCRATCH/W4", "VALUE")

In words:

  • use curvature to flatten DEM
  • find 'depth' of water over flatten DEM
  • find local maximums of depth and use them as sinks for watershed analysis

Picture below shows semi-transparent watersheds over original DEM hillshade: enter image description here

  • 2
    Thanks, exactly what I was looking for!
    – Illundra
    Mar 9 '20 at 7:28
  • Interesting surface. Melting ice? Do you mind telling what is it please?
    – FelixIP
    Mar 11 '20 at 4:22
  • Almost - some icebergs do have those forms as well yes, it is generated through waterflow, however, in this case it's limestone. See also: cambridge.org/core/journals/journal-of-fluid-mechanics/article/… for a short description.
    – Illundra
    Mar 11 '20 at 13:12

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.