# Commits

committed 699278d

almost finished

# demfunctions.py

` #!/usr/bin/env python`
` #-*- encoding:utf-8 -*-`
` `
`+'''`
`+This module is developed to extract waters information from DEM format file.`
`+`
`+Use at your own risk.`
`+`
`+Any suggestin is welcome`
`+'''`
`+__authoer__ = 'luoboiqingcai'`
`+__contact__ = 'sf.cumt@gmail.com'`
`+__licence__ = 'MIT'`
`+`
` import math`
` import copy`
`+import logging`
`+`
`+logger = logging.getLogger()`
`+logger.setLevel(logging.WARNING)`
` `
` def fill_sinks(data):`
`     '''`
`-    填充洼地`
`+    find and fill the sinks`
`     '''`
`     r_max = len(data)`
`     c_max = len(data[0])`
`             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]`
`+                logger.debug('x0')`
`+                logger.debug(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]`
`+                logger.debug('x1')`
`+                logger.debug(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]`
`+                logger.debug('x2')`
`+                logger.debug(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]`
`+                logger.debug('x3')`
`+                logger.debug(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]`
`+                logger.debug('x4')`
`+                logger.debug(data[r+1][c-1])`
`             if c-1>=0:`
`                 s5 = data[r][c] <= data[r][c-1]`
`-                print 'x5'`
`-                print data[r][c-1]`
`+                logger.debug('x5')`
`+                logger.debug(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]`
`+                logger.debug('x6')`
`+                logger.debug(data[r-1][c-1])`
`             if r-1>=0:`
`                 s7 = data[r][c] <= data[r-1][c]`
`-                print 'x7'`
`-                print data[r-1][c]`
`+                logger.debug('x7')`
`+                logger.debug(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 '==================================='`
`+                logger.debug('===================================')`
`+                logger.debug('center point:r:%s,c:%s'%(r,c))`
`+                logger.debug('===================================')`
`                 len_ = 0`
`                 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 '___'`
`+                    logger.debug('___')`
`+                    logger.debug('center point:r:%s,c:%s'%(r,c))`
`+                    logger.debug('len_:%s'%len_)`
`+                    logger.debug('___')`
`                     len_ += 1`
`                     for i in range(len_+1):`
`                         if r+len_ <= r_max-1 and c+i <=c_max-1:`
`                             tem = [r-i,c-len_]`
`                             if tem not in border:`
`                                 border.append(tem)`
`-                    print 'border:%s'%border`
`+                    logger.debug('border:%s'%border)`
`                     if len(border) < 8*len_:`
`-                        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]`
`+                        logger.debug('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`
`+                        logger.debug('which node:%s'%node)`
`                         loc0 = [node[0]+1,node[1]+1]`
`                         loc1 = [node[0],node[1]+1]`
`                         loc2 = [node[0]-1,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`
`+                        logger.debug('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`
`+                                logger.debug('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`
`+                                logger.debug('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`
`+                                logger.debug('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`
`+                                logger.debug('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`
`+                                logger.debug('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`
`+                                logger.debug('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`
`+                                logger.debug('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`
`+                                logger.debug('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`
`                         # mark belonging to sink area`
`                         marks.append(node)`
`                     if len(hidden_outset) > 0:`
`-                        print 'current hidden_outset number:%s'%len(hidden_outset)`
`+                        logger.debug('current hidden_outset number:%s'%len(hidden_outset))`
`                     else:`
`                         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`
`+                    logger.debug('hidden_outset:%s'%hidden_outset)`
`+                    logger.debug('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:`
`                         break`
`                 if len(border) == 8*len_:`
`                     if hidden_outset_min > data[r][c]:`
`-                        print "sink"`
`+                        logger.debug("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"`
`+                        logger.debug("plane")`
`     return data`
` `
` def lift(data,p):`
`     '''`
`-    抬升平地`
`+    lift plane which includes plane as the result of filling the sinks and the original plane`
`     '''`
`     r_max = len(data)`
`     c_max = len(data[0])`
`                 modified = True`
`     if modified == False:`
`         return lifted`
`-    print 'life return'`
`+    logger.debug('life return')`
`     return lift(lifted,p)`
` `
`-def get_vect_martrix(data,x,y):`
`+def get_vect_matrix(data,x,y):`
`     '''`
`-    方向代码阵列`
`+    generate orientation matrix`
`     x: x distance of each node on column`
`     y: y distance of each node on row`
`     '''`
`                 v7 = (data[r][c]-data[r-1][c])/y`
`             which = 0`
`             v_max = max([v0,v1,v2,v3,v4,v5,v6,v7])`
`-            print 'v_max:%s'%v_max`
`+            logger.debug('v_max:%s'%v_max)`
`             for i in [v0,v1,v2,v3,v4,v5,v6,v7]:`
`                 if v_max == i:`
`                     break`
` `
` def trackpath(data,r,c,r_ax,c_ax):`
`     '''`
`-    跟踪水流方向`
`+    tracing the river orientation`
`+    Given a node location this function will return a list reprenting the water flow from this node to to the end of dem file.`
`     data: vect_matrix`
`     '''`
`     if r > r_ax-1 or r < 0 or c>c_ax-1 or c < 0:`
` `
` def count_water(data):`
`     '''`
`-    水流线等级`
`+    return catchment area of each node in the matrix`
`     data: vect_matrix`
`     '''`
`     r_max = len(data) # exactly len(data) - 1`
` `
` def river_map(data,threshold=5):`
`     '''`
`-    河流的0-1表示图`
`+    return the bitmap in which 1 means catchment area of this node is larger than threshold and can be considered along the river`
`     '''`
`     r_max = len(data) # exactly len(data) - 1`
`     c_max = len(data[0]) # exactly len(data[0]) - 1`
` `
` def water_grad(cross_sections,paths,r_max,c_max):`
`     '''`
`-    determine the water grad on the basis of paths generated from river_paths function`
`+    determine the water grade on the basis of paths generated from river_paths function`
`     '''`
`     grads = [[0 for i in range(c_max)] for j in range(r_max)]`
`     for x in range(len(paths)):`

# test.py

`-from demfunctions import (fill_sinks,lift,get_vect_martrix,count_water,river_map,river_paths,water_grad)`
`+from demfunctions import (fill_sinks,lift,get_vect_matrix,count_water,river_map,river_paths,water_grad)`
` import pprint`
` `
` #data_const = [[1,1,3,4,5,5,7,5,4],`
` #            [1,1,2,3,3,4,5,7,8]]`
` #    data_modified = fill_sinks(data)`
` #    data_lifted = lift(data_modified,1e-5)`
`-#    return get_vect_martrix(data_lifted,1,1)`
`+#    return get_vect_matrix(data_lifted,1,1)`
` #`
` #print test_get_vect_matrix()`
` `
` #            [1,1,2,3,3,4,5,7,8]]`
` #    data_modified = fill_sinks(data)`
` #    data_lifted = lift(data_modified,1e-5)`
`-#    vect_martrix = get_vect_martrix(data_lifted,1,1)`
`-#    return count_water(vect_martrix)`
`+#    vect_matrix = get_vect_matrix(data_lifted,1,1)`
`+#    return count_water(vect_matrix)`
` #print test_count_water()`
` #`
` #def test_river_map():`
` #            [1,1,2,3,3,4,5,7,8]]`
` #    data_modified = fill_sinks(data)`
` #    data_lifted = lift(data_modified,1e-5)`
`-#    vect_martrix = get_vect_martrix(data_lifted,1,1)`
`-#    return river_map(count_water(vect_martrix),5)`
`+#    vect_matrix = get_vect_matrix(data_lifted,1,1)`
`+#    return river_map(count_water(vect_matrix),5)`
` #`
` #pprint.pprint(test_river_map())`
` `
` #            [1,1,2,3,3,4,5,7,8]]`
` #    data_modified = fill_sinks(data)`
` #    data_lifted = lift(data_modified,1e-5)`
`-#    vect_martrix = get_vect_martrix(data_lifted,1,1)`
`-#    map_ = river_map(count_water(vect_martrix),5)`
`-#    return river_paths(map_,vect_martrix)`
`+#    vect_matrix = get_vect_matrix(data_lifted,1,1)`
`+#    map_ = river_map(count_water(vect_matrix),5)`
`+#    return river_paths(map_,vect_matrix)`
` #`
` #pprint.pprint(test_river_paths())`
` `
`     c_max = len(data[0])`
`     data_modified = fill_sinks(data)`
`     data_lifted = lift(data_modified,1e-5)`
`-    vect_martrix = get_vect_martrix(data_lifted,1,1)`
`-    map_ = river_map(count_water(vect_martrix),5)`
`-    path_dict = river_paths(map_,vect_martrix)`
`+    vect_matrix = get_vect_matrix(data_lifted,1,1)`
`+    map_ = river_map(count_water(vect_matrix),5)`
`+    path_dict = river_paths(map_,vect_matrix)`
`     #return path_dict`
`     #return water_grad(**path_dict)`
`     return water_grad(path_dict['cross_sections'],path_dict['paths'],r_max,c_max)`