# dem_waters_extractor / demfunctions.py

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150``` ```#!/usr/bin/env python #-*- encoding:utf-8 -*- 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])): 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 # window size default to 2*2 print '===================================' print 'center point:r:%s,c:%s'%(r,c) print '===================================' len_ = 1 marks = [[r,c]] hidden_outset_min = None while True: border = [] print '___' print 'len_:%s'%len_ print '___' for i in range(len_+1): if r+len_ <= r_max-1 and c+i <=c_max-1: tem = [r+len_,c+i] if tem not in border: border.append(tem) if r+len_ <= r_max-1 and c-i >=0: tem= [r+len_,c-i] if tem not in border: border.append(tem) if r-len_ >= 0 and c+i <=c_max-1: tem = [r-len_,c+i] if tem not in border: border.append(tem) if r-len_ >= 0 and c-i >=0: tem = [r-len_,c-i] if tem not in border: border.append(tem) if c+len_ <= c_max-1 and r+i <=r_max-1: tem = [r+i,c+len_] if tem not in border: border.append(tem) if c+len_ <= c_max-1 and r-i >=0: tem = [r-i,c+len_] if tem not in border: border.append(tem) if c-len_ >= 0 and r+i <=r_max-1: tem = [r+i,c-len_] if tem not in border: border.append(tem) if c-len_ >=0 and r-i >=0: tem = [r-i,c-len_] if tem not in border: 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 hidden_outset = [] for node in border: print 'which node:%s'%node loc0 = [node[0]+1,node[1]+1] loc1 = [node[0],node[1]+1] loc2 = [node[0]-1,node[1]+1] loc3 = [node[0]-1,node[1]] loc4 = [node[0]-1,node[1]-1] loc5 = [node[0],node[1]-1] loc6 = [node[0]+1,node[1]-1] loc7 = [node[0]+1,node[1]] t0,t1,t2,t3,t4,t5,t6,t7 = [None for i in range(8)] print 'marks:%s'%marks if loc0 not in marks: if loc0[0]>=0 and loc0[1]>=0 and loc0[0]<=r_max-1 and loc0[1]<=c_max-1: print 'loc0:%s'%loc0 t0 = data[loc0[0]][loc0[1]] if loc1 not in marks: if loc1[0]>=0 and loc1[1]>=0 and loc1[0]<=r_max-1 and loc1[1]<=c_max-1: print 'loc1:%s'%loc1 t1 = data[loc1[0]][loc1[1]] if loc2 not in marks: if loc2[0]>=0 and loc2[1]>=0 and loc2[0]<=r_max-1 and loc2[1]<=c_max-1: print 'loc2:%s'%loc2 t2 = data[loc2[0]][loc2[1]] if loc3 not in marks: if loc3[0]>=0 and loc3[1]>=0 and loc3[0]<=r_max-1 and loc3[1]<=c_max-1: print 'loc3:%s'%loc3 t3 = data[loc3[0]][loc3[1]] if loc4 not in marks: if loc4[0]>=0 and loc4[1]>=0 and loc4[0]<=r_max-1 and loc4[1]<=c_max-1: print 'loc4:%s'%loc4 t4 = data[loc4[0]][loc4[1]] if loc5 not in marks: if loc5[0]>=0 and loc5[1]>=0 and loc5[0]<=r_max-1 and loc5[1]<=c_max-1: print 'loc5:%s'%loc5 t5 = data[loc5[0]][loc5[1]] if loc6 not in marks: if loc6[0]>=0 and loc6[1]>=0 and loc6[0]<=r_max-1 and loc6[1]<=c_max-1: print 'loc6:%s'%loc6 t6 = data[loc6[0]][loc6[1]] if loc7 not in marks: if loc7[0]>=0 and loc7[1]>=0 and loc7[0]<=r_max-1 and loc7[1]<=c_max-1: print 'loc7:%s'%loc7 t7 = data[loc7[0]][loc7[1]] if True in [data[node[0]][node[1]]>i for i in [t0,t1,t2,t3,t4,t5,t6,t7] if i!=None]: # is a hidden output node hidden_outset.append(node) # mark belonging to sink area marks.append(node) if len(hidden_outset) > 0: #enlarge window print 'current hidden_outset number:%s'%len(hidden_outset) else: len_ += 1 continue hidden_outset_min = hidden_outset[0] for node in (hidden_outset): if hidden_outset_min > data[node[0]][node[1]]: hidden_outset_min = data[node[0]][node[1]] print 'hidden_outset:%s'%hidden_outset print 'hidden_outset_min:%s'%hidden_outset_min next_border_e = 0 # escape while loop for node in border: if data[node[0]][node[1]] < hidden_outset_min: len_ += 1 next_border_e = 1 break if next_border_e == 0: break print "sink or plane has to extend again" if hidden_outset_min > data[r][c]: print "sink" else: print "plane" ```