Skip to main content
Source Link
Aquamarine
  • 2k
  • 9
  • 14

This python script which you can run from the console with create points where 2 lines intersect and one or both does not have a vertex at the intersection points. It does not look for self-intersection.

# create point layer with points where either intersecting line or both do not have a vertex

from qgis.core import (
    QgsProject, QgsGeometry, QgsFeature, QgsPointXY, QgsSpatialIndex
)

LAYER_NAME = "pipes"
pipes = QgsProject.instance().mapLayersByName(LAYER_NAME)[0]

epsg = pipes.crs().authid()
out = QgsVectorLayer(f"Point?crs={epsg}", f"{LAYER_NAME}_ints_one_or_none_vertices", "memory")
prov = out.dataProvider()
out.updateFields()

# Build vertex dictionary: { feature_id : set((x,y),...) }
vertex_dict = {}

for f in pipes.getFeatures():
    verts = set()
    for v in f.geometry().vertices():
        verts.add((v.x(), v.y()))
    vertex_dict[f.id()] = verts


index = QgsSpatialIndex(pipes.getFeatures())

added = 0

for f in pipes.getFeatures():
    geomA = f.geometry()
    bbox = geomA.boundingBox()
    candidates = index.intersects(bbox)

    for cid in candidates:

        # skip self-intersections
        if cid == f.id():
            continue

        # ensure pair only processed once
        if cid < f.id():
            continue

        f2 = pipes.getFeature(cid)
        geomB = f2.geometry()

        inter = geomA.intersection(geomB)
        if inter.isEmpty():
            continue

        pts = inter.asMultiPoint() if inter.isMultipart() else [inter.asPoint()]

        for p in pts:
            pt = (p.x(), p.y())

            # Count how many of the two lines have a vertex at this point
            has_vertex_A = pt in vertex_dict[f.id()]
            has_vertex_B = pt in vertex_dict[cid]
            count_vertices = has_vertex_A + has_vertex_B


            # 0 -> no vertex on either line (add)
            # 1 -> exactly one line has a vertex (add)
            # 2 -> both lines have a vertex (skip)
            if count_vertices == 2:
                continue

            # Add the point
            feat = QgsFeature(out.fields())
            feat.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(p)))
            prov.addFeature(feat)
            added += 1

QgsProject.instance().addMapLayer(out)
print(f"Done — added {added} intersection points (where 0 or 1 lines have vertices).")