3

I have a collection of contour LineStrings that all either close where they started or start at an edge and end at another edge. Is is possible to convert these to Polygons for purposes of shading? Each LineString has an associated value (I would like to maintain a min/max value on each Polygon).

This functionality is described here in ArcGIS: https://support.esri.com/en/technical-article/000012180

Example FeatureCollection:

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "type": "LineString",
        "coordinates": [
          [-60, 36],
          [-55, 34],
          [-46, 35],
          [-43, 36 ]
        ]
      }
    },
    {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "type": "LineString",
        "coordinates": [
          [-60, 32],
          [-55, 29],
          [-45, 30],
          [-43, 32 ]
        ]
      }
    },
    {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "type": "LineString",
        "coordinates": [
          [-60, 39],
          [-55, 37],
          [-47, 37],
          [-43, 39 ]
        ]
      }
    },
    {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "type": "LineString",
        "coordinates": [
          [-52, 32],
          [-53, 31],
          [-52, 31],
          [-51, 31],
          [-51, 31],
          [-50, 32],
          [-50, 32],
          [-51, 32],
          [-52, 32]
        ]
      }
    }
  ]
}

enter image description here

I am looking to specify some bounding box and convert these lines to polygons with a result like this:

enter image description here

4
  • If this is a one time conversion this is trivial in any GIS system, for example in ArcMap/ArcPro it is the Feature to Polygon tool.
    – Hornbydd
    Commented Oct 10, 2022 at 15:22
  • Thanks, @Hornbydd. I figured these tools can do this. I'm looking for a compact script to run this on many, many collections.
    – stewart715
    Commented Oct 10, 2022 at 15:40
  • If you have access to ArcPro and the tool I suggest works for you, why not call that tool in a model or a python script, no need to call shapely.
    – Hornbydd
    Commented Oct 10, 2022 at 15:42
  • I can use ArcPro if absolutely necessary but I would prefer a standalone script if possible. Doesn't necessarily have to be shapely.
    – stewart715
    Commented Oct 10, 2022 at 15:59

1 Answer 1

4
+50

For your example to work I added some height values to each line. Additionally I created a sample polygon which intersects with the lines. In QGIS you can than use multiple process-tools in a python script (which you can also use stand-alone):

  1. Split polygon by lines native:splitwSithlines
  2. Join attributes from lines to newly created polygons native:joinattributesbylocation
  3. Create temporary layer with unique geometries native:deleteduplicategeometries
  4. Iterate over the temporary layer and get all height values stored in a dict per feature.
  5. Extract the min and max values
# Access layers
linesLayer = QgsProject.instance().mapLayersByName('sample_lines')[0]
polygonLayer = QgsProject.instance().mapLayersByName('clipping_polygon')[0]

# Split polygon by lines
result1 = processing.run("native:splitwSithlines", {
    'INPUT': polygonLayer,
    'LINES': linesLayer,
    'OUTPUT':'TEMPORARY_OUTPUT'
})

# Join attributes from lines to newly created polygons
result2 = processing.run("native:joinattributesbylocation", {
    'INPUT': result1['OUTPUT'],
    'JOIN': linesLayer,
    'PREDICATE':[0],
    'JOIN_FIELDS':[],
    'METHOD':0,
    'DISCARD_NONMATCHING':False,
    'PREFIX':'',
    'OUTPUT':'TEMPORARY_OUTPUT'
})

# Create temporary layer with unique geometries
result3 = processing.run("native:deleteduplicategeometries", {
    'INPUT': result2['OUTPUT'],
    'OUTPUT':'TEMPORARY_OUTPUT'
})

# Iterate over the temporary layer and get all height values stored in a set per feature
dict = {}
for feature1 in result3['OUTPUT'].getFeatures():
    dict[feature1] = set([float(feature1["height"])])
    for feature2 in result2['OUTPUT'].getFeatures():
        if feature1.geometry().isGeosEqual(feature2.geometry()):
            dict[feature1].add(float(feature2["height"]))

# Create final layer 
resultLayer = QgsVectorLayer("Polygon", "result", "memory")
resultLayerDataProvider = resultLayer.dataProvider()
resultLayerDataProvider.addAttributes([ 
    QgsField("min", QVariant.Double),
    QgsField("max", QVariant.Double)
])
resultLayer.startEditing()

# Iterate over all entries in dict and get min, max values
for feature, values in dict.items():
    fet = QgsFeature()
    fet.setGeometry(feature.geometry())
    
    minimum = min(values)
    maximum = max(values)
    
    fet.setAttributes([minimum, maximum])
    resultLayerDataProvider.addFeature(fet)
    resultLayer.commitChanges()

# Add layer to project
QgsProject.instance().addMapLayer(resultLayer)

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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