5

I have contour lines (interval 0.1m), derived from a DEM (1m).

I used the thinning ratio to find only features, which are like circles.

I've imported these features into a Spatialite database.

The image below shows two groups of contour features (there are thousands of them). The left one shows the contour features of a small depression, the right one of a small hill.

enter image description here

My purpose is to find these feature groups, union them (or find the outer ring) and classify them to depressions and hills.

Sample of data:

ID;ELEV;thin_ratio
1;200;0.8

My first idea is to use the ST_Contains function to check if there is a feature with smaller or higher elevation then the feature which overlays and wrap it into a loop.

This query counts the contour features which are within the outer ring of a group. Excluding the outer ring: a.id!=b.id)

SELECT a.id, count(b.id) AS count FROM circle_test a, circle_test b
  WHERE ST_Contains(ST_BuildArea(a.geom),ST_BuildArea(b.geom))
    AND a.id!=b.id
      GROUP BY a.id

The other idea is to intersect them with a Depth in Sink raster (shown as blue in the image), processed with WhiteBox. The raster fills up sinks to a plane surface. The problem is that there are big sinks with both depressions and hills within.

Solutions with Spatialite or Python are favored.

  • Do you want to iterate over all the features or over that ones which have a specified thinning ratio? Is it possibile to have a sample file like the one you posted in the question? – mgri Jan 7 '17 at 13:52
  • And, finally, you only need to separate the sinks from the hills or you want to also group them (i.e.: hill1; hill2; sink1; etc.)? – mgri Jan 7 '17 at 14:07
  • I also want to group them. – Stefan Jan 9 '17 at 12:51
3

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):

enter image description here

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:

enter image description here

which contains the following attribute table:

enter image description here

  • Thank you for your answer. I can test your code in 3 weeks. – Stefan Jan 30 '17 at 19:19
  • I have tried this and I get this error: "Algorithm [Unnamed algorithm] starting... min() arg is an empty sequence See log for more details". Any ideas? Worth mentioning is that my data has 20 m elevation between contour lines but that shouldn't change anything? – oskarlin Feb 27 '17 at 13:26
  • 1
    @OskarKarlin thanks for testing it! I'll give you a feedback as soon as possible! – mgri Feb 27 '17 at 13:35
  • I should maybe say that I run this as a processing script. – oskarlin Feb 27 '17 at 13:35
  • 1
    Cheers. Regarding the crs I tried reprojecting it to another projection which adds crs and still it didn't work. Regarding another question, I already wrote one before I saw this question: gis.stackexchange.com/questions/230068/… – oskarlin Feb 28 '17 at 7:42

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.