GISHydroNXT System Documentation
#*********************************************************************************
# Author: UMD
# Date: 24-07-2018
# Modified: n/a
# Classes: WriteSubAreaLandUseDistribution()
# Functions: hydro.lu_description() ; hydro.openbrowser()
# Modules: arcpy
# Comments: n/a
#*********************************************************************************
class WriteSubAreaLandUseDistribution(object):
"""Implementation for GISHydroNXT_addin.button12 (Button)"""
def __init__(self):
self.enabled = False
self.checked = False
def onClick(self):
arcpy.env.scratchWorkspace = scratchfolder
arcpy.env.workspace = optfolder
# read datetime time stamp single time (so it should be here outside loop)
now = datetime.now()
month = now.strftime("%B")
day = now.strftime("%d")
year = now.strftime("%Y")
# write initial lines of sub-basin composition outside for loop
# (concatenate all sub-areas basin compositions)
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 + "" "\n"
# create a folder for basin composition files
sub_basincomp = optfolder + "/sub_basincomp"
if not os.path.exists(sub_basincomp):
os.makedirs(sub_basincomp)
fc = optfolder + "/subshed.shp"
desc = arcpy.Describe(fc)
rows = arcpy.SearchCursor(fc)
for idx,row in enumerate(rows):
aPoly = row.getValue(desc.shapefieldname)
setExtent = aPoly.extent
arcpy.env.extent = setExtent
# create mask raster for each sub-basin
basingrid = optfolder + "/basingrid"
arcpy.MakeFeatureLayer_management(fc, "layer" + str(idx), ' "FID" = ' + str(idx))
mask = arcpy.Clip_management(basingrid, "", optfolder + "/aux_folder/mask" + str(idx), "layer" + str(idx),"0", "ClippingGeometry")
## grids and shapefile used for sub-basin land distribution
lu = optfolder + "/landuse"
soil = optfolder + "/Soils"
nhd = r"" + Directory + "/data/maryland/nhd_streamsm.shp"
lu_ext = optfolder + "/sub_basincomp/lu_ext"
soil_ext = optfolder + "/sub_basincomp/soil_ext"
lu_out = arcpy.sa.Times(lu,mask)
lu_out.save(lu_ext + str(idx))
soil_out = arcpy.sa.Times(soil,mask)
soil_out.save(soil_ext + str(idx))
# 1] intersect sub-basins and NHD streams to obtain string of stream names in each sub-basin
# 2] obtain name strings
arcpy.env.addOutputsToMap = False
nhd_sub = optfolder + "/sub_basincomp/nhd" + str(idx) + ".shp"
arcpy.Intersect_analysis([nhd,"layer" + str(idx)],nhd_sub,"ALL","#","INPUT")
nhd_sc = arcpy.SearchCursor(nhd_sub,"","","NAME","")
nhd_streams = []
for s in nhd_sc:
nhd_names = s.getValue("NAME")
nhd_streams.append(nhd_names)
# unicode to string conversion and then replace " " in list to "Unidentified Reach"
nhd_streams_tmp = [x.encode("UTF8") for x in nhd_streams]
nhd_streams = ["Unidentified Reach" if i == " " else i for i in nhd_streams_tmp] # replace "" in a list
## convert all clipped rasters to polygons
lu_poly = optfolder + "/sub_basincomp/lu_poly" + str(idx) + ".shp"
soil_poly = optfolder + "/sub_basincomp/soil_poly" + str(idx) + ".shp"
arcpy.env.addOutputsToMap = False
arcpy.RasterToPolygon_conversion(lu_ext + str(idx),lu_poly,"NO_SIMPLIFY","VALUE")
arcpy.env.addOutputsToMap = False
arcpy.RasterToPolygon_conversion(soil_ext + str(idx),soil_poly,"NO_SIMPLIFY","VALUE")
## intersect land use and soil to prepare two polygons: "lu_soil" and "lu_cn"
arcpy.env.addOutputsToMap = False
lu_soil = optfolder + "/sub_basincomp/lu_soil" + str(idx) + ".shp"
arcpy.Intersect_analysis([lu_poly,soil_poly],lu_soil,"ALL","#","INPUT")
## dissolve above intersected polygon
arcpy.env.addOutputsToMap = False
lu_soil_diss = optfolder + "/sub_basincomp/lu_soil_diss" + str(idx) + ".shp"
arcpy.Dissolve_management(lu_soil,lu_soil_diss,"GRIDCODE;GRIDCODE_1","#","MULTI_PART","DISSOLVE_LINES")
## add filed to both of above dissolved polygons and compute area in acres
arcpy.AddField_management(lu_soil_diss,"area","FLOAT",15,4)
arcpy.CalculateField_management(lu_soil_diss,"area","!shape.area@acres!","PYTHON")
# prepre a list of lu codes to feed into lu_description function in order to obtain matching descriptions list
lu_match = []
sc = arcpy.SearchCursor(lu_ext + str(idx),"","","VALUE","")
for i in sc:
v = i.getValue("VALUE")
lu_match.append(v)
# create list of lists with zeroes
soil_acre_lists = [[0,0,0,0] for i in range(len(lu_match))]
# preapre a list of soil acreage using lu_match list
lc_soil_diss = []
soil_lc_diss = []
lc_soil_aa = []
sr = arcpy.SearchCursor(lu_soil_diss,"","","GRIDCODE;GRIDCODE_1;area","")
for s in sr:
lc = s.getValue("GRIDCODE")
lc_soil_diss.append(lc)
sc = s.getValue("GRIDCODE_1")
soil_lc_diss.append(sc)
aa = s.getValue("area")
lc_soil_aa.append(round(aa,2))
for n,lu in enumerate(lu_match):
for l,s,a in zip(lc_soil_diss,soil_lc_diss,lc_soil_aa):
if l == lu:
soil_acre_lists[n][int(s)-1] = a
# prepare matching list of lu description using lu codes from lu raster of watershed
if landuse == "NLCD 2011 Landuse":
if hyd == "Fair":
lut_file = r"" + Directory + "/data/mdfiles/lookup/nlcdlookupfair.txt"
elif hyd == "Good":
lut_file = r"" + Directory + "/data/mdfiles/lookup/nlcdlookupgood.txt"
elif hyd == "Poor":
lut_file = r"" + Directory + "/data/mdfiles/lookup/nlcdlookuppoor.txt"
if landuse == "NLCD 2006 Landuse":
if hyd == "Fair":
lut_file = r"" + Directory + "/data/mdfiles/lookup/nlcdlookupfair.txt"
elif hyd == "Good":
lut_file = r"" + Directory + "/data/mdfiles/lookup/nlcdlookupgood.txt"
elif hyd == "Poor":
lut_file = r"" + Directory + "/data/mdfiles/lookup/nlcdlookuppoor.txt"
if landuse == "NLCD 2001 Landuse":
if hyd == "Fair":
lut_file = r"" + Directory + "/data/mdfiles/lookup/nlcdlookupfair.txt"
elif hyd == "Good":
lut_file = r"" + Directory + "/data/mdfiles/lookup/nlcdlookupgood.txt"
elif hyd == "Poor":
lut_file = r"" + Directory + "/data/mdfiles/lookup/nlcdlookuppoor.txt"
if landuse == "1997 MOP Landuse":
if hyd == "Fair":
lut_file = r"" + Directory + "/data/mdfiles/lookup/andlookupfair.txt"
elif hyd == "Good":
lut_file = r"" + Directory + "/data/mdfiles/lookup/andlookupgood.txt"
elif hyd == "Poor":
lut_file = r"" + Directory + "/data/mdfiles/lookup/andlookuppoor.txt"
if landuse == "2002 MOP Landuse":
if hyd == "Fair":
lut_file = r"" + Directory + "/data/mdfiles/lookup/andlookupfair.txt"
elif hyd == "Good":
lut_file = r"" + Directory + "/data/mdfiles/lookup/andlookupgood.txt"
elif hyd == "Poor":
lut_file = r"" + Directory + "/data/mdfiles/lookup/andlookuppoor.txt"
if landuse == "2010 MOP Landuse":
if hyd == "Fair":
lut_file = r"" + Directory + "/data/mdfiles/lookup/andlookupfair.txt"
elif hyd == "Good":
lut_file = r"" + Directory + "/data/mdfiles/lookup/andlookupgood.txt"
elif hyd == "Poor":
lut_file = r"" + Directory + "/data/mdfiles/lookup/andlookuppoor.txt"
if landuse == "2002 MD/DE Landuse":
if hyd == "Fair":
lut_file = r"" + Directory + "/data/mdfiles/lookup/mddelookupfair.txt"
elif hyd == "Good":
lut_file = r"" + Directory + "/data/mdfiles/lookup/mddelookupgood.txt"
elif hyd == "Poor":
lut_file = r"" + Directory + "/data/mdfiles/lookup/mddelookuppoor.txt"
if landuse == "Ultimate Landuse":
if hyd == "Fair":
lut_file = r"" + Directory + "/data/mdfiles/lookup/zoninglookupfair.txt"
elif hyd == "Good":
lut_file = r"" + Directory + "/data/mdfiles/lookup/zoninglookupgood.txt"
elif hyd == "Poor":
lut_file = r"" + Directory + "/data/mdfiles/lookup/zoninglookuppoor.txt"
if landuse == "MRLC Landuse":
if hyd == "Fair":
lut_file = r"" + Directory + "/data/mdfiles/lookup/mrlclookupfair.txt"
elif hyd == "Good":
lut_file = r"" + Directory + "/data/mdfiles/lookup/mrlclookupgood.txt"
elif hyd == "Poor":
lut_file = r"" + Directory + "/data/mdfiles/lookup/mrlclookuppoor.txt"
if landuse == "1970s USGS Landuse":
if hyd == "Fair":
lut_file = r"" + Directory + "/data/mdfiles/lookup/usgslookupfair.txt"
elif hyd == "Good":
lut_file = r"" + Directory + "/data/mdfiles/lookup/usgslookupgood.txt"
elif hyd == "Poor":
lut_file = r"" + Directory + "/data/mdfiles/lookup/usgslookuppoor.txt"
# run hydro function to obtain land use description of categories present in watershed
lu_desc = hydro.lu_description(lut_file,lu_match)
# Perform two tasks:
# a) take "lu_desc" and "soil_acre_lists" and concatenate them
# b) format according to basin composition file in legacy version
land_soil_area = ""
width=30
for ld,sg in zip(lu_desc,soil_acre_lists):
land_soil_area = land_soil_area + "{: <{}}".format(ld, width) + str(sg[0]).rjust(10) + str(sg[1]).rjust(10) + str(sg[2]).rjust(10) + str(sg[3]).rjust(10)
land_soil_area = land_soil_area + "" "\n"
#*******************************************************************************************************
# Text file string variables
#*******************************************************************************************************
datastring = datastring + "" "\n"
datastring = datastring + "Landuse and Soil Distributions for: Sub-Area %s" "\n" %(str(idx+1)) # 1/2/2018: because we want sub-basin index to start from 1 so "idx+1"
datastring = datastring + "" "\n"
datastring = datastring + "Streams located in this sub-area:"
datastring = datastring + "" "\n"
# add stream names
for sn in nhd_streams:
datastring = datastring + str(sn)
datastring = datastring + "" "\n"
datastring = datastring + "" "\n"
datastring = datastring + "Distribution of Landuse by Soil Group" "\n"
datastring = datastring + "" "\n"
datastring = datastring + "Acres on Indicated Soil Group".rjust(66)
datastring = datastring + "" "\n"
datastring = datastring + "Land Use".rjust(8) + "A-Soil".rjust(32) + "B-Soil".rjust(10) + "C-Soil".rjust(10) + "D-Soil".rjust(10)
datastring = datastring + "" "\n"
datastring = datastring + "" "\n"
datastring = datastring + land_soil_area
# sum list of lists separately and cat at the end of lu description
total_area = [sum(i) for i in zip(*soil_acre_lists)]
datastring = datastring + "{: <{}}".format("Total Area:", width) + str(total_area[0]).rjust(10) + str(total_area[1]).rjust(10) + str(total_area[2]).rjust(10) + str(total_area[3]).rjust(10)
datastring = datastring + "" "\n"
datastring = datastring + "" "\n"
datastring = datastring + "Distribution of Land Use and Curve Numbers Used" "\n"
datastring = datastring + "" "\n"
datastring = datastring + "Land Use".rjust(8) + "Acres".rjust(34) + "Percent".rjust(10) + "A".rjust(4) + "B".rjust(4) + "C".rjust(4) + "D".rjust(4)
datastring = datastring + "" "\n"
datastring = datastring + "" "\n"
# loop over land use, related total acreage, percent of land covered by this lu category, and A-B-C-D curve numbers
curve_num = []
for l in lu_match:
with open(lut_file, "r") as f:
next(f)
for line in f:
luc = line.split("\t")[0]
if int(l) == int(luc):
temp=[]
A = line.split("\t")[2] # CN A
temp.append(A)
B = line.split("\t")[3] # CN B
temp.append(B)
C = line.split("\t")[4] # CN C
temp.append(C)
D = line.split("\t")[5] # CN D
temp.append(D)
curve_num.append(temp)
# sum areas for each sub-list individually
acres = [sum(i) for i in soil_acre_lists]
total_all = sum(total_area)
percent = [round(float(x/total_all)*100,2) for x in acres]
for l,a,p,cn in zip(lu_desc,acres,percent,curve_num):
datastring = datastring + "{: <{}}".format(l, width) + str(a).rjust(12) + str(p).rjust(10) + str(cn[0]).rjust(4) + str(cn[1]).rjust(4) + str(cn[2]).rjust(4) + str(cn[3]).rjust(4)
datastring = datastring + "" "\n"
#*******************************************************************************************************
# 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 == "lu_ext":
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "soil_ext":
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "lu_poly":
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "soil_poly":
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "lu_soil":
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "lu_soil_diss":
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "layer" + str(idx):
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "mask" + str(idx):
arcpy.mapping.RemoveLayer(df, lyr)
arcpy.RefreshTOC()
arcpy.RefreshActiveView()
#*******************************************************************************************************
# write strings to basin stat text file.
# Message box containing datastring as message
#*******************************************************************************************************
defFN = optfolder + "/sub_basincomp.txt"
compfile = open(defFN,"w")
compfile.write(datastring)
compfile.close()
#*******************************************************************************************************
# open "basincomp" file in text editor
#*******************************************************************************************************
hydro.openbrowser(defFN)
#*******************************************************************************************************
# turn peak discharge OFF
#*******************************************************************************************************
button12.enabled = False
arcpy.env.extent = "MAXOF"