一个用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())