I am working on LiDAR data from GIS and R and I want to extract landforms from the DEM. I have found a method to obtain local minima points of my LiDAR data. However, I would like to obtain the geometry of each sink. How can I do that? I have seen that Gradient Vector Flow algorithm (on matlab) and active contour can be used for it. Are there solutions to use it on GIS or R? Or another method?


I'm not sure whether your LiDAR data are raw (i.e. a point cloud) or have already been interpolated to a raster DEM, but if it is the former, you'll need to interpolate the data to a raster DEM first. The next thing you need to do is to perform a 'Depression Filling' operation. Most GIS will have a tool for this task. Once you've created your filled DEM simply minus the original DEM from the filled DEM and you'll have a raster of 'depth in sink', i.e. the depression morphology map.

Here is an example of a DEM of an area in Quebec, Canada, which contains several large depressions that are the result of open-pit mines:

enter image description here

This is the DEM that results from a depression filling operation performed on these same data:

enter image description here

And lastly, here is the result of differencing the original DEM from the filled:

enter image description here

Each grid cell tells you the depth, in this case in metres, within a depression. You'll notice that most grid cells have a zero (white) value and that only those cells associated with a depression have values larger than zero (colours other than white that are associated with the colour gradient).

Many GIS will actually offer a tool that combines the filling and differencing operations in one convenient tool. For example, in Whitebox GAT, the open-source GIS that I develop, there is a tool called 'Depth In Sink' and I believe that there are corresponding tools in other GIS as well (SAGA?). It is a fairly common terrain analysis operation to have to perform.


In the comments below, you mention that you know that some of your depressions are very circular in shape and that when you use the approach described above, the resulting depression boundaries are not very circular. One thing that you might consider doing is thresholding the DEM of difference resulting from the above analysis. For example, you might reclass the raster such that all grid cells with a depth-in-sink value less than the RMSE of your DEM equal to zero. My guess is that the deeper interior of the features are more circular shaped than the lip of the depression, particularly when you account for the DEM error. Naturally this will have the disadvantage of reducing the overall depressional area, but may well improve the depression boundary outlines. In the example that I gave above, here is the result of thresholding the DEM of difference to identify the cells with depth-in-sink values greater than 15 m (the approximate DEM error):

enter image description here

Notice how this better captures the depressions' boundary shape due to the extended lip of the depressions in the original DEM.

  • Thank you for your answer. Actually, I have already done this step of my work. I went even further by extracting local minima points. However, I am not happy with the results of the border extraction. So I would like to use those local minima points in order to draw the contour and to obtain accurate boundaries. – rdmato33 Jun 5 '15 at 15:17
  • @rdmato33 By local minima, do you mean the depression bottom cells? If so, how do you propose using these to draw the depression borders? Also, the approach described above can be used to draw the depression boundaries to the level of accuracy provided by the limit of the DEM data. There is no better way to draw the depression boundaries from the DEM than this approach. A depression's interior is defined by the height of its outlet and by the interior geometry defined by the DEM. – WhiteboxDev Jun 5 '15 at 15:37
  • The contour that you describe would need to be derived from the DEM data itself. This is a further abstraction (resulting from the generalization or smoothing of the contour line) and will necessarily involve a reduction in accuracy compared with the raster polygons defined by the above approach. – WhiteboxDev Jun 5 '15 at 15:40
  • Yes, I mean the depression bottom cells. I have seen in some articles (especially this one : scholarcommons.usf.edu/cgi/…) that it is possible to create a vector flow around the depression bottom and draw a curve around it. Then, they obtain the geometry of the sink. I am not happy with the method described above (perhaps because I want to detect small depressions of 1 meter of diameter and 50 cm of depth, so the smoothing operation does not work well - just for information, my LiDAR is 20 cm of accuracy). – rdmato33 Jun 5 '15 at 15:52
  • 1
    Yes, that worked. Okay, you say that 'the smoothing op does not work well' for you, but the method that I outline does not involve smoothing, whereas the method described in the paper you link to certainly does smooth. I can't figure out why you would prefer a method that 1) interpolates flow direction (which is a model with limited accuracy), then 2) interpolates contours to the flow lines (which is a model with limited accuracy) over a method that maps the geometry directly from the topography of the DEM...especially if your depressions are so small. – WhiteboxDev Jun 5 '15 at 16:00

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.