水文模型流域划分

这周需要利用新生成的流向进行流域划分，所有流入同一条河流的地块单元划分到同一流域中。由于处理回路的问题，有些最后通过沟渠的水最后流向沟渠点而没有流到河流中，但沟渠点的就是属于河流单元，所以要将这一部分计入河流的流域中。开始的思路是先将直接流入河流的地块划分出来，对于剩下的没有划分流域的地块，找到这些地块中与河流相交的地块，之后顺着这个地块递推。判断空间关系时用touches（相邻），结果判断相邻时找到的是河流，不是地块。不清楚为什么，没有解决。第二个思路，利用沟渠点，在dit_pot.shp图层中找到downslope是rive_id的沟渠点，沟渠点的FID加上‘p’对应找到流向中的边，由边开始递推，成功。代码如下`

# -*- coding: utf-8 -*-
"""
Created on Thu Aug  6 21:11:58 2020@author: lzc
"""from shapely.geometry.polygon import Polygon
from osgeo import gdal, osr
from osgeo import ogr
import networkx as nxdef find_upslope(g, down_id, unit_set):for item in g.adj[down_id]:unit_set.add(item)find_upslope(g, item, unit_set)
#找到河流的沟渠点，加上‘p'存入列表中       
def find_dit_basin(dit_filename, river_id):ds = ogr.Open(dit_filename, True)layer = ds.GetLayer()ft = layer.GetNextFeature()dit_list = []while ft is not None:field_id = ft.GetFID()downslope = ft.GetField('downslope')if downslope == river_id:dit_id = 'p' + str(field_id)dit_list.append(dit_id)             ft = layer.GetNextFeature()return dit_list
#根据流向构建网络g
def find_land_upslope(dit_list, flow_filename):g = nx.DiGraph()  ds = ogr.Open(flow_filename, True)layer = ds.GetLayer()ft = layer.GetNextFeature()   while ft is not None:field_id = ft.GetFID()up = ft.GetField('up')down = ft.GetField('down') g.add_node(up)g.add_node(down)g.add_edge(down, up)    ft = layer.GetNextFeature()print( len(g.nodes), len(g.edges) )return g#如果流向中有终点是在沟渠点构成的列表则递推寻找这一路流向，存入集合中
def make(g, dit_list, flow_filename):ds = ogr.Open(flow_filename, True)layer = ds.GetLayer()initialise_set = set()ft = layer.GetNextFeature() while ft is not None:down = ft.GetField('down')if str(down) in dit_list:unit_set = set()find_upslope(g, down, unit_set)initialise_set = initialise_set.union(unit_set)print(initialise_set)ft = layer.GetNextFeature()return initialise_set#对找到的这一路流向在地块单元中对’basin'字段赋值
def modify_lumap(lu_filename, initialise_set, river_id):ds = ogr.Open(lu_filename, True)layer = ds.GetLayer()ft = layer.GetNextFeature()while ft is not None:field_id = ft.GetFID()if str(field_id) in initialise_set:print(field_id)ft.SetField('basin', river_id)layer.SetFeature(ft)#create  setft = layer.GetNextFeature()ds.Destroy()if __name__ == '__main__':river_id = 92749dit_shp = r'C:\Users\lzc\Desktop\1\test\dit_pot.shp'flow_filename = r'C:\Users\lzc\Desktop\1\test\all_flow.shp'lu_filename   = r'C:\Users\lzc\Desktop\1\test\lud0415_Merge.shp'dit_list = find_dit_basin(dit_shp, river_id)print(dit_list)g = find_land_upslope(dit_list, flow_filename)initialise_set = make(g, dit_list, flow_filename)modify_lumap(lu_filename, initialise_set, river_id)

主要用到的是networkx这个包和一个递归。

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！