3

I have line and polygon layers. I want to make a list of points' coordinates (ends or begginings of lines), which are touching the polygons. The coordinates which I need are selected as yellow.

enter image description here

With the code I have right now I get coordinates of the whole line. How to change it? Here is the part of the code and visualisation of my result:

enter image description here

lineFeatures = lines.getFeatures()

lakes=[]

for lineFeat in lineFeatures:
    lineGeom=lineFeat.geometry()
    polygonFeatures=polygons.getFeatures()
    for polygonFeat in polygonFeatures:
        polyGeom=polygonFeat.geometry()
        if lineGeom.touches(polyGeom): 
            lakes.append(lineGeom.asPolyline()[0])
            lakes.append(lineGeom.asPolyline()[1])

The points in the pictures are there just for visualization. In my project, I have only lines and polygons.

I am working with QGIS and python.


From one of the answers, I found a solution. The code below do exactly what I wanted.

lakes=[]

lineFeatures = [feat for feat in lines.getFeatures()]
polyFeatures = [feat for feat in polygons.getFeatures()]

for lineFeat in lineFeatures:
    for polyFeat in polyFeatures:
        if polyFeat.geometry().intersects(lineFeat.geometry()):
            geom = polyFeat.geometry().intersection(lineFeat.geometry())
            line_intersection = geom.asPolyline()
            if line_intersection:
                lakes.append(line_intersection[0])
                lakes.append(line_intersection[1])
  • Your code works for me. Perhaps two aspects: 1) The line 'for polygonFeat in ...' must be indented. 2) Are the geometries actually above each other or a little off? – eftas Jul 7 '16 at 14:26
  • 1) Yes, it is indented in the original version, i just made a mistake with copying. I changed it already in my question. 2) The lines are snapped to polygons, so they are actually above each other. I know also that this code works, but it just not I wanted to obtain and I do not have how to change it :) – sweet.sugar.cola Jul 8 '16 at 8:11
3

Try out with 'intersects' and 'intersection' methods. Next code produces only two points.

mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

feats_line = [ feat for feat in layers[0].getFeatures() ]

feats_poly = [ feat for feat in layers[1].getFeatures() ]

for feat_l in feats_line:
    for feat_p in feats_poly:
        if feat_p.geometry().intersects(feat_l.geometry()):
            geom = feat_p.geometry().intersection(feat_l.geometry())
            print geom.asPolyline()[0]
            print geom.asPolyline()[1]

I tested it with next simple hypothetical situation:

enter image description here

and it works. Points printed at the Python Console were:

(386524,4.45324e+06)
(410342,4.44987e+06)

Editing Note:

With 4 lines:

enter image description here

running this code:

mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

feats_line = [ feat for feat in layers[0].getFeatures() ]

feats_poly = [ feat for feat in layers[1].getFeatures() ]

my_list_points = []

for feat_l in feats_line:
    for feat_p in feats_poly:
        if feat_p.geometry().intersects(feat_l.geometry()):
            geom = feat_p.geometry().intersection(feat_l.geometry())
            my_list_points.append(geom.asPolyline()[0])
            my_list_points.append(geom.asPolyline()[1])

print my_list_points

produces a list with 8 points:

[(387941,4.45518e+06), (410116,4.45158e+06), (385087,4.45128e+06), (411027,4.44466e+06), (383102,4.44575e+06), (408687,4.43823e+06), (392714,4.4348e+06), (405952,4.43338e+06)]
  • This is exactly what I was looking for! Thanks, now it is working! I just have one question: why in my case (4 lines intersecting polygon), the code creates 4 geoms, but only one of them contains points? – sweet.sugar.cola Jul 8 '16 at 11:51
  • Perhaps, the append to the list is not well codified. I slightly modified my code to include it and it works well (see my editing note). – xunilk Jul 9 '16 at 9:02
  • I still get the error: my_list_points.append(geom.asPolyline()[0]) IndexError: list index out of range. I solved it as in the edit up in my question, but I just wanted to know how this function works that it is giving me empty lists... – sweet.sugar.cola Jul 11 '16 at 8:46
1

You need to iterate through the segments of the line to find the one that touches/intersect the polygons (as in Get the vertices on a LineString either side of a Point)

enter image description here

I show here the results with one line and one polygon as example

# transform the line to a list of points
linegeom = liengeom.asPolyline() 
# extract the segments seg_start, seg_end of the line
for seg_start, seg_end in zip(lineGeom, lineGeom[1:]):
   line_start = QgsPoint(seg_start)
   line_end = QgsPoint(seg_end) 
   segment = QgsGeometry.fromPolyline([line_start,line_end])
   if segment.touches(polyGeom): # or if segment.intersects(polyGeom): 
       # segment is the result with line_start and line_end

enter image description here

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.