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
| import arcpy import copy from arcpy import env from arcpy.sa import *
arcpy.CheckOutExtension("Spatial") env.workspace = "D:\\NonSourceInundation\\result" env.overwriteOutput = True
dem="D:\\NonSourceInundation\\dem.tif"
cellXSzie = float(arcpy.GetRasterProperties_management(dem, "CELLSIZEX")[0]) cellYSzie = float(arcpy.GetRasterProperties_management(dem, "CELLSIZEY")[0]) maxdem=float(arcpy.GetRasterProperties_management(dem,"MAXIMUM")[0]) mindem=float(arcpy.GetRasterProperties_management(dem,"MINIMUM")[0])
def getValueSum(raster): sum=0 rstArray = arcpy.RasterToNumPyArray(raster) rows, cols = rstArray.shape for rowNum in xrange(rows): for colNum in xrange(cols): v=rstArray.item(rowNum, colNum) sum=sum+v return sum
def getVandRaster(cellxsize,cellysize,h): inRaster = Raster(dem) deepRaster = Con(inRaster<=h, h-inRaster, 0) deepRaster2 = SetNull(deepRaster == 0,deepRaster) Vraster=deepRaster*cellxsize*cellysize Vraster2=Con(IsNull(Vraster),0,Vraster) V=getValueSum(Vraster2) if V>=1921911.658: deepRaster2.save("Level"+str(h)+"m_V"+str(V)+"m3.tif") print "Level"+str(h)+"m_V"+str(V)+"m3.tif" return 1 else: return 0
waterLevel=mindem
step=0.01 while waterLevel<maxdem: if getVandRaster(cellXSzie,cellYSzie,waterLevel)==1: break else: waterLevel=waterLevel+step
|