一個用python寫的從數字高程格式文件(DEM)中提取水系的模塊

https://bitbucket.org/luoboiqingcai/dem_waters_extractor



本模塊的功能是從dem文件中提取出水系。



模塊經過了基本測試，對《基於數字高程模型的水系提取算法》(地理學與國土研究2000.11;周貴雲,劉瑜,鄔倫)中提到的9*9維度的矩陣進行計算，結果與該文中所得結果徹底一致。



寫這個模塊主要是用於學習目的。因爲不少實際的問題沒有考慮進來，若是用此模塊處理的dem文件的維度比較大或dem文件所描述的地型爲盆地或描述區域內有大的水坑或平底峽谷，程序可能要跑很是長的時候。



本模塊參考了《基於數字高程模型的水系提取算法》(地理學與國土研究2000.11;周貴雲,劉瑜,鄔倫)提到的方法。



下面是測試矩陣：



[[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]]

如下是測試模塊中各個函數的測試代碼：

from demfunctions import (fill_sinks,lift,get_vect_matrix,count_water,river_map,river_paths,water_grad)
import pprint

def test_fill_sinks():
    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] ]
    return fill_sinks(data)

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_modified,0.1)

print test_lift()

def test_get_vect_matrix():
    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_modified,1e-5)
    return get_vect_matrix(data_lifted,1,1)

print test_get_vect_matrix()

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_modified,1e-5)
    vect_matrix = get_vect_matrix(data_lifted,1,1)
    return count_water(vect_matrix)
print test_count_water()

def test_river_map():
    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_modified,1e-5)
    vect_matrix = get_vect_matrix(data_lifted,1,1)
    return river_map(count_water(vect_matrix),5)

pprint.pprint(test_river_map())

def test_river_paths():
    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_modified,1e-5)
    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())

def test_water_grad():
    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]]
    r_max = len(data)
    c_max = len(data[0])
    data_modified = fill_sinks(data)
    data_lifted = lift(data_modified,1e-5)
    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)

pprint.pprint(test_water_grad())