Wednesday 23 May 2018

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


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


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...