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:
This is the DEM that results from a depression filling operation performed on these same data:
And lastly, here is the result of differencing the original DEM from the filled:
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):
Notice how this better captures the depressions' boundary shape due to the extended lip of the depressions in the original DEM.