python中做种子蔓延分析

wet510

近期需要用到种子蔓延分析，做基于DEM的有源洪水淹没分析。简单即是下：有源淹没是指水流收到地表起伏特征的影响，某些地方即使处于低洼处，也可能由于地形阻隔而不会被淹没。也就是说不仅要考虑哪些地方地势低还要考虑这些地方是否联通。而种子蔓延算法指给定空间中一点，从它周围的四个方向或8个方向，考虑：一周围点的地势是否低于种子点；二地势低于种子点的点是否联通；
简单说这些，我找到网上已有的代码：http://cadesign.cn/bbs/thread-742-1-1.html
先在此感谢作者，但由于是草稿，未做优化，在使用中出现了问题。主要问题即是，代码运行效率低，很长时间没有运行结果。现将代码复制如下，请各位大神帮帮忙：

1. import arcpy
   
2. from arcpy import env
   
3. from arcpy.sa import *
   
4. import copy
   
5. 
   
6. #original_para
   
7. worksapce_path=r"D:\GISForDesign_2\data"
   
8. submerge_dem=r"D:\GISForDesign_2\data\DemMerge_Clip.tif"
   
9. submerge_region=r"D:\GISForDesign_2\data\submergesAreaSmall.shp"
   
10. seed=r"D:\GISForDesign_2\data\seed.shp"
    
11. env.overwriteOutput=True
    
12. 
    
13. #environment
    
14. env.workspace=worksapce_path
    
15. env.overwriteOutput=True
    
16. 
    
17. submergeAreaName=worksapce_path+r"\submergeArea.tif"
    
18. submergeArea=arcpy.Clip_management(submerge_dem,"",submergeAreaName,submerge_region,"","ClippingGeometry")
    
19. 
    
20. floodLevel=1200
    
21. 
    
22. #arcpy.SetProgressor("step", "submergeCal...",0, fc_count, 1)
    
23. 
    
24. elevationCell={}
    
25. #rowCount=arcpy.GetRasterProperties_management(submergeArea,"ROWCOUNT")
    
26. #columnCount=arcpy.GetRasterProperties_management(submergeArea,"COLUMNCOUNT")
    
27. cellXSzie=arcpy.GetRasterProperties_management(submergeArea,"CELLSIZEX")
    
28. cellYSzie=arcpy.GetRasterProperties_management(submergeArea,"CELLSIZEY")
    
29. leftPosition=arcpy.GetRasterProperties_management(submergeArea,"LEFT")
    
30. topPosition=arcpy.GetRasterProperties_management(submergeArea,"TOP")
    
31. bottomPosition=arcpy.GetRasterProperties_management(submergeArea,"BOTTOM")
    
32. 
    
33. rstArray=arcpy.RasterToNumPyArray(arcpy.Raster(submergeArea))
    
34. rows,cols=rstArray.shape
    
35. for rowNum in xrange(rows):
    
36.     for colNum in xrange(cols):
    
37.         elevationCell[(rowNum,colNum)]=[rstArray.item(rowNum,colNum),False]
    
38. 
    
39. 
    
40. for row in arcpy.da.SearchCursor(seed,["SHAPE@XY"]):
    
41.     seedXY=row[0]
    
42. del row
    
43. 
    
44. #seedValue=arcpy.GetCellValue_management(submergeArea,"%s %s" % seedXY)
    
45. #originalCordi=arcpy.ExportRasterWorldFile_management (submergeArea)
    
46. 
    
47. 
    
48. 
    
49. XNum=(seedXY[0]-float(leftPosition[0]))//float(cellXSzie[0])
    
50. YNum=(seedXY[1]-float(topPosition[0]))//float(cellYSzie[0])
    
51. seedCoordi=(int(abs(XNum)),int(abs(YNum)))
    
52. 
    
53. 
    
54. def floodCal(floodList,floodLevel):
    
55.     #if len(floodList)>0:         
    
56.     x=floodList[0][0]
    
57.     y=floodList[0][1]
    
58.     elevationCellCopy[floodList[0]][1]=True
    
59.     floodList.pop(0)
    
60.     for i in range(x-1,x+2):
    
61.         for j in range(y-1,y+2):
    
62.             if ((i,j) in  positionKeys and elevationCellCopy[(i,j)][1]==False and elevationCellCopy[(i,j)][0]<=floodLevel):
    
63.                     elevationCellCopy[(i,j)][1]=True
    
64.                     floodList.append((i,j))
    
65.      
    
66. '''                  
    
67. def recursion(floodList,floodLevel):
    
68.     if len(floodList)==0:
    
69.         return 1
    
70.     else:
    
71.         floodCal(floodList,floodLevel)
    
72.         return recursion(floodList,floodLevel)
    
73. '''
    
74. 
    
75. #import sys
    
76. #sys.setrecursionlimit(1000)
    
77. 
    
78. elevationCellCopy=copy.deepcopy(elevationCell)
    
79. positionKeys=elevationCellCopy.keys()
    
80. floodList=[]
    
81. if elevationCell[seedCoordi][0]<=floodLevel:
    
82.     floodList.append(seedCoordi)  
    
83. 
    
84. while True:
    
85.     if (len(floodList)==0):break
    
86.     floodCal(floodList,floodLevel)   
    
87. 
    
88. for rowNum in xrange(rows):
    
89.     for colNum in xrange(cols):
    
90.         rstArray[(rowNum,colNum)]=elevationCellCopy[(rowNum,colNum)][1]
    
91. 
    
92. 
    
93. 
    
94. resultSubmerge=arcpy.NumPyArrayToRaster(rstArray,arcpy.Point(leftPosition[0],bottomPosition[0]),float(cellXSzie[0]),float(cellYSzie[0]))
    
95. 
    
96. 
    
97. 
    
98. #recursion(floodList,floodLevel)                  
    
99.                      
    
100. 
     
101. 
     
102. 
     
103. '''
     
104. for i in range(x-1,x+2):
     
105.     for j in range(y-1,y+2):
     
106.         if ((i,j) in  positionKeys and elevationCellCopy[(i,j)][1]==False and elevationCellCopy[(i,j)][0]<=floodLevel):
     
107.             print((i,j),elevationCellCopy[(i,j)])
     
108. 
     
109. 
     
110. 
     
111. FalseRaster=arcpy.sa.Con(submergeArea,0,0,"Value>=0")
     
112. FalseRaster=submergeArea+FalseRaster
     
113. 
     
114. demToPointsName=worksapce_path+r"\demToPoints.shp"
     
115. demToPoint=arcpy.RasterToPoint_conversion(submergeArea,demToPointsName,"Value")
     
116. '''

复制代码

邺下放歌

不错{:eb315:}{:eb315:}

931404656

{:eb502:}{:eb502:}{:eb502:}{:eb502:}