Counting number of lines connected to point?

I have a water network which and I need to install elbows and Tees (T). When 2 lines connect to a point, I want to install an elbow, when 3 lines connect, I want to install a Tee.

How do I get ArcGIS or QGIS to determine how many lines are connected to a point and designate an elbow or a tee based on that criteria? • What is your data structure like? Separate layers for points and lines? – Simbamangu Feb 16 '14 at 9:08
• Gene is right, make pair then use a dfs algorithm to know how many edges. – user27003 Feb 16 '14 at 14:45

It is not a problem of geometry (buffer, etc.) but a simple problem of Graph Theory (nodes and connected edges) and Python has many modules for dealing with Graphs (PythonGraphApi, Graphs in Python, Python Patterns - Implementing Graphs, ...)

In your problem, you don't need a point layer, you only need the points/nodes of the polyline. They are the vertex/nodes of a Graph and the segments between two nodes are the edges of the same Graph. And you want to find the nodes with a degree 2 and 3 (number of edges incident to the vertex = when 2 lines connect to a point, I want to install an elbow, when 3 lines connect, I want to install a Tee)

The layer: The corresponding Graph where we can see directly the desired points (2 and 3 lines connect - figure made with the NetworkX module): I will use here the NetworkX module (pure Python module). The module can directly open shapefiles to generates a networkx.Graph (Networkx: GIS shapefile) but I do so in the Python console of QGIS:

We need to extract all the segments of the layer (from point to point): the LineString must be iterated over pairs to divide the polyline in segments.

def pairs(list):
'''generator to iterate over pairs in a list '''
for i in range(1, len(list)):
yield list[i-1], list[i]

Creation of the Graph:

import networkx as nx
# select the polyline layer
layer = qgis.utils.iface.activeLayer()
# create a Graph
G = nx.Graph()
# add the nodes and the edges to the Graph = segment point i -> point i+1
for features in layer.getFeatures():
segment = features.geometry().asPolyline()
# add the two end points of a segment as two nodes and an edge (point1 -> point2) to the Graph
for seg_start, seg_end in pair(segment):

Now we have the nodes and the edges of the Graph:

print G.nodes() # = nodes of the polyline in (x,y) form
[(198133,93258.8), (198757,93232.4), (197592,94296.2), (197579,93232.4), (198761,94274.3), (198146,94287.5)]
print G.edges() # = all the segments of the polyline
[((198133,93258.8), (198146,94287.5)), ((198757,93232.4), (198761,94274.3)), ((197592,94296.2), (197579,93232.4)), ((197592,94296.2), (198146,94287.5)), ((198761,94274.3), (198146,94287.5))]

We search the nodes which the degree is > 1

for i in G.nodes():
if G.degree(i) > 1:
print QgsPoint(i), G.degree(i)
(197592,94296.2) 2
(198761,94274.3) 2
(198146,94287.5) 3

And we obtain the desired points:

Result: It depends on Your data structure and precision but I assume that You have separate layers for points and lines. If so, this is relatively easy with Python (either PyQGIS or arcpy) - example with pyqgis could look like this:

#point_layer - vector layer with points
#line_layer - vector layer with lines

iter_p=point_layer.getFeatures() #features from layer
iter_l=list(line_layer.getFeatures()) #changing iterator to list to avoid exhaustion

connections=0 #counter for number of intersections

for p in iter_p:
geom_p=p.geometry().buffer(10,3) #create a small buffer around points, 10 is an arbitrary value, change it according to precision of Your data
for l in iter_l:
geom_l=l.geometry()
if geom_p.intersects(geom_l):
connections=connections+1 # if intersection is true increase counter by one
if connections==2:
#here a set of instructions for doing something - for example changing attrubute value for Tee, elbow or something else
#etc...for each significant value of intersections
connections=0 #reseting counter
• Just checking the buffer isn't the best way. The lines might not connect to the points. – Nathan W Feb 15 '14 at 23:10
• That's true but I've assumed that this is a more real-life problem and data structure is unknown - then it is safer to assume errors etc. - hence very small buffer. If the data is a topologically correct network then Gene answer is better – mrz Feb 17 '14 at 10:43

In Arcgis, the most simple way is to use "spatial join" : this will give you the count of lines that interesect your point. In order to make sure that your points are exactly located at the intersections, an integrate could be necessary.

On the other hand, you should be aware of your topology in order to know what you are counting. You assume here that all your lines are split at point location, but you should make sure of that. Otherwise, you could have points that intersect 2 lines on a Tee. A more robust method would thus consist in creating buffers around each point and counting the number of point intersections between your lines and your buffer, not just the number of lines that intersect.

If you are using ArcMap and have VBA installed then you could adapt the following code. Equally you could translate most of this code in VB .Net

This example selects dangling polylines but if you explore the code you will see that the key test is the line If pTopologyNode.Degree = 1 Then which is testing the valency of the node, so adapting this loop will give you the count you desire.

Public Sub SelectDanglingPolylines()
' Description: Takes a polyline dataset and select all dangling polylines.
'
' Requirements: You need to be in edit mode and have added the layer to a MAP TOPOLOGY,
' also polyline layer must be first in TOC.
'
' Limitations: Large datasets take a long time to build the cache and may even fail.
'
' Author: Duncan Hornby
' Created: 11/12/2011
'

' Get map and then first layer, must be of polyline type
Dim pMXDocument As IMxDocument
Set pMXDocument = ThisDocument
Dim pMap As IMap
Set pMap = pMXDocument.FocusMap
Dim pLayer As ILayer
Set pLayer = pMap.Layer(0)
Dim pFeatureLayer As IFeatureLayer
Set pFeatureLayer = pLayer
Dim pFeatureClass As IFeatureClass
Set pFeatureClass = pFeatureLayer.FeatureClass
If pFeatureClass.ShapeType <> esriGeometryPolyline Then
MsgBox "This code works only with polylines!", vbExclamation, "Wrong data type at layer 0"
Exit Sub
End If

' Get editor and topology extension
Dim pEditor As IEditor
Dim pID As New UID
Dim pTopologyExtension As ITopologyExtension
Dim pTEID As New UID
pID = "esriEditor.editor"
Set pEditor = Application.FindExtensionByCLSID(pID)
pTEID = "esriEditorExt.TopologyExtension"
Set pTopologyExtension = Application.FindExtensionByCLSID(pTEID)
If pTopologyExtension.CurrentTopology Is Nothing Then Exit Sub

' Get a MAP topology not a geodatabase topology
Dim pMapTopology As IMapTopology
If TypeOf pTopologyExtension.CurrentTopology Is IMapTopology Then
Set pMapTopology = pTopologyExtension.MapTopology
Else
' Not a Map Topology
Exit Sub
End If

' This is the colection that FID are added to
Dim aColl As New Collection

' Build cache
Application.StatusBar.Message(0) = "Building MAP TOPOLOGY cache, this can take a long time on large datasets..."
DoEvents
Dim pGeoDataset As IGeoDataset
Set pGeoDataset = pFeatureClass
Dim pEnvelope As IEnvelope
Set pEnvelope = pGeoDataset.Extent
pMapTopology.Cache.Build pEnvelope, False

' Identify dangling nodes and add polyline FID to collection
Application.StatusBar.Message(0) = "Identifying dangling nodes..."
DoEvents
Dim pEnumTopologyParent As IEnumTopologyParent
Dim pTopologyNode As ITopologyNode
Dim pEnumTopologyNode As IEnumTopologyNode
Set pEnumTopologyNode = pMapTopology.Cache.Nodes
pEnumTopologyNode.Reset
Set pTopologyNode = pEnumTopologyNode.Next
While Not pTopologyNode Is Nothing
If pTopologyNode.Degree = 1 Then
' As this has 1 degree it has only 1 parent polyline
Set pEnumTopologyParent = pTopologyNode.Parents
pEnumTopologyParent.Reset
aColl.Add (pEnumTopologyParent.Next.m_FID) 'Adds polyline FID to collection
End If
Set pTopologyNode = pEnumTopologyNode.Next
Wend

' Use collection to select polylines
Application.StatusBar.Message(0) = "Selecting polylines..."
DoEvents
Dim pFeatureSelection As IFeatureSelection
Set pFeatureSelection = pFeatureLayer
Dim X As Variant
For Each X In aColl