Stream order calculation in Python
from osgeo import ogr
import matplotlib.pyplot as plt
from pprint import pprint
def find_head_lines(lines):
head_idx = []
num_lines = len(lines)
for i in range(num_lines):
line = lines[i]
first_point = line[0]
has_upstream = False
for j in range(num_lines):
if j == i:
continue
line = lines[j]
last_point = line[len(line)-1]
if first_point == last_point:
has_upstream = True
if not has_upstream:
head_idx.append(i)
return head_idx
def find_next_line(curr_idx, lines):
num_lines = len(lines)
line = lines[curr_idx]
last_point = line[len(line)-1]
next_idx = None
for i in range(num_lines):
if i == curr_idx:
continue
line = lines[i]
first_point = line[0]
if last_point == first_point:
next_idx = i
break
return next_idx
def find_sibling_line(curr_idx, lines):
num_lines = len(lines)
line = lines[curr_idx]
last_point = line[len(line)-1]
sibling_idx = None
for i in range(num_lines):
if i == curr_idx:
continue
line = lines[i]
last_point2 = line[len(line)-1]
if last_point == last_point2:
sibling_idx = i
break
return sibling_idx
def find_prev_lines(curr_idx, lines):
prev_idx = []
# Homework 9: can you do this?
return prev_idx
shapefile_path = "p:/courses/python/streamnetwork.shp"
file = ogr.Open(shapefile_path, 1)
lyr = file.GetLayer(0)
num_features = lyr.GetFeatureCount()
lines = []
for i in range(num_features):
feat = lyr.GetFeature(i)
geom = feat.geometry()
line = geom.GetPoints()
lines.append(line)
head_idx = find_head_lines(lines)
# assign stream order 1 to these lines at head_idx
for idx in head_idx:
print("Head line", idx)
curr_idx = idx
while curr_idx is not None:
print(curr_idx)
curr_idx = find_next_line(curr_idx, lines)
print("Sibling:", find_sibling_line(10, lines))