6

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
  • 2
    look at intersecting two shapefiles from Python or command line for example – gene Oct 30 '16 at 16:01
  • @gene, I updated my post. – Heinz Oct 31 '16 at 4:38
  • 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) – gene Nov 1 '16 at 17:17
4

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

  • 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. – Heinz Oct 31 '16 at 13:27
  • look above in New – gene Oct 31 '16 at 16:49
  • 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. – Heinz Oct 31 '16 at 17:07
  • How many edges with my script in New ? – gene Oct 31 '16 at 17:12
  • I use the same line shapefile in gis.stackexchange.com/a/215880 and ran your script in New, and I got 12 edges. – Heinz Nov 1 '16 at 6:59

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.