GISHydroNXT System Documentation
#*********************************************************************************
# Author: UMD
# Date: 24-07-2018
# Modified: n/a
# Classes: BasinStatistics()
# Functions: hydro.channelslope() ; hydro.slopegrid() ; hydro.SSURGOPct()
# hydro.GetLUCountNLCD() ; hydro.GetLUCountAnderson()
# hydro.GetLUCountMDDE() ; hydro.GetLUCountUltimate()
# hydro.GetLUCountMRLC() ; hydro.GetLUCountUSGS()
# hydro.GetImpCountNLCDFair() ; hydro.GetImpCountNLCDGood()
# hydro.GetImpCountNLCDPoor() ; hydro.GetImpCountAndersonFair()
# hydro.GetImpCountAndersonGood() ; hydro.GetImpCountAndersonPoor()
# hydro.GetImpCountMDDEFair() ; hydro.GetImpCountMDDEGood()
# hydro.GetImpCountMDDEPoor() ; hydro.GetImpCountUltimateFair()
# hydro.GetImpCountUltimateGood() ; hydro.GetImpCountUltimatePoor()
# hydro.GetImpCountMRLCFair() ; hydro.GetImpCountMRLCGood()
# hydro.GetImpCountMRLCPoor() ; hydro.GetImpCountUSGSFair()
# hydro.GetImpCountUSGSGood() ; hydro.GetImpCountUSGSPoor()
# hydro.SoilPct()
# Modules: arcpy ; datetime ; webbrowser
# Comments: n/a
#*********************************************************************************
class BasinStatistics(object):
"""Implementation for GISHydroNXT_addin.button3 (Button)"""
def __init__(self):
self.enabled = False
self.checked = False
def onClick(self):
arcpy.env.scratchWorkspace = scratchfolder
arcpy.env.workspace = optfolder
button1.enabled = False
button3.enabled = False
opath = optfolder
pythonaddins.MessageBox("Basin Statistics calculation will continue. Please wait for analysis output""\n"
"to appear at the interface in the form of a message box.", "Basin Statistics Analysis", 0)
#*******************************************************************************************************
# Warning messages
#*******************************************************************************************************
Impwarntext = """
******************************************************
IMPERVIOUS AREA IN WATERSHED EXCEEDS 10%.
Calculated discharges from USGS Regression
Equations may not be appropriate.
******************************************************
"""
provwarntext = """
******************************************************
Watershed is within 5km of physiographic
province boundary. You should consider
sensitivity of discharges to region location.
******************************************************
"""
limewarntext = """
******************************************************
Watershed is within 1km of underlying limestone
geology. You should consider sensitivity
of discharges to percent limestone calculated.
******************************************************
"""
#*******************************************************************************************************
# Get outlet coordinates and prepare masked grids for calculations
#*******************************************************************************************************
outletcell = optfolder + "/outletcell"
rast = arcpy.Raster(outletcell)
global xoutlet
xoutlet = float(xoutletstring)
global youtlet
youtlet = float(youtletstring)
basingrid = optfolder + "/basingrid"
dirgrid = arcpy.sa.Times(optfolder + "/flowdir_dem",basingrid)
elevgrid = arcpy.sa.Times(optfolder + "/dem",basingrid)
lantype = arcpy.sa.Times(optfolder + "/landuse",basingrid)
# Get basingrid count [number of pixels]
shedtab = arcpy.SearchCursor(basingrid,"","","Count","")
for row in shedtab:
basinarea = row.getValue("Count")
del row
#*******************************************************************************************************
# Compute channel and land slope
#*******************************************************************************************************
theslope = hydro.channelslope(dirgrid,elevgrid,opath) # already converted into feet/mile
theslope_feet = float(theslope/5280.0)
global thechannelslope # 1/18/2018: global variable created for use in discharge comparison
thechannelslope = theslope_feet
maxlength = float(hydro.maxlength) # already converted into miles for use in channel slope function
src = r"" + Directory + "/data/slope_calc/calcslope.exe"
slopegrid = hydro.slopegrid(dirgrid,elevgrid,src,opath)
landsloperesult = arcpy.GetRasterProperties_management(slopegrid,"MEAN")
landslopevalue = float(landsloperesult.getOutput(0))
global landslope
landslope = float(landslopevalue)/3.2808 # modified: 03/19/2014 (divided by 3.2808)
#*******************************************************************************************************
# Determine province of outlet location
#*******************************************************************************************************
outletcell = arcpy.sa.Int(optfolder + "/outletcell")
outletcell.save(optfolder + "/outcell")
arcpy.RasterToPolygon_conversion(optfolder + "/outcell", optfolder + "/outletpoly","NO_SIMPLIFY","COUNT")
prov = r"" + Directory + "/data/maryland/Mdprov.shp"
outletpoly = optfolder + "/outletpoly.shp"
outpoint = optfolder + "/outletpoint.shp"
arcpy.Clip_analysis(prov, outletpoly, outpoint)
ProvTab = arcpy.SearchCursor(outpoint,"","","PROVINCE","")
for row in ProvTab:
if row.getValue("PROVINCE") == "A":
provstring = "Appalachian Plateaus and Allegheny Ridges"
elif row.getValue("PROVINCE") == "B":
provstring = "Blue Ridge and Great Valley"
elif row.getValue("PROVINCE") == "P":
provstring = "Piedmont"
elif row.getValue("PROVINCE") == "W":
provstring = "Western Coastal Plain"
elif row.getValue("PROVINCE") == "E":
provstring = "Eastern Coastal Plain"
else:
row.getValue("PROVINCE") == "Unknown"
provstring == "Unknown"
del row
#*******************************************************************************************************
# Get percent soil types from SSURGO soil dataset
#*******************************************************************************************************
ssurgo = arcpy.sa.Times(optfolder + "/ssurgo",basingrid)
pct = hydro.SSURGOPct(basinarea,ssurgo)
pctSoil = map(float,pct)
pctAsoil = float(pctSoil[0])
pctBsoil = float(pctSoil[1])
pctCsoil = float(pctSoil[2])
pctDsoil = float(pctSoil[3])
pctWsoil = float(pctSoil[4])
#*******************************************************************************************************
# Get LU count for Urban, Nil, Forest, and Storage -- it will be different for
# MOP, MD/DE, MRLC, and USGS
#*******************************************************************************************************
if landuse == "NLCD 2011 Landuse":
count = hydro.GetLUCountNLCD(basinarea,lantype)
LUcount = map(float,count)
UrbPct = float(LUcount[0])
FC = float(LUcount[2])
ST = float(LUcount[3])
elif landuse == "NLCD 2006 Landuse":
count = hydro.GetLUCountNLCD(basinarea,lantype)
LUcount = map(float,count)
UrbPct = float(LUcount[0])
FC = float(LUcount[2])
ST = float(LUcount[3])
elif landuse == "NLCD 2001 Landuse":
count = hydro.GetLUCountNLCD(basinarea,lantype)
LUcount = map(float,count)
UrbPct = float(LUcount[0])
FC = float(LUcount[2])
ST = float(LUcount[3])
elif landuse == "2010 MOP Landuse":
count = hydro.GetLUCountAnderson(basinarea,lantype)
LUcount = map(float,count)
UrbPct = float(LUcount[0])
FC = float(LUcount[2])
ST = float(LUcount[3])
elif landuse == "2002 MOP Landuse":
count = hydro.GetLUCountAnderson(basinarea,lantype)
LUcount = map(float,count)
UrbPct = float(LUcount[0])
FC = float(LUcount[2])
ST = float(LUcount[3])
elif landuse == "1997 USGS Landuse":
count = hydro.GetLUCountAnderson(basinarea,lantype)
LUcount = map(float,count)
UrbPct = float(LUcount[0])
FC = float(LUcount[2])
ST = float(LUcount[3])
elif landuse == "2002 MD/DE Landuse":
count = hydro.GetLUCountMDDE(basinarea,lantype)
LUcount = map(float,count)
UrbPct = float(LUcount[0])
FC = float(LUcount[2])
ST = float(LUcount[3])
elif landuse == "Ultimate Landuse":
count = hydro.GetLUCountUltimate(basinarea,lantype)
LUcount = map(float,count)
UrbPct = float(LUcount[0])
FC = float(LUcount[2])
ST = float(LUcount[3])
elif landuse == "MRLC Landuse":
count = hydro.GetLUCountMRLC(basinarea,lantype)
LUcount = map(float,count)
UrbPct = float(LUcount[0])
FC = float(LUcount[2])
ST = float(LUcount[3])
else: #landuse == "1970s USGS Landuse":
count = hydro.GetLUCountUSGS(basinarea,lantype)
LUcount = map(float,count)
UrbPct = float(LUcount[0])
FC = float(LUcount[2])
ST = float(LUcount[3])
#*******************************************************************************************************
# Get Impervious count -- count will vary depending upon choice of input Landuse and Hyd condition
#*******************************************************************************************************
# For Landuse: "NLCD 2011 Landuse"
if landuse == "NLCD 2011 Landuse":
if hyd == "Fair":
Impcount = hydro.GetImpCountNLCDFair(lantype)
elif hyd == "Good":
Impcount = hydro.GetImpCountNLCDGood(lantype)
elif hyd == "Poor":
Impcount = hydro.GetImpCountNLCDPoor(lantype)
# For Landuse: "NLCD 2006 Landuse"
if landuse == "NLCD 2006 Landuse":
if hyd == "Fair":
Impcount = hydro.GetImpCountNLCDFair(lantype)
elif hyd == "Good":
Impcount = hydro.GetImpCountNLCDGood(lantype)
elif hyd == "Poor":
Impcount = hydro.GetImpCountNLCDPoor(lantype)
# For Landuse: "NLCD 2001 Landuse"
if landuse == "NLCD 2001 Landuse":
if hyd == "Fair":
Impcount = hydro.GetImpCountNLCDFair(lantype)
elif hyd == "Good":
Impcount = hydro.GetImpCountNLCDGood(lantype)
elif hyd == "Poor":
Impcount = hydro.GetImpCountNLCDPoor(lantype)
# For Landuse: "2010 MOP Landuse", "2002 MOP Landuse", and "1997 USGS Landuse"
if landuse == "2010 MOP Landuse":
if hyd == "Fair":
Impcount = hydro.GetImpCountAndersonFair(lantype)
elif hyd == "Good":
Impcount = hydro.GetImpCountAndersonGood(lantype)
elif hyd == "Poor":
Impcount = hydro.GetImpCountAndersonPoor(lantype)
if landuse == "2002 MOP Landuse":
if hyd == "Fair":
Impcount = hydro.GetImpCountAndersonFair(lantype)
elif hyd == "Good":
Impcount = hydro.GetImpCountAndersonGood(lantype)
elif hyd == "Poor":
Impcount = hydro.GetImpCountAndersonPoor(lantype)
if landuse == "1997 MOP Landuse":
if hyd == "Fair":
Impcount = hydro.GetImpCountAndersonFair(lantype)
elif hyd == "Good":
Impcount = hydro.GetImpCountAndersonGood(lantype)
elif hyd == "Poor":
Impcount = hydro.GetImpCountAndersonPoor(lantype)
# For Landuse: "2002 MD/DE Landuse"
if landuse == "2002 MD/DE Landuse":
if hyd == "Fair":
Impcount = hydro.GetImpCountMDDEFair(lantype)
elif hyd == "Good":
Impcount = hydro.GetImpCountMDDEGood(lantype)
elif hyd == "Poor":
Impcount = hydro.GetImpCountMDDEPoor(lantype)
# For Landuse: "Ultiamte Landuse"
if landuse == "Ultimate Landuse":
if hyd == "Fair":
Impcount = hydro.GetImpCountUltimateFair(lantype)
elif hyd == "Good":
Impcount = hydro.GetImpCountUltimateGood(lantype)
elif hyd == "Poor":
Impcount = hydro.GetImpCountUltimatePoor(lantype)
# For Landuse: "MRLC Landuse"
if landuse == "MRLC Landuse":
if hyd == "Fair":
Impcount = hydro.GetImpCountMRLCFair(lantype)
elif hyd == "Good":
Impcount = hydro.GetImpCountMRLCGood(lantype)
elif hyd == "Poor":
Impcount = hydro.GetImpCountMRLCPoor(lantype)
# For Landuse: "1970s USGS Landuse"
if landuse == "1970s USGS Landuse":
if hyd == "Fair":
Impcount = hydro.GetImpCountUSGSFair(lantype)
elif hyd == "Good":
Impcount = hydro.GetImpCountUSGSGood(lantype)
elif hyd == "Poor":
Impcount = hydro.GetImpCountUSGSPoor(lantype)
global IA
IA = float((Impcount/basinarea)* 100)
#*******************************************************************************************************
# Get Limestone percent count
#*******************************************************************************************************
arcpy.env.extent = "MAXOF"
limestonem = r"" + Directory + "/data/maryland/limestonem.shp"
LIcnt = 0
limegrid = arcpy.sa.ExtractByMask("basingrid",limestonem)
limegrid.save(optfolder + "/limegrid")
arcpy.BuildRasterAttributeTable_management(optfolder + "/limegrid","Overwrite")
with arcpy.da.SearchCursor(optfolder + "/limegrid","Count") as rows:
for row in rows:
LIcnt += row[0] or 0
global LI
LI = float((float(LIcnt)/basinarea)* 100) # 10-23-2013
LI = float((float(LIcnt)/basinarea)* 100) # 10-23-2013
#*******************************************************************************************************
# can"t add unit error message -- data is tested and have linear units
#*******************************************************************************************************
cellsize = rast.meanCellWidth
cellsq = cellsize* cellsize
global areami2
areami2 = float((basinarea* cellsq) / 2588881) # conversion into sq miles
#*******************************************************************************************************
# Get basi relief [it is difference of mean elevation and outlet elevation]
#*******************************************************************************************************
elev1 = arcpy.GetRasterProperties_management(elevgrid,"MEAN")
mean_elev = float(elev1.getOutput(0))
global basinrelief
basinrelief = float(mean_elev - outletelev) # Assuming it is already converted into feets
#*******************************************************************************************************
# Average CN number using above outgrid depending upon user choice of HydCon
#*******************************************************************************************************
cnGrid = arcpy.sa.Times(optfolder + "/curveNumber",basingrid)
avgCN_val = arcpy.GetRasterProperties_management(cnGrid,"MEAN") # figure out outgrid source
global avgCN
avgCN = float(avgCN_val.getOutput(0))
#*******************************************************************************************************
# percent soil types except SSURGO
#*******************************************************************************************************
if soil == "STATSGO Soils":
tempsoilgrid = arcpy.sa.Times(basingrid, optfolder + "/soilG_a")
pctAR = arcpy.GetRasterProperties_management(tempsoilgrid, "MEAN")
pctAsoilR = float(pctAR.getOutput(0))
tempsoilgrid = arcpy.sa.Times(basingrid, optfolder + "/soilG_b")
pctBR = arcpy.GetRasterProperties_management(tempsoilgrid, "MEAN")
pctBsoilR = float(pctBR.getOutput(0))
tempsoilgrid = arcpy.sa.Times(basingrid, optfolder + "/soilG_c")
pctCR = arcpy.GetRasterProperties_management(tempsoilgrid, "MEAN")
pctCsoilR = float(pctCR.getOutput(0))
tempsoilgrid = arcpy.sa.Times(basingrid, optfolder + "/soilG_d")
pctDR = arcpy.GetRasterProperties_management(tempsoilgrid, "MEAN")
pctDsoilR = float(pctDR.getOutput(0))
else:
# cell count from SSURGO or RAGAN soil dataset
#*****************************************************************************
ssurgan = arcpy.sa.Times(optfolder + "/soils",basingrid)
pctR = hydro.SoilPct(basinarea,ssurgan)
## pctSoil = map(int,pctR)
pctSoil = map(float,pctR)
pctAR = float(pctSoil[0])
pctAsoilR = float((pctAR/basinarea)* 100)
pctBR = float(pctSoil[1])
pctBsoilR = float((pctBR/basinarea)* 100)
pctCR = float(pctSoil[2])
pctCsoilR = float((pctCR/basinarea)* 100)
pctDR = float(pctSoil[3])
pctDsoilR = float((pctDR/basinarea)* 100)
"""
The following code calculates the Time of Concentration. If multiple provinces
are involved, tc is weighted average of area of watershed in each province.
More correct would be to perform weighted average based on length of channel
in each province. This modification will be performed at a later time. (GEM - 12/01/99)
"""
#*******************************************************************************************************
# Create Zonal stats table based on shedpoly and Mdprov"s Province field
# and add fields
#*******************************************************************************************************
Mdprov = r"" + Directory + "/data/maryland/Mdprov.shp"
theVTab = optfolder + "/theVTab.dbf"
#*******************************************************************************************************
# don"t add "theVTab" to TOC -- it will change list by drawing order to list
# by source which will prohibit addition of new layers
#*******************************************************************************************************
arcpy.sa.ZonalStatisticsAsTable(Mdprov,"PROVINCE","basingrid","theVTab","DATA","ALL")
arcpy.DeleteField_management("theVTab","ZONE_CODE;MIN;MAX;RANGE;MEAN;STD;SUM;VARIETY;MAJORITY;MINORITY;MEDIAN")
addFieldNameList = ["Q1.25", "Q1.50", "Q1.75", "Q2", "Q5", "Q10", "Q25", "Q50", "Q100", "Q200","Q500"]
for each in addFieldNameList:
arcpy.AddField_management("theVTab",each,"FLOAT",10,3)
#*******************************************************************************************************
# Create regionlist and regionarea and declare them as global for use in
# Thomas Discharge script (is it used in Thomas Discharge script?)
#*******************************************************************************************************
sumarea = 0
theVTab = arcpy.SearchCursor("theVTab","","","Count","")
for each in theVTab:
count = each.getValue("Count")
sumarea = sumarea+count
sumArea = sumarea
del each
regionlist = []
regionarea = []
breakstring = ""
theVTab = arcpy.SearchCursor("theVTab","","","Province;Count","")
for row in theVTab:
theProv = row.getValue("Province")
theArea = float(row.getValue("Count"))
areapercent = float((theArea/sumArea)* 100)
areapercent = "{0:.2f}".format(areapercent)
regionlist.append(theProv)
regionarea.append(areapercent)
if row.getValue("Province") == "A":
breakstring = breakstring + " -Appalachian Plateaus and Allegheny Ridges %s percent of area" "\n" %(areapercent)
elif row.getValue("Province") == "B":
breakstring = breakstring + " -Blue Ridge and Great Valley %s percent of area" "\n" %(areapercent)
elif row.getValue("Province") == "P":
breakstring = breakstring + " -Piedmont %s percent of area" "\n" %(areapercent)
elif row.getValue("Province") == "W":
breakstring = breakstring + " -Western Coastal Plain %s percent of area" "\n" %(areapercent)
elif row.getValue("Province") == "E":
breakstring = breakstring + " -Eastern Coastal Plain %s percent of area" "\n" %(areapercent)
del row
#*******************************************************************************************************
# Compute Time of Concentration:
# 1] W.O. Thomas, Jr. Equation [tc]
# 2] SCS Lag equation * 1.67 [lagtime]
#*******************************************************************************************************
sumtc = 0
theVTab = arcpy.SearchCursor("theVTab","","","Province;Count","")
for row in theVTab:
theProv = row.getValue("Province")
theArea = row.getValue("Count")
if row.getValue("Province") == "A":
temptc = 0.133* ((maxlength)**(0.475))* ((theslope)**(-0.187))* ((101 - FC)**(-0.144))* ((101 - IA)**(0.861))* ((ST + 1)**(0.154))* ((10)**(0.194))
elif row.getValue("Province") == "W":
temptc = 0.133* ((maxlength)**(0.475))* ((theslope)**(-0.187))* ((101 - FC)**(-0.144))* ((101 - IA)**(0.861))* ((ST + 1)**(0.154))* ((10)**(0.366))
elif row.getValue("Province") == "E":
temptc = 0.133* ((maxlength)**(0.475))* ((theslope)**(-0.187))* ((101 - FC)**(-0.144))* ((101 - IA)**(0.861))* ((ST + 1)**(0.154))* ((10)**(0.366))
else:
temptc = 0.133* ((maxlength)**(0.475))* ((theslope)**(-0.187))* ((101 - FC)**(-0.144))* ((101 - IA)**(0.861))* ((ST + 1)**(0.154))
sumtc = sumtc + (temptc* theArea)
del row
global tc
tc = (sumtc/basinarea)
#*******************************************************************************************************
# Calculate lagtime
#*******************************************************************************************************
global lagtime
lagtime = ((np.float64(100* ((maxlength* 5280)**(0.8))* (((1000/avgCN) - 9)**(0.7))) / (1900* ((abs(landslope)* 100)**(0.5)))) / 60)
#*******************************************************************************************************
# Calculate Mean Annual Precipitation
#*******************************************************************************************************
arcpy.env.scratchWorkspace = scratchfolder
arcpy.env.workspace = optfolder
mapstpm = r"" + Directory + "/data/maryland/map/mapstpm"
prec_grid = r"" + Directory + "/data/prec/p2-24m"
arcpy.env.cellSize = str(cellsize)
basingrid_p = arcpy.Raster(basingrid)
analysis_extent = basingrid_p.extent
arcpy.env.extent = analysis_extent
maprecbasin = arcpy.sa.Times(mapstpm,basingrid_p) # Make sure basingrid has value 1 otherwise all precip will be 0
theprec = arcpy.sa.Times(prec_grid,basingrid_p)
precavg = arcpy.GetRasterProperties_management(maprecbasin,"MEAN")
precavg = float(precavg.getOutput(0))
avgprec = arcpy.GetRasterProperties_management(theprec,"MEAN")
avgprec = float(avgprec.getOutput(0))
global maprec
maprec = float(precavg/(1000* 2.54))
global p2yr
p2yr = float(avgprec/1000)
#*******************************************************************************************************
# format precision before text file string settings
#*******************************************************************************************************
maxlength = "{0:.2f}".format(maxlength)
theslope = "{0:.1f}".format(theslope)
theslope_feet = "{0:.5f}".format(theslope_feet)
landslope = "{0:.3f}".format(landslope)
pctAsoil = "{0:.1f}".format(pctAsoil)
pctBsoil = "{0:.1f}".format(pctBsoil)
pctCsoil = "{0:.1f}".format(pctCsoil)
pctDsoil = "{0:.1f}".format(pctDsoil)
pctWsoil = "{0:.1f}".format(pctWsoil)
UrbPct = "{0:.1f}".format(UrbPct)
FC = "{0:.1f}".format(FC)
ST = "{0:.1f}".format(ST)
IA = "{0:.1f}".format(IA)
LI = "{0:.1f}".format(LI)
areami2 = "{0:.2f}".format(areami2)
basinrelief = "{0:.2f}".format(basinrelief)
avgCN = "{0:.1f}".format(avgCN)
pctAsoilR = "{0:.1f}".format(pctAsoilR)
pctBsoilR = "{0:.1f}".format(pctBsoilR)
pctCsoilR = "{0:.1f}".format(pctCsoilR)
pctDsoilR = "{0:.1f}".format(pctDsoilR)
tc = "{0:.2f}".format(tc)
lagtime = "{0:.2f}".format(lagtime)
maprec = "{0:.2f}".format(maprec)
p2yr = "{0:.2f}".format(p2yr)
#*******************************************************************************************************
# Basin statistics analysis date
#*******************************************************************************************************
now = datetime.now()
month = now.strftime("%B")
day = now.strftime("%d")
year = now.strftime("%Y")
#*******************************************************************************************************
# Text file string variables
#*******************************************************************************************************
datastring = ""
datastring = datastring + "GISHydro Release Version Date: %s" "\n" %(Modifieddt)
datastring = datastring + "Project Name: %s" %(proj)
datastring = datastring + "" "\n"
datastring = datastring + "Analysis Date: %s %s, %s " "\n" %(month,day,year)
datastring = datastring + "Data Selected:" "\n"
datastring = datastring + " DEM Coverage: %s" "\n" %(demtot)
datastring = datastring + " Land Use Coverage: %s" "\n" %(landuse)
datastring = datastring + " Soil Coverage: %s" "\n" %(soil)
datastring = datastring + " Hydrologic Condition: %s" "\n" %(hyd)
datastring = datastring + " Impose NHD stream Locations: Yes" "\n"
datastring = datastring + " Outlet Easting: %s m (MD Stateplane, NAD 1983)" "\n" %(int(xoutlet))
datastring = datastring + " Outlet Northing: %s m (MD Stateplane, NAD 1983)" "\n" %(int(youtlet))
datastring = datastring + "Findings: " "\n"
datastring = datastring + " Outlet Location: %s" "\n" %(provstring)
datastring = datastring + " Outlet State: Maryland" "\n"
datastring = datastring + " Drainage Area %s square miles" "\n" %(areami2)
datastring = datastring + breakstring
datastring = datastring + "" "\n"
datastring = datastring + " Channel Slope: %s feet/mile (%s feet/feet)" "\n" %(theslope,theslope_feet)
datastring = datastring + " Land Slope: %s feet/feet" "\n" %(landslope)
datastring = datastring + " Urban Area (percent): %s " "\n" %(UrbPct)
datastring = datastring + " Impervious Area (percent): %s " "\n" %(IA)
#*******************************************************************************************************
# Print out Impervious area warning message
# *** warning message is included in for loop despite the fact that technically it could be printed
# twice. Since both "Appalachian Plateau" and "Eastern Coastal Plain" are far apart so it is
# impossible to have that big watershed while doing analysis with GISHydroNXT
#*******************************************************************************************************
theVTab = arcpy.SearchCursor("theVTab","","","Province","")
for row in theVTab:
if (row.getValue("Province") == "A") or (row.getValue("Province") == "E"):
if float(IA) >= 10:
pythonaddins.MessageBox("IMPERVIOUS AREA IN WATERSHED EXCEEDS 10%." "\n"
"Calculated discharges from USGS Regression" "\n"
"Equations may not be appropriate.", "Warning", 0)
datastring = datastring + Impwarntext
else:
datastring = datastring + ""
datastring = datastring + ""
#*******************************************************************************************************
# Close to boundary condition for provinces -- Near tool isn"t available with
# basic level license therefore a more crude method was emplyed here. It can
# be improved in future by using "arcpy.Geometry()" tool to get distance
#*******************************************************************************************************
arcpy.Buffer_analysis(optfolder + "/watershed.shp",optfolder + "/wats_prov.shp","5000","#","#","ALL","FID")
province = r"" + Directory + "/data/maryland/provlines.shp"
wats_prov = optfolder + "/wats_prov.shp"
prov_int = optfolder + "/prov_int"
arcpy.Intersect_analysis([province,wats_prov],prov_int,"ALL","#","INPUT")
prov_cursor = arcpy.SearchCursor(optfolder + "/prov_int.shp","","","FID","")
prov = prov_cursor.next()
if prov != None:
datastring = datastring + provwarntext
#*******************************************************************************************************
# Close to boundary condition for limestone -- Near tool isn"t available with
# basic level license therefore a more crude method was emplyed here. It can
# be improved in future by using "arcpy.Geometry()" tool to get distance
#*******************************************************************************************************
arcpy.Buffer_analysis(optfolder + "/watershed.shp",optfolder + "/wats_lime.shp","1000","#","#","ALL","FID")
wats_lime = optfolder + "/wats_lime.shp"
lime_int = optfolder + "/lime_int"
arcpy.Intersect_analysis([limestonem,wats_lime],lime_int,"ALL","#","INPUT")
lime_cursor = arcpy.SearchCursor(optfolder + "/lime_int.shp","","","FID","")
lime = lime_cursor.next()
if lime != None:
datastring = datastring + limewarntext
#*******************************************************************************************************
# Continued -- Text file string variables
#*******************************************************************************************************
datastring = datastring + "" "\n"
datastring = datastring + " Time of Concentration: %s hours [W.O. Thomas, Jr. Equation]" "\n" %(tc)
datastring = datastring + " Time of Concentration: %s hours [From SCS Lag Equation * 1.67]" "\n" %(lagtime)
datastring = datastring + " Longest Flow Path: %s miles" "\n" %(maxlength)
datastring = datastring + " Basin RelieF: %s feet" "\n" %(basinrelief)
datastring = datastring + " Average CN: %s" "\n" %(avgCN)
datastring = datastring + " Forest Cover (percent): %s" "\n" %(FC)
datastring = datastring + " Storage (percent): %s" "\n" %(ST)
datastring = datastring + " Limestone (percent): %s" "\n" %(LI)
datastring = datastring + " Selected Soils Data Statistics Percent:" "\n"
datastring = datastring + " A Soils: %s" "\n" %(pctAsoilR)
datastring = datastring + " B Soils: %s" "\n" %(pctBsoilR)
datastring = datastring + " C Soils: %s" "\n" %(pctCsoilR)
datastring = datastring + " D Soils: %s" "\n" %(pctDsoilR)
datastring = datastring + " SSURGO Soils Data Statistics Percent (used in Regression Equations):" "\n"
datastring = datastring + " A Soils: %s" "\n" %(pctAsoil)
datastring = datastring + " B Soils: %s" "\n" %(pctBsoil)
datastring = datastring + " C Soils: %s" "\n" %(pctCsoil)
datastring = datastring + " D Soils: %s" "\n" %(pctDsoil)
datastring = datastring + " 2-Year,24-hour Prec.: %s inches" "\n" %(p2yr)
datastring = datastring + " Mean Annual Prec.: %s inches" "\n" %(maprec)
#*******************************************************************************************************
# turn layers ON/OFF in current data frame
#*******************************************************************************************************
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]
layers = arcpy.mapping.ListLayers(mxd, "", df)
for lyr in layers:
if lyr.name == "landslope":
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "ssurgo":
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "limegrid":
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "wats_prov":
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "prov_int":
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "wats_lime":
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "lime_int":
arcpy.mapping.RemoveLayer(df, lyr)
arcpy.RefreshTOC()
arcpy.RefreshActiveView()
#*******************************************************************************************************
# write strings to basin stat text file.
# Message box containing datastring as message
#*******************************************************************************************************
defFN = optfolder + "/basinstat.txt"
statfile = open(defFN,"w")
statfile.write(datastring)
statfile.close()
#*******************************************************************************************************
# open "basinstat" file in text editor
#*******************************************************************************************************
hydro.openbrowser(defFN)
#*******************************************************************************************************
# turn peak discharge OFF
#*******************************************************************************************************
button1.enabled = True
button4.enabled = True
arcpy.env.extent = "MAXOF"