Return to Question

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:

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:

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

and with Fiona and networkx, I also got some edges in this shapefile (stored in numpy array) (How to get lines and nodes around the confluence point in a network system (line shapefile)?):

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)

I want to:

1. Add a column (type:float) containing 1 in the attribute table of tc_line

2. In the table, find out rows (line segments), which contain edges listed above

3. Update value in the new-added column of selected rows

But now I got stuck in stage 2 described above. I wanted to find out lines (rows in the attribute table) containing edges above by setting those edges as condition statement, but I don't know how to construct the condition statement:

import itertools
from shapely.geometry import Point, LineString
import numpy

G = nx.Graph()
for line in fiona.open('tc_line.shp'):
    ...

I am newbie to use Fiona and networkx, and how to solve this and achieve my objectives? I am glad to receive some tips.

I am working with python 2.7.12, Fiona 1.7.0 and networkx 1.11 under Win10 64bits.

UPDATE1

I tried to work out stage 2, and I took intersecting two shapefiles from Python or command line as reference and wrote the following script:

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

#copy and create new shapefile
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_pi.shp'):
    for seg_start, seg_end in itertools.izip(line['geometry']['coordinates'],line['geometry']['coordinates'][1:]):
        G.add_edge(seg_start, seg_end)

#store edges around confluence point in numpy array
conflu = []
edge_node = []
for node in G.nodes_iter():
    if G.degree(node) > 2:
        for edge in G.edges(node):
            conflu.append(edge[0])
            edge_node.append(edge[1])
conflu_set = numpy.array(conflu)            
edge_node_set = numpy.array(edge_node)
edge_set = numpy.hstack((conflu_set, edge_node_set))
print edge_set

#find rows (lines) contain edges
i = 0
for line in fiona.open('tc_pi.shp'):
    while i < edge_set.shape[0]:
        print i
        if shape(line['geometry']).contains(LineString([edge_set[i][:2], edge_set[i][2:]])):    
            print line
        i = i + 1

The edge_set stored coordinates of edges around confluence point:

[[  260586.25967407  2736302.81560516   260586.25967407  2736273.93140412]
 [  260586.25967407  2736302.81560516   260586.25967407  2736331.69980621]
 [  260586.25967407  2736302.81560516   260615.14387512  2736302.81560516]
 [  252209.84136963  2744072.66568756   252238.72557068  2744072.66568756]
 [  252209.84136963  2744072.66568756   252209.84136963  2744101.54988861]
 [  252209.84136963  2744072.66568756   252209.84136963  2744043.78148651]
 [  263041.41676331  2727262.06067658   263041.41676331  2727233.17647553]
 [  263041.41676331  2727262.06067658   262983.64836121  2727319.82907868]
 [  263041.41676331  2727262.06067658   263070.30096436  2727262.06067658]
 [  270782.38264465  2733241.09029389   270753.4984436   2733269.97449494]
 [  270782.38264465  2733241.09029389   270811.2668457   2733241.09029389]
 [  270782.38264465  2733241.09029389   270782.38264465  2733212.20609284]]

and these edges were found from the same line shapefile, but the script above outputs only one record, and it was strange:

0
1
2
3
{'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), (u'pi', 1.0)])}
4
5
6
7
8
9
10
11

I also tried if shape(line['geometry']).intersects(LineString([edge_set[i][:2], edge_set[i][2:]])):, but it outputted even weird result:

0
1
2
3
{'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), (u'pi', 1.0)])}
4
{'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), (u'pi', 1.0)])}
5
{'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), (u'pi', 1.0)])}
6
7
8
9
10
11

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?