Finding out and editing specific rows of attribute table in line shapefile using Python?

Question

I got a line shapefile tc_line with 166 rows (line segments), and the figure below was its atttribute table:

enter image description here

Originally, I had to manually select lines (these lines were lines upstream and around confluence point, as showed below) of the shapefile in QGIS and edited values in their rows in the attribute table. The figure below showed those lines (I used invert selection to mark them as little blue segments) and their rows in attribute table:

enter image description here

I want to make the whole process automatic, and recently, with Fiona and networkx, I got edges around confluence point in this shapefile (How to get lines and nodes around the confluence point in a network system (line shapefile)?), and after some works, I obtained the following numpy array storing:

1. coordinates of upstream edges (one row represents one edge) around confluence points (the figure above showed there were 4 confluences, so I got total of 8 coordinates of upstream edges (confluence_x,confluence_y in col.1 and col.2; node_x,node_y in col.3 and col.4))
2. values I want to add in the rows containing those 8 edges in the attribute table (col.6)

Here's the numpy array:

[[ 2.60586260e+05 2.73630282e+06 2.60586260e+05 2.73627393e+06
 2.41920287e+08 5.77421597e-01]
 [ 2.60586260e+05 2.73630282e+06 2.60615144e+05 2.73630282e+06
 1.77046181e+08 4.22578403e-01]
 [ 2.52209841e+05 2.74407267e+06 2.52238726e+05 2.74407267e+06
 5.04274178e+08 9.99998346e-01]
 [ 2.52209841e+05 2.74407267e+06 2.52209841e+05 2.74404378e+06
 8.34297070e+02 1.65444856e-06]
 [ 2.63041417e+05 2.72726206e+06 2.63041417e+05 2.72723318e+06
 1.89988632e+08 9.10935725e-01]
 [ 2.63041417e+05 2.72726206e+06 2.63070301e+05 2.72726206e+06
 1.85756243e+07 8.90642751e-02]
 [ 2.70782383e+05 2.73324109e+06 2.70753498e+05 2.73326997e+06
 5.63667787e+07 9.99955598e-01]
 [ 2.70782383e+05 2.73324109e+06 2.70782383e+05 2.73321221e+06
 2.50289121e+03 4.44016873e-05]]

So now all I need to do is using python to find out rows (line segments), which contain edges listed above, in the attribute table, and then update values in these rows. Obviously that the 8 edges were definitely contained by lines in the shapefile, so I wrote the following script:

import sys, os
import fiona
import networkx as nx
import itertools
import numpy
from shapely.geometry import Point, LineString, mapping, shape
G = nx.Graph()
for line in fiona.open('tc_line.shp'):
 for seg_start, seg_end in itertools.izip(line['geometry']['coordinates'],line['geometry']['coordinates'][1:]):
 G.add_edge(seg_start, seg_end)
for node in G.nodes_iter():
 if G.degree(node) > 2:
 for edge in G.edges(node):
 print edge
abc = numpy.genfromtxt('abc.csv', delimiter = ',')
i = 0
edges = [LineString(edge) for node,edge in itertools.product(G.nodes_iter(), G.edges(node)) if G.degree(node) > 2]
for line in fiona.open('tc_line.shp'):
 for edge in edges:
 while i < abc.shape[0]:
 print i
 print abc[i][:2], abc[i][2:4]
 if shape(line['geometry']).contains(LineString([abc[i][:2], abc[i][2:4]])): 
 print line
 i = i + 1

I read the numpy array from abc.csv, but the script above outputs only one record, and it was strange:

0
[ 260586.25967407 2736302.81560516] [ 260586.25967407 2736273.93140412]
1
[ 260586.25967407 2736302.81560516] [ 260615.14387512 2736302.81560516]
2
[ 252209.84136963 2744072.66568756] [ 252238.72557068 2744072.66568756]
{'geometry': {'type': 'LineString', 'coordinates': [(252238.72557067798, 2744072.6656875624), (252209.84136962818, 2744072.6656875624)]}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'cat_', 73L), (u'value', 64L), (u'label', None)])}
3
[ 252209.84136963 2744072.66568756] [ 252209.84136963 2744043.78148651]
4
[ 263041.41676331 2727262.06067658] [ 263041.41676331 2727233.17647553]
5
[ 263041.41676331 2727262.06067658] [ 263070.30096436 2727262.06067658]
6
[ 270782.38264465 2733241.09029389] [ 270753.4984436 2733269.97449494]
7
[ 270782.38264465 2733241.09029389] [ 270782.38264465 2733212.20609284]

How to use python to find out rows (line segments), which contain coordinates of the edges listed above, in the attribute table, and then update values in these rows?

---

UPDATE#1

I also tried some codes provided by gene:

import sys, os
from shapely.geometry import mapping, shape
import fiona
import networkx as nx
import itertools
from shapely.geometry import Point, LineString
import numpy
# with fiona.open('tc_line.shp', 'r') as input:
# schema = input.schema.copy()
# input_crs = input.crs
# schema['properties']['pi'] = 'float'
# with fiona.open('tc_pi.shp', 'w', 'ESRI Shapefile', schema, input_crs) as output:
# for elem in input:
# elem['properties']['pi']= 1
# output.write({'properties':elem['properties'],'geometry': mapping(shape(elem['geometry']))})
G = nx.Graph()
for line in fiona.open('tc_line.shp'):
 for seg_start, seg_end in itertools.izip(line['geometry']['coordinates'],line['geometry']['coordinates'][1:]):
 G.add_edge(seg_start, seg_end)
for node in G.nodes_iter():
 if G.degree(node) > 2:
 for edge in G.edges(node):
 print 'edge:', edge
edges = [LineString(edge) for node,edge in itertools.product(G.nodes_iter(), G.edges(node)) if G.degree(node) > 2]
lines = [line for line in fiona.open('tc_line.shp') for edge in edges if shape(line['geometry']).contains(edge)]
nodes = [node for node,edge in itertools.product(G.nodes_iter(), G.edges(node)) if G.degree(node) > 2]
print len(edges)
print len(lines)
print len(nodes)

Outputs were very strange and contradictory: On the same basis of G.degree(node) > 2, the first returned 12 and the rest 3 statements returned 8.

edge: ((260586.25967407227, 2736302.815605165), (260586.25967407227, 2736273.931404115))
edge: ((260586.25967407227, 2736302.815605165), (260586.25967407227, 2736331.699806215))
edge: ((260586.25967407227, 2736302.815605165), (260615.1438751221, 2736302.815605165))
edge: ((252209.84136962818, 2744072.6656875624), (252238.72557067798, 2744072.6656875624))
edge: ((252209.84136962818, 2744072.6656875624), (252209.84136962818, 2744101.5498886122))
edge: ((252209.84136962818, 2744072.6656875624), (252209.84136962818, 2744043.7814865126))
edge: ((263041.4167633059, 2727262.060676576), (263041.4167633059, 2727233.1764755263))
edge: ((263041.4167633059, 2727262.060676576), (262983.6483612063, 2727319.8290786757))
edge: ((263041.4167633059, 2727262.060676576), (263070.3009643557, 2727262.060676576))
edge: ((270782.38264465425, 2733241.0902938857), (270753.49844360445, 2733269.9744949355))
edge: ((270782.38264465425, 2733241.0902938857), (270811.26684570406, 2733241.0902938857))
edge: ((270782.38264465425, 2733241.0902938857), (270782.38264465425, 2733212.206092836))
8
8
8

• gene: look at intersecting two shapefiles from Python or command line for example
• Heinz: @gene, I updated my post.
• gene: The correct solution is 12, look in my answer in a modified New . Work directly with LineString and not another format (you don't need NumPy here for example)

Answer 1

In your script:

1. Why do you work with nodes and not directly edges (= LineString and you want the predicate LineString.contains(LineString)) ?
2. Why Numpy here, you can use the simple zip command (edge_set = zip(conflu, edge_node)) ?
3. You forgot a for loop (for lines and for edges

With the example of How to get lines and nodes around the confluence point in a network system (line shapefile)? using itertools.product and list comprehension

import itertools
# with list comprehension
edges = [LineString(edge) for node,edge in itertools.product(G.nodes_iter(), G.edges(node)) if G.degree(node) > 2]
for line in fiona.open('question_stac.shp'):
 for edge in edges:
 if shape(line['geometry']).contains(edge): 
 print line
{'geometry': {'type': 'LineString', 'coordinates': [(89.43117977528091, 7.619382022471937), (139.57162921348316, 1.7205056179775546), (173.49016853932585, -33.18117977528088), (172.50702247191012, -96.5941011235955), (172.50702247191012, -160.49859550561797), (172.50702247191012, -226.86095505617976), (173.49016853932585, -289.7823033707865), (114.00983146067416, -348.27949438202245)]}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'id', None)])}
{'geometry': {'type': 'LineString', 'coordinates': [(89.43117977528091, 7.619382022471937), (139.57162921348316, 1.7205056179775546), (173.49016853932585, -33.18117977528088), (172.50702247191012, -96.5941011235955), (172.50702247191012, -160.49859550561797), (172.50702247191012, -226.86095505617976), (173.49016853932585, -289.7823033707865), (114.00983146067416, -348.27949438202245)]}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'id', None)])}
{'geometry': {'type': 'LineString', 'coordinates': [(172.50702247191012, -160.49859550561797), (236.90308988764048, -161.4817415730337), (330.79353932584274, -255.86376404494382), (357.8300561797753, -227.35252808988764)]}, 'type': 'Feature', 'id': '1', 'properties': OrderedDict([(u'id', None)])}

And the solution of the problem in two lines

edges = [LineString(edge) for node,edge in itertools.product(G.nodes_iter(), G.edges(node)) if G.degree(node) > 2]
lines = [line for line in fiona.open('question_stac.shp') for edge in edges if shape(line['geometry']).contains(edge)]

If you want the nodes, simply use

nodes = [node for node,edge in itertools.product(G.nodes_iter(), G.edges(node)) if G.degree(node) > 2]

New

I convert your results in a shapefile (the edge list) and directly in a Graph

import networkx as nx
G = nx.read_shp('edges.shp')
from shapely.geometry import LineString
for node in G.nodes_iter():
 if G.degree(node) > 2:
 for edge in G.edges(node):
 print LineString(edge)
LINESTRING (252209.8413696282 2744072.665687562, 252238.725570678 2744072.665687562)
LINESTRING (252209.8413696282 2744072.665687562, 252209.8413696282 2744101.549888612)
LINESTRING (252209.8413696282 2744072.665687562, 252209.8413696282 2744043.781486513)
LINESTRING (263041.4167633059 2727262.060676576, 263041.4167633059 2727233.176475526)
LINESTRING (263041.4167633059 2727262.060676576, 262983.6483612063 2727319.829078676)
LINESTRING (263041.4167633059 2727262.060676576, 263070.3009643557 2727262.060676576)
LINESTRING (270782.3826446543 2733241.090293886, 270753.4984436044 2733269.974494935)
LINESTRING (270782.3826446543 2733241.090293886, 270811.2668457041 2733241.090293886)
LINESTRING (270782.3826446543 2733241.090293886, 270782.3826446543 2733212.206092836)
LINESTRING (260586.2596740723 2736302.815605165, 260586.2596740723 2736273.931404115)
LINESTRING (260586.2596740723 2736302.815605165, 260586.2596740723 2736331.699806215)
LINESTRING (260586.2596740723 2736302.815605165, 260615.1438751221 2736302.815605165)

Number of nodes involved

nodes_edges = [G.edges(i) for i in G.nodes_iter() if G.degree(i) > 2]
print len(nodes_edges)
4

Number total of edges

print(len([edge for i,j in enumerate(nodes_edges) for edge in nodes_edges[i]]))
12

Therefore forgot the [LineString(edge) for node,edge in itertools.product(G.nodes_iter(), G.edges(node)) if G.degree(node) > 2] solution for your shapefile.

In conclusion

1) you create a NetworkX Graph with Fiona
2) you use NetworkX to find the edges in a Shapely format (no intermediate formats as Numpy arrays, unnecessary and you don't need to save an intermediate file, only a simple list -> [LineString(edge) for i,j in enumerate(nodes_edges) for edge in nodes_edges[i]])
3) you use Shapely to find the lines/egdes predicate

• Heinz: I don't fully understand your answer, but as I tried edges = [LineString(edge) for node,edge in itertools.product(G.nodes_iter(), G.edges(node)) if G.degree(node) > 2], it printed only 8 edges, not 12, which was the number in my edge_set and also my expectation.
• gene: look above in New
• Heinz: I checked the last set of edges, edge: ((270782.38264465425, 2733241.0902938857), (270753.49844360445, 2733269.9744949355)) edge: ((270782.38264465425, 2733241.0902938857), (270811.26684570406, 2733241.0902938857)) edge: ((270782.38264465425, 2733241.0902938857), (270782.38264465425, 2733212.206092836)), and these edges still intersects.
• gene: How many edges with my script in New ?
• Heinz: I use the same line shapefile in gis.stackexchange.com/a/215880 and ran your script in New, and I got 12 edges.