6

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?

enter image description here

  • 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
6

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:

enter image description here

The corresponding Graph where we can see directly the desired points (2 and 3 lines connect - figure made with the NetworkX module):

enter image description here

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):
        G.add_edges_from([(seg_start, seg_end)])

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:

enter image description here

4

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
  • 1
    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
2

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.

0

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
 pFeatureSelection.SelectionSet.Add CLng(X)
 Next
 pMXDocument.ActiveView.PartialRefresh esriViewGeoSelection, Nothing, pEnvelope
 Application.StatusBar.Message(0) = ""
 DoEvents
End Sub

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.