# HG changeset patch # User luoboiqingcai # Date 1350984881 -28800 # Node ID 958aa69f1b9ae669eeeed1e1785aae99b26f0191 # Parent ac8ef80687724c49ec714e9d34f8734bc49dab75 alpha0.1 diff --git a/demfunctions.py b/demfunctions.py --- a/demfunctions.py +++ b/demfunctions.py @@ -1,39 +1,63 @@ #!/usr/bin/env python #-*- encoding:utf-8 -*- +import math + def fill_sinks(data): r_max = len(data) c_max = len(data[0]) for r in xrange(len(data)): for c in xrange(len(data[r])): + if r == 0 or r == r_max-1 or c == 0 or c == c_max-1: + continue s0,s1,s2,s3,s4,s5,s6,s7 = [True for i in range(8)] # all the border is highter than this node if r-1>0 and c+1<=c_max-1: s0 = data[r][c] <= data[r-1][c+1] + print 'x0' + print data[r-1][c+1] if c+1<=c_max-1: s1 = data[r][c] <= data[r][c+1] + print 'x1' + print data[r][c+1] if r+1<=r_max-1 and c+1<=c_max-1: s2 = data[r][c] <= data[r+1][c+1] + print 'x2' + print data[r+1][c+1] if r+1<=r_max-1: s3 = data[r][c] <= data[r+1][c] + print 'x3' + print data[r+1][c] if r+1<=r_max-1 and c-1>0: s4 = data[r][c] <= data[r+1][c-1] + print 'x4' + print data[r+1][c-1] if c-1>0: s5 = data[r][c] <= data[r][c-1] + print 'x5' + print data[r][c-1] if r-1>0 and c-1>0: s6 = data[r][c] <= data[r-1][c-1] + print 'x6' + print data[r-1][c-1] if r-1>0: s7 = data[r][c] <= data[r-1][c] + print 'x7' + print data[r-1][c] if False not in [s0,s1,s2,s3,s4,s5,s6,s7]: # sink # window size default to 2*2 print '===================================' print 'center point:r:%s,c:%s'%(r,c) print '===================================' + if r == 8 and c ==8: + print [s0,s1,s2,s3,s4,s5,s6,s7] len_ = 1 marks = [[r,c]] hidden_outset_min = None + border = [] while True: border = [] print '___' + print 'center point:r:%s,c:%s'%(r,c) print 'len_:%s'%len_ print '___' for i in range(len_+1): @@ -71,8 +95,8 @@ border.append(tem) print 'border:%s'%border if len(border) == 0: - print 'the end of dem for this center node.\n this node is assumed to be a sink' - return + print 'center node: %s.\n the end of dem for this center node.\n this node is assumed to be a biggggggg sink'%[r,c] + break hidden_outset = [] for node in border: print 'which node:%s'%node @@ -143,8 +167,143 @@ break if next_border_e == 0: break - print "sink or plane has to extend again" - if hidden_outset_min > data[r][c]: - print "sink" + if len(border) == 8*len_: + if hidden_outset_min > data[r][c]: + print "sink" + for node in marks: + if data[node[0]][node[1]] < hidden_outset_min: + data[node[0]][node[1]] = hidden_outset_min + else: + print "plane" + elif len(border) > 0: + print 'sink on the border of the dem file' else: - print "plane" + print 'too large sink' + return data + +def lift(data,data_modified,p): + r_max = len(data) + c_max = len(data[0]) + for r in xrange(len(data)): + for c in xrange(len(data[r])): + s0,s1,s2,s3,s4,s5,s6,s7 = [True for i in range(8)] # all the border is highter than this node + if r-1>0 and c+1<=c_max-1: + s0 = data[r][c] <= data[r-1][c+1] + if c+1<=c_max-1: + s1 = data[r][c] <= data[r][c+1] + if r+1<=r_max-1 and c+1<=c_max-1: + s2 = data[r][c] <= data[r+1][c+1] + if r+1<=r_max-1: + s3 = data[r][c] <= data[r+1][c] + if r+1<=r_max-1 and c-1>0: + s4 = data[r][c] <= data[r+1][c-1] + if c-1>0: + s5 = data[r][c] <= data[r][c-1] + if r-1>0 and c-1>0: + s6 = data[r][c] <= data[r-1][c-1] + if r-1>0: + s7 = data[r][c] <= data[r-1][c] + if False not in [s0,s1,s2,s3,s4,s5,s6,s7]: # sink + data_modified[r][c] += p + print 'life return' + return data_modified + +def get_vect_martrix(data,x,y): + ''' + x: x distance of each node on column + y: y distance of each node on row + ''' + r_max = len(data) # exactly len(data) - 1 + c_max = len(data[0]) # exactly len(data[0]) - 1 + vect_martix = [[0 for c in range(c_max)] for r in range(r_max)] + for r in xrange(len(data)): + for c in xrange(len(data[r])): + s0,s1,s2,s3,s4,s5,s6,s7 = [True for i in range(8)] # all the border is highter than this node + v0,v1,v2,v3,v4,v5,v6,v7 = [0 for i in range(8)] + if r-1>0 and c+1<=c_max-1: + s0 = data[r][c] <= data[r-1][c+1] + v0 = (data[r][c]-data[r-1][c+1])/math.pow(x**2+y**2,0.5) + if c+1<=c_max-1: + s1 = data[r][c] <= data[r][c+1] + v1 = (data[r][c]-data[r][c+1])/x + if r+1<=r_max-1 and c+1<=c_max-1: + s2 = data[r][c] <= data[r+1][c+1] + v2 = (data[r][c]-data[r+1][c+1])/math.pow(x**2+y**2,0.5) + if r+1<=r_max-1: + s3 = data[r][c] <= data[r+1][c] + v3 = (data[r][c]-data[r+1][c])/y + if r+1<=r_max-1 and c-1>0: + s4 = data[r][c] <= data[r+1][c-1] + v4 = (data[r][c]-data[r+1][c-1])/math.pow(x**2+y**2,0.5) + if c-1>0: + s5 = data[r][c] <= data[r][c-1] + v5 = (data[r][c]-data[r][c-1])/x + if r-1>0 and c-1>0: + s6 = data[r][c] <= data[r-1][c-1] + v6 = (data[r][c]-data[r-1][c-1])/math.pow(x**2+y**2,0.5) + if r-1>0: + s7 = data[r][c] <= data[r-1][c] + v7 = (data[r][c]-data[r-1][c])/y + if False not in [s0,s1,s2,s3,s4,s5,s6,s7]: # sink + pass + which = 0 + v_max = max([v0,v1,v2,v3,v4,v5,v6,v7]) + for i in [v0,v1,v2,v3,v4,v5,v6,v7]: + if v_max == i: + break + which += 1 + if which == None: + raise 'error' + vect_martix[r][c] = which + return vect_martix + +def count_water(data): + r_max = len(data) # exactly len(data) - 1 + c_max = len(data[0]) # exactly len(data[0]) - 1 + water_matrix = [[0 for c in range(c_max)] for r in range(r_max)] + def trackpath(r,c): + #if r > r_max-1 or r < 0 or c>c_max-1 or c < 0: + # return [] + if r>5 or r<0 or c>5 or c<0: + return[] + node = data[r][c] + path = [[r,c],] + if node == 0: + temp = trackpath(r-1,c+1) + if len(temp) >0: + path.extend(temp) + if node == 1: + temp = trackpath(r,c+1) + if len(temp) >0: + path.extend(temp) + if node == 2: + temp = trackpath(r+1,c+1) + if len(temp) >0: + path.extend(temp) + if node == 3: + temp = trackpath(r+1,c) + if len(temp) >0: + path.extend(temp) + if node == 4: + temp = trackpath(r+1,c-1) + if len(temp) >0: + path.exten(temp) + if node == 5: + temp = trackpath(r,c-1) + if len(temp) >0: + path.extend(temp) + if node == 6: + temp = trackpath(r-1,c-1) + if len(temp) >0: + path.extend() + if node == 7: + temp = trackpath(r-1,c) + if len(temp) >0: + path.extend(temp) + return path + for r in xrange(len(data)): + for c in xrange(len(data)): + path = trackpath(r,c) + for node in path: + water_matrix[node[0]][node[1]] += 1 + return water_matrix diff --git a/test.py b/test.py --- a/test.py +++ b/test.py @@ -1,8 +1,72 @@ -from demfunctions import fill_sinks +from demfunctions import (fill_sinks,lift,get_vect_martrix,count_water) +#data_const = [[1,1,3,4,5,5,7,5,4], +# [1,2,4,4,4,4,6,5,3], +# [4,4,3,4,3,3,6,7,5], +# [3,3,2,3,2,2,4,6,7], +# [1,2,2,2,1,2,3,4,5], +# [3,3,2,1,1,2,3,6,6], +# [2,3,3,1,1,2,5,7,8], +# [1,2,3,3,3,3,4,8,6], +# [1,1,2,3,3,4,5,7,8] +# ] def test_fill_sinks(): - data = [[2,2,2,3,3],[3,2,1,1,3],[2,1,1,1,3],[2,2,2,2,3]] + data = [[1,1,3,4,5,5,7,5,4], + [1,2,4,4,4,4,6,5,3], + [4,4,3,4,3,3,6,7,5], + [3,3,2,3,2,2,4,6,7], + [1,2,2,2,1,2,3,4,5], + [3,3,2,1,1,2,3,6,6], + [2,3,3,1,1,2,5,7,8], + [1,2,3,3,3,3,4,8,6], + [1,1,2,3,3,4,5,7,8] ] fill_sinks(data) -test_fill_sinks() +print test_fill_sinks() +#def test_lift(): +# data = [[1,1,3,4,5,5,7,5,4], +# [1,2,4,4,4,4,6,5,3], +# [4,4,3,4,3,3,6,7,5], +# [3,3,2,3,2,2,4,6,7], +# [1,2,2,2,1,2,3,4,5], +# [3,3,2,1,1,2,3,6,6], +# [2,3,3,1,1,2,5,7,8], +# [1,2,3,3,3,3,4,8,6], +# [1,1,2,3,3,4,5,7,8]] +# data_modified = fill_sinks(data) +# return lift(data,data_modified,1e-5) +# +#print test_lift() + +#def test_get_vect_martrix(): +# data = [[1,1,3,4,5,5,7,5,4], +# [1,2,4,4,4,4,6,5,3], +# [4,4,3,4,3,3,6,7,5], +# [3,3,2,3,2,2,4,6,7], +# [1,2,2,2,1,2,3,4,5], +# [3,3,2,1,1,2,3,6,6], +# [2,3,3,1,1,2,5,7,8], +# [1,2,3,3,3,3,4,8,6], +# [1,1,2,3,3,4,5,7,8]] +# data_modified = fill_sinks(data) +# data_lifted = lift(data,data_modified,1e-5) +# return get_vect_martrix(data_lifted,1,1) +# +#print test_get_vect_martrix() + +#def test_count_water(): +# data = [[1,1,3,4,5,5,7,5,4], +# [1,2,4,4,4,4,6,5,3], +# [4,4,3,4,3,3,6,7,5], +# [3,3,2,3,2,2,4,6,7], +# [1,2,2,2,1,2,3,4,5], +# [3,3,2,1,1,2,3,6,6], +# [2,3,3,1,1,2,5,7,8], +# [1,2,3,3,3,3,4,8,6], +# [1,1,2,3,3,4,5,7,8]] +# data_modified = fill_sinks(data) +# data_lifted = lift(data,data_modified,1e-5) +# vect_martrix = get_vect_martrix(data_lifted,1,1) +# return count_water(vect_martrix) +#print test_count_water()