GIS二分法淹没分析

发表于 更新于 分类于 Waline： 本文字数： 236
基于arcpy实现的的GIS水文分析的二分淹没工具，用于计算每个分区的水高和淹没面积、水深等。

本篇内容：

清理文件的时候突然发现自己之前还帮别人写过这个GIS水文分析的二分淹没工具，分享一下。

原理：二分查找

  1. 确定上下界。
  2. 计算mid。(上界 + 下界) / 2
  3. 判断中间点是否满足淹没条件。如果满足,则把中间点设为新的上界;如果不满足,则把中间点设为新的下界。
  4. 重复步骤2-3,直到上下界足够接近,达到所需的精度。
  5. 此时的上界或下界就是近似的淹没高度。

下面是脚本，可以封装到arctoolbox里面。😄

计算水高

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
import arcpy
import os
import re
from arcpy import env
import math

arcpy.CheckOutExtension("3D")
env.workspace = os.getenv("TEMP")

# 加载分区数据
zones = arcpy.GetParameterAsText(0)

# 原始DEM数据
dem = arcpy.GetParameterAsText(1)

output_dir=arcpy.GetParameterAsText(2)

precision = arcpy.GetParameterAsText(3)
# 填洼，看你填不填
# filled_dem = arcpy.sa.Fill(dem)

# 流向栅格
# flow_dir = arcpy.sa.FlowDirection(dem)

# 获取dem的栅格信息
desc = arcpy.Describe(dem)
extent = desc.Extent
cell_size = desc.meanCellWidth

# 创建与dem范围相同、值全为0的栅格
flooded_mask = arcpy.sa.CreateConstantRaster(0, "FLOAT", cell_size, extent)
# 读取分区水量

with arcpy.da.SearchCursor(zones, ["FID", "水量","SHAPE@"]) as cursor:
    for row in cursor:
        oid = row[0]
        if oid<=1188:
            continue
        shape=row[2]
        volume = row[1]
        # 二分法逼近
        clipped_dem = arcpy.sa.ExtractByMask(dem, shape,analysis_extent=extent)
        arcpy.CalculateStatistics_management(clipped_dem)
        upper = arcpy.GetRasterProperties_management(dem, 'MAXIMUM').getOutput(0)
        try:
            lower = arcpy.GetRasterProperties_management(clipped_dem, 'MINIMUM').getOutput(0)
        except:
            upper = 0
            lower = 0
            print("无法获取分区{}的统计数据，请检查分区是否过小或存在其他问题。".format(oid))

        upper = float(upper)
        lower = float(lower)
        h = (upper + lower) / 2
        volume_pattern = re.compile(r'体积=\s*([0-9.]+)')
        this_volume = 0
        while math.fabs(this_volume - volume) > float(precision):
            h = (upper + lower) / 2
            # 计算淹没区域体积
            flooded_volume = arcpy.ddd.SurfaceVolume(clipped_dem, "", "BELOW", str(h), "1", str(cell_size))

            Messages= flooded_volume.getMessages()
            # 没有匹配到volume的情况,需要处理
            match = volume_pattern.search(flooded_volume.getMessages())
            if match:
                this_volume = float(match.group(1))
            else:
                print("未找到体积值")
            if this_volume > volume:
                upper = h
            else:
                lower = h
            if this_volume == 0:
                print("无法计算该分区表面体积：")
                h=0
                break
        print("分区{}水位:{}".format(oid, h))

淹没的深度

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
import re
import arcpy
import os
import re
from arcpy import env

arcpy.CheckOutExtension("3D")
env.workspace = os.getenv("TEMP")
# 加载分区数据
zones = arcpy.GetParameterAsText(0)

# 原始DEM数据
dem = arcpy.GetParameterAsText(1)

output_dir=arcpy.GetParameterAsText(2)

desc = arcpy.Describe(dem)
extent = desc.Extent
cell_size = desc.meanCellWidth

# 创建与dem范围相同、值全为0的栅格
flooded_mask = arcpy.sa.CreateConstantRaster(0, "FLOAT", cell_size, extent)
# 读取分区水量

with open(output_dir+'wh.txt', encoding='utf-8') as f:
    text = f.read()

pattern = r'分区(\d+)水位:(\d+\.?\d*)'

partitions = {}
for match in re.finditer(pattern, text):
    if "无法" in match.group(0):
        continue
    partition = int(match.group(1))
    level = float(match.group(2))
    partitions[partition] = level
count = 7
with arcpy.da.SearchCursor(zones, ["FID", "水量","SHAPE@"]) as cursor:
    for row in cursor:
        oid = row[0]
        shape=row[2]
        h=partitions[oid]
        clipped_dem = arcpy.sa.ExtractByMask(dem, shape,analysis_extent=extent)

        mask=arcpy.sa.Con(arcpy.sa.Raster(clipped_dem) < h, h ,0)
        # mask_dem = arcpy.sa.ExtractByMask(clipped_dem, mask)
        # 将NoData赋值为0
        mask = arcpy.sa.Con(arcpy.sa.IsNull(mask), 0,mask)
        flooded_mask += mask
        # 分步保存
        if oid%100==0:
            count+=1
            mask_dem =arcpy.sa.Con(arcpy.sa.Raster(flooded_mask) > 0, flooded_mask-dem, 0)
            arcpy.CopyRaster_management(mask_dem, output_dir+str(count)+'.tif')
            flooded_mask = arcpy.sa.CreateConstantRaster(0, "FLOAT", cell_size, extent)
# print(partitions)

结果（水高+淹没区域深度）

最后的结果是一个记录各分区水面高度的txt文件和淹没区域深度的tif文件，可以根据自己的要求微调。

水高

淹没区域深度