I asked this question on the osgeo forum too, and some very helpful people replied. Would just like to list out the gist of the answers so that it'd be helpful to anyone who'd search for a contour algorithm later.
DEM or DTM files contain geographic features. Contours aren't geographic features; they're cartographic artifacts. The quality of the geographic features are important.The contour accuracy can never be more accurate than the source data. GDAL isn't the worst or the best. It gives an algorithmic result. Worst or best lies in your data, the processing chain of the data (interpolation/smoothing of DEM etc.) and your intended use.
The only case where route finding where wrong elevation can be life critical would be beach travel beneath sea cliffs and aircraft terrain avoidance.
As for interpolation, there's also the method of using TINS (intersecting the triangular 'faces'). Interpolation doesn't necessarily lead to incorrect results. DEM data is mostly already interpolated.
http://osgeo-org.1560.n6.nabble.com/Most-optimal-algorithm-for-contour-correctness-td5034148.html
I propose a solution with PyQGIS (probably not the most efficient one, but it should work for you).
My idea is focused on the using of the QgsSpatialIndex() and the recurring to spatial queries: because of this, I temporarily convert the input contours to polygons and then do the opposite (the final output is a linestring layer again). Moreover, you didn't answer to my comment about the thinning ratio, so it isn't used in this solution (but it would be easily considered with a simple edit of the code).
Starting from this sample contours (the one on the left is a depression, while the ones in the middle and on the right are hills):
and using this code:
##circles=vector line
##output=output vector
from qgis.core import *
from qgis.PyQt.QtCore import QVariant
layer = processing.getObject(circles)
crs = layer.crs().toWkt()
# Create the output layer
outLayer = QgsVectorLayer('Polygon?crs='+ crs, 'poly' , 'memory')
prov = outLayer.dataProvider()
fields = layer.pendingFields() # Fields from the input layer
fields.append(QgsField('TYPE', QVariant.String, '', 20, 0)) # Name for the new field in the output layer
fields.append(QgsField('GROUP', QVariant.String, '', 20, 0)) # Name for the new field in the output layer
prov.addAttributes(fields) # Add input layer fields to the outLayer
outLayer.updateFields()
# Convert the input contours to polygons
line_to_poly = processing.runalg("qgis:linestopolygons", layer, None)
poly_layer = processing.getObject(line_to_poly['OUTPUT'])
all_features = {} # This dictionary contains all the input features
index = QgsSpatialIndex() # Spatial index
for ft in poly_layer.getFeatures():
index.insertFeature(ft)
all_features[ft.id()] = ft
i=1
j=1
for feat in poly_layer.getFeatures():
inAttr = feat.attributes() # Input attributes
inGeom = feat.geometry() # Input geometry
idsList = index.intersects(inGeom.boundingBox())
if len(idsList) > 0:
temp_group = {}
for id in idsList:
testGeom = all_features[id].geometry()
# Check if the current feature is within or contains the returned id: ...
if inGeom.contains(testGeom) or inGeom.within(testGeom):
# ... if yes, populate the temp_group dictionary for further processing
temp_group[id] = all_features[id]["ELEV"]
if temp_group:
min_val =min(temp_group, key=lambda key: temp_group[key]) # Get the minimum value from temp_group
testGeom1 = all_features[min_val].geometry()
max_val =max(temp_group, key=lambda key: temp_group[key]) # Get the maximum value from temp_group
testGeom2 = all_features[max_val].geometry()
if testGeom1.within(testGeom2):
type = 'Depression'
group = 'D_%s' %i
flag = 1
else:
type = 'Hill'
group = 'H_%s' %j
flag = 0
inAttr.append(type)
inAttr.append(group)
if group.startswith('D'):
i += 1
else:
j += 1
for value in temp_group:
# Delete the features just processed from the spatial index ( it avoids the repeating of the same operation on the other features in the same group)
index.deleteFeature(all_features[value])
outGeom = QgsFeature()
outGeom.setAttributes(inAttr) # Output attributes
outGeom.setGeometry(inGeom) # Output geometry
prov.addFeatures([outGeom]) # Output feature
# Add the layer to the Layers panel
QgsMapLayerRegistry.instance().addMapLayer(outLayer)
# Convert the output to lines
poly_to_line = processing.runalg("qgis:polygonstolines",outLayer,output)
# Add the layer to the Layers panel
QgsMapLayerRegistry.instance().removeMapLayer(outLayer)
I obtain the desired output:
which contains the following attribute table:
Best Answer
It sounds like you have the original surface model so to identify the depressions use a sink algorithm. Here is a link to filling sinks using QGIS.
I have not used this tool in QGIS, only something similar in ArcGIS so I am not sure of the output. In any case, you could subtract the sink output surface from the original surface. Any greater than zero pixel would be a depression.