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:

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:

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

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