GISHydroNXT System Documentation
"""
This is hydro library containing necessary functions to compute basin
statistics and variables as required by hydro menu in GISHydro2000. It also use global
variable definition to populate data into memory for access across functions and scripts.
"""
#*******************************************************************************************************
# Author: UMD
# Date: 24-07-2018
# Modified: n/a
# Purpose: To house all required functions
# Calls: None
#*******************************************************************************************************
import arcpy
from arcpy import env
import os
import shutil
from arcpy.sa import *
from datetime import datetime
import numpy as np
from numpy import copy
import time
import pythonaddins
import functools
import threading
import webbrowser
#*******************************************************************************************************
# environmental settings are only for statsgo data -- it can"t create grids in memory and therefore
# are required to be referenced to the current working directory for function execution
#*******************************************************************************************************
arcpy.env.overwriteOutput = True
#*******************************************************************************************************
# Calculate channel slope
#*******************************************************************************************************
def channelslope(dirgrid,elevgrid,opath):
arcpy.env.workspace = opath
arcpy.env.scratchWorkspace = opath + "/aux_folder"
arcpy.env.addOutputsToMap = False
global maxlength
# check names and number of input files as function input arguments carefully
dGrid = arcpy.sa.FlowLength(dirgrid, "DOWNSTREAM", "") # saved all rasters at opath, to extract raster properties
uGrid = arcpy.sa.FlowLength(dirgrid, "UPSTREAM", "")
sGrid = dGrid + uGrid
sGrid.save(opath + "/maxlength")
maxlen = arcpy.GetRasterProperties_management(uGrid, "MAXIMUM") # 05/02/2013: changed name from "maxlength" to "maxlen
maxlen = float(maxlen.getOutput(0))
desc = arcpy.Describe(dirgrid)
cellSizeVal = desc.meanCellWidth
q = cellSizeVal
tolerance = (0.1)* q
Lpath = arcpy.sa.Con(sGrid > (maxlen - tolerance), arcpy.sa.Minus(1,sGrid))
Lpath = arcpy.sa.IsNull(Lpath)
long_path = arcpy.sa.Con(Lpath == 0, 1, 0)
lpath_upgrid = arcpy.sa.Times(long_path,uGrid)
maxlen_09 = float(0.9*maxlen)
maxlen_15 = float(0.15*maxlen)
lpath_09 = arcpy.sa.Con(lpath_upgrid < maxlen_09, lpath_upgrid, 0)
lpath_1085 = arcpy.sa.Con(lpath_09 > maxlen_15, lpath_09, 0)
long_1085 = arcpy.sa.SetNull(lpath_1085,1,"VALUE = 0")
long_elev = arcpy.sa.Times(long_1085,elevgrid)
maxlength = float(maxlen* (0.000621371))
elevmin = arcpy.GetRasterProperties_management(long_elev, "MINIMUM")
elevmax = arcpy.GetRasterProperties_management(long_elev, "MAXIMUM")
min_elev = float(elevmin.getOutput(0))
max_elev = float(elevmax.getOutput(0))
channelslope = float((max_elev) - (min_elev))/(maxlength* (0.75)) # units: xxx feet/mile
return channelslope
#*******************************************************************************************************
# Calculate land slope
#*******************************************************************************************************
def slopegrid(flowdir,elevgrid,src,opath): # what is the use of basingrid in this function?
arcpy.env.workspace = opath
arcpy.env.scratchWorkspace = opath + "/aux_folder"
arcpy.env.addOutputsToMap = False
arcpy.env.extent = "flowdir"
path = opath + "/slope_calc" # changed optfolder to opath in following lines as well
os.mkdir(path, 0755)
dirfile = "temp1.asc"
demfile = "temp2.asc"
slopefile = "temp3.asc"
thegisdir = opath + "/slope_calc/calcslope.txt"
datafilename = open(thegisdir, "w")
datafilename.write(str(demfile) + "\n" + str(dirfile) + "\n" + str(slopefile))
datafilename.close()
dirGrid = arcpy.sa.Reclassify("flowdir","Value",
"1 1;2 8;4 7;8 6;16 5;32 4;64 3;128 2","DATA")
arcpy.RasterToASCII_conversion(dirGrid, opath + "/slope_calc/temp1.asc")
arcpy.RasterToASCII_conversion(elevgrid, opath + "/slope_calc/temp2.asc")
dst = opath + "/slope_calc"
shutil.copy2(src, dst)
os.chdir(opath + "/slope_calc")
os.system(opath + "/slope_calc/calcslope.exe")
time.sleep(4)
arcpy.ASCIIToRaster_conversion(opath + "/slope_calc/temp3.asc",
opath + "/slope_calc/landslope", "FLOAT") # modified to get float output
slopegrid = arcpy.sa.Con(opath + "/slope_calc/landslope" > 0,
opath + "/slope_calc/landslope", 0) # 06-03-2013: slope value was going below
# 0 which isn"t correct. Condition is
# added to at least have 0.01
return slopegrid
#*******************************************************************************************************
# Calculate curve number for fair soil condition using RAGAN or SSURGO soil data using NLCD tables
#*******************************************************************************************************
def nlcdlookupfair(inLU,inSoil):
lu = arcpy.RasterToNumPyArray(inLU)
sc = arcpy.RasterToNumPyArray(inSoil)
desc = arcpy.Describe(inLU)
cellSize = desc.meanCellHeight
extent = desc.Extent
pnt = arcpy.Point(extent.XMin,extent.YMin)
table = {(11,1):100,(11,2):100,(11,3):100,(11,4):100,
(21,1):49,(21,2):69,(21,3):79,(21,4):84,
(22,1):61,(22,2):76,(22,3):84,(22,4):88,
(23,1):68,(23,2):80,(23,3):86,(23,4):89,
(24,1):81,(24,2):88,(24,3):91,(24,4):93,
(31,1):77,(31,2):86,(31,3):91,(31,4):94,
(41,1):36,(41,2):60,(41,3):73,(41,4):79,
(42,1):36,(42,2):60,(42,3):73,(42,4):79,
(43,1):36,(43,2):60,(43,3):73,(43,4):79,
(52,1):35,(52,2):56,(52,3):70,(52,4):77,
(71,1):49,(71,2):69,(71,3):79,(71,4):84,
(81,1):49,(81,2):69,(81,3):79,(81,4):84,
(82,1):70,(82,2):80,(82,3):87,(82,4):90,
(90,1):100,(90,2):100,(90,3):100,(90,4):100,
(95,1):100,(95,2):100,(95,3):100,(95,4):100}
curveNumber = np.ones(lu.shape, dtype = np.int)
for k,v in table.iteritems():
curveNumber[((lu == k[0]) & (sc == k[1]))] = v
cnArray = arcpy.NumPyArrayToRaster(curveNumber,pnt,cellSize,cellSize)
return cnArray
def nlcdlookupgood(inLU,inSoil):
lu = arcpy.RasterToNumPyArray(inLU)
sc = arcpy.RasterToNumPyArray(inSoil)
desc = arcpy.Describe(inLU)
cellSize = desc.meanCellHeight
extent = desc.Extent
pnt = arcpy.Point(extent.XMin,extent.YMin)
table = {(11,1):100,(11,2):100,(11,3):100,(11,4):100,
(21,1):39,(21,2):61,(21,3):74,(21,4):80,
(22,1):54,(22,2):70,(22,3):80,(22,4):85,
(23,1):61,(23,2):75,(23,3):83,(23,4):87,
(24,1):77,(24,2):85,(24,3):90,(24,4):92,
(31,1):77,(31,2):86,(31,3):91,(31,4):94,
(41,1):30,(41,2):55,(41,3):70,(41,4):77,
(42,1):30,(42,2):55,(42,3):70,(42,4):77,
(43,1):30,(43,2):55,(43,3):70,(43,4):77,
(52,1):30,(52,2):48,(52,3):65,(52,4):73,
(71,1):39,(71,2):61,(71,3):74,(71,4):80,
(81,1):39,(81,2):61,(81,3):74,(81,4):80,
(82,1):67,(82,2):78,(82,3):85,(82,4):89,
(90,1):100,(90,2):100,(90,3):100,(90,4):100,
(95,1):100,(95,2):100,(95,3):100,(95,4):100}
curveNumber = np.ones(lu.shape, dtype = np.int)
for k,v in table.iteritems():
curveNumber[((lu == k[0]) & (sc == k[1]))] = v
cnArray = arcpy.NumPyArrayToRaster(curveNumber,pnt,cellSize,cellSize)
return cnArray
def nlcdlookuppoor(inLU,inSoil):
lu = arcpy.RasterToNumPyArray(inLU)
sc = arcpy.RasterToNumPyArray(inSoil)
desc = arcpy.Describe(inLU)
cellSize = desc.meanCellHeight
extent = desc.Extent
pnt = arcpy.Point(extent.XMin,extent.YMin)
table = {(11,1):100,(11,2):100,(11,3):100,(11,4):100,
(21,1):68,(21,2):79,(21,3):86,(21,4):89,
(22,1):76,(22,2):84,(22,3):89,(22,4):91,
(23,1):79,(23,2):86,(23,3):91,(23,4):92,
(24,1):88,(24,2):91,(24,3):94,(24,4):95,
(31,1):77,(31,2):86,(31,3):91,(31,4):94,
(41,1):45,(41,2):66,(41,3):77,(41,4):83,
(42,1):45,(42,2):66,(42,3):77,(42,4):83,
(43,1):45,(43,2):66,(43,3):77,(43,4):83,
(52,1):48,(52,2):67,(52,3):77,(52,4):83,
(71,1):68,(71,2):79,(71,3):86,(71,4):89,
(81,1):68,(81,2):79,(81,3):86,(81,4):89,
(82,1):72,(82,2):81,(82,3):88,(82,4):91,
(90,1):100,(90,2):100,(90,3):100,(90,4):100,
(95,1):100,(95,2):100,(95,3):100,(95,4):100}
curveNumber = np.ones(lu.shape, dtype = np.int)
for k,v in table.iteritems():
curveNumber[((lu == k[0]) & (sc == k[1]))] = v
cnArray = arcpy.NumPyArrayToRaster(curveNumber,pnt,cellSize,cellSize)
return cnArray
#*******************************************************************************************************
# Calculate curve number for fair soil condition using RAGAN or SSURGO soil data using Anderson tables
#*******************************************************************************************************
def andlookupfair(inLU,inSoil):
lu = arcpy.RasterToNumPyArray(inLU)
sc = arcpy.RasterToNumPyArray(inSoil)
desc = arcpy.Describe(inLU)
cellSize = desc.meanCellHeight
extent = desc.Extent
pnt = arcpy.Point(extent.XMin,extent.YMin)
table = {(10,1):68,(10,2):80,(10,3):86,(10,4):89,
(11,1):61,(11,2):76,(11,3):84,(11,4):88,
(12,1):68,(12,2):80,(12,3):86,(12,4):89,
(13,1):81,(13,2):88,(13,3):91,(13,4):93,
(14,1):91,(14,2):94,(14,3):95,(14,4):96,
(15,1):81,(15,2):88,(15,3):91,(15,4):93,
(16,1):74,(16,2):84,(16,3):89,(16,4):91,
(17,1):77,(17,2):86,(17,3):91,(17,4):94,
(18,1):49,(18,2):69,(18,3):79,(18,4):84,
(20,1):70,(20,2):80,(20,3):87,(20,4):90,
(21,1):70,(21,2):80,(21,3):87,(21,4):90,
(22,1):49,(22,2):69,(22,3):79,(22,4):84,
(23,1):43,(23,2):65,(23,3):76,(23,4):82,
(24,1):59,(24,2):74,(24,3):82,(24,4):86,
(25,1):70,(25,2):80,(25,3):87,(25,4):90,
(40,1):36,(40,2):60,(40,3):73,(40,4):79,
(41,1):36,(41,2):60,(41,3):73,(41,4):79,
(42,1):36,(42,2):60,(42,3):73,(42,4):79,
(43,1):36,(43,2):60,(43,3):73,(43,4):79,
(44,1):35,(44,2):56,(44,3):70,(44,4):77,
(50,1):100,(50,2):100,(50,3):100,(50,4):100,
(60,1):100,(60,2):100,(60,3):100,(60,4):100,
(70,1):77,(70,2):86,(70,3):91,(70,4):94,
(71,1):77,(71,2):86,(71,3):91,(71,4):94,
(72,1):77,(72,2):86,(72,3):91,(72,4):94,
(73,1):77,(73,2):86,(73,3):91,(73,4):94,
(80,1):98,(80,2):98,(80,3):98,(80,4):98,
(191,1):70,(191,2):80,(191,3):87,(191,4):90,
(192,1):36,(192,2):60,(192,3):73,(192,4):79,
(241,1):59,(241,2):74,(241,3):82,(241,4):86,
(242,1):59,(242,2):74,(242,3):82,(242,4):86}
curveNumber = np.ones(lu.shape, dtype = np.int)
for k,v in table.iteritems():
curveNumber[((lu == k[0]) & (sc == k[1]))] = v
cnArray = arcpy.NumPyArrayToRaster(curveNumber,pnt,cellSize,cellSize)
return cnArray
def andlookupgood(inLU,inSoil):
lu = arcpy.RasterToNumPyArray(inLU)
sc = arcpy.RasterToNumPyArray(inSoil)
desc = arcpy.Describe(inLU)
cellSize = desc.meanCellHeight
extent = desc.Extent
pnt = arcpy.Point(extent.XMin,extent.YMin)
table = {(10,1):61,(10,2):75,(10,3):83,(10,4):87,
(11,1):54,(11,2):70,(11,3):80,(11,4):85,
(12,1):61,(12,2):75,(12,3):83,(12,4):87,
(13,1):77,(13,2):85,(13,3):90,(13,4):92,
(14,1):89,(14,2):92,(14,3):94,(14,4):95,
(15,1):81,(15,2):88,(15,3):91,(15,4):93,
(16,1):69,(16,2):80,(16,3):86,(16,4):89,
(17,1):77,(17,2):86,(17,3):91,(17,4):94,
(18,1):39,(18,2):61,(18,3):74,(18,4):80,
(20,1):67,(20,2):78,(20,3):85,(20,4):89,
(21,1):67,(21,2):78,(21,3):85,(21,4):89,
(22,1):39,(22,2):61,(22,3):74,(22,4):80,
(23,1):32,(23,2):58,(23,3):72,(23,4):79,
(24,1):59,(24,2):74,(24,3):82,(24,4):86,
(25,1):67,(25,2):78,(25,3):85,(25,4):89,
(40,1):30,(40,2):55,(40,3):70,(40,4):77,
(41,1):30,(41,2):55,(41,3):70,(41,4):77,
(42,1):30,(42,2):55,(42,3):70,(42,4):77,
(43,1):30,(43,2):55,(43,3):70,(43,4):77,
(44,1):30,(44,2):48,(44,3):65,(44,4):73,
(50,1):100,(50,2):100,(50,3):100,(50,4):100,
(60,1):100,(60,2):100,(60,3):100,(60,4):100,
(70,1):77,(70,2):86,(70,3):91,(70,4):94,
(71,1):77,(71,2):86,(71,3):91,(71,4):94,
(72,1):77,(72,2):86,(72,3):91,(72,4):94,
(73,1):77,(73,2):86,(73,3):91,(73,4):94,
(80,1):83,(80,2):89,(80,3):92,(80,4):94,
(191,1):67,(191,2):78,(191,3):85,(191,4):89,
(192,1):30,(192,2):55,(192,3):70,(192,4):77,
(241,1):59,(241,2):74,(241,3):82,(241,4):86,
(242,1):59,(242,2):74,(242,3):82,(242,4):86}
curveNumber = np.ones(lu.shape, dtype = np.int)
for k,v in table.iteritems():
curveNumber[((lu == k[0]) & (sc == k[1]))] = v
cnArray = arcpy.NumPyArrayToRaster(curveNumber,pnt,cellSize,cellSize)
return cnArray
def andlookuppoor(inLU,inSoil):
lu = arcpy.RasterToNumPyArray(inLU)
sc = arcpy.RasterToNumPyArray(inSoil)
desc = arcpy.Describe(inLU)
cellSize = desc.meanCellHeight
extent = desc.Extent
pnt = arcpy.Point(extent.XMin,extent.YMin)
table = {(10,1):79,(10,2):86,(10,3):91,(10,4):92,
(11,1):76,(11,2):84,(11,3):89,(11,4):91,
(12,1):79,(12,2):86,(12,3):91,(12,4):92,
(13,1):88,(13,2):81,(13,3):94,(13,4):95,
(14,1):94,(14,2):95,(14,3):96,(14,4):97,
(15,1):90,(15,2):93,(15,3):95,(15,4):95,
(16,1):83,(16,2):89,(16,3):92,(16,4):94,
(17,1):77,(17,2):86,(17,3):91,(17,4):94,
(18,1):68,(18,2):79,(18,3):86,(18,4):89,
(20,1):72,(20,2):81,(20,3):88,(20,4):91,
(21,1):72,(21,2):81,(21,3):88,(21,4):91,
(22,1):68,(22,2):79,(22,3):86,(22,4):89,
(23,1):57,(23,2):73,(23,3):82,(23,4):86,
(24,1):59,(24,2):74,(24,3):82,(24,4):86,
(25,1):72,(25,2):81,(25,3):88,(25,4):91,
(40,1):45,(40,2):66,(40,3):77,(40,4):83,
(41,1):45,(41,2):66,(41,3):77,(41,4):83,
(42,1):45,(42,2):66,(42,3):77,(42,4):83,
(43,1):45,(43,2):66,(43,3):77,(43,4):83,
(44,1):48,(44,2):67,(44,3):77,(44,4):83,
(50,1):100,(50,2):100,(50,3):100,(50,4):100,
(60,1):100,(60,2):100,(60,3):100,(60,4):100,
(70,1):77,(70,2):86,(70,3):91,(70,4):94,
(71,1):77,(71,2):86,(71,3):91,(71,4):94,
(72,1):77,(72,2):86,(72,3):91,(72,4):94,
(73,1):77,(73,2):86,(73,3):91,(73,4):94,
(80,1):98,(80,2):98,(80,3):98,(80,4):98,
(191,1):72,(191,2):81,(191,3):88,(191,4):91,
(192,1):45,(192,2):66,(192,3):77,(192,4):83,
(241,1):59,(241,2):74,(241,3):82,(241,4):86,
(242,1):59,(242,2):74,(242,3):82,(242,4):86}
curveNumber = np.ones(lu.shape, dtype = np.int)
for k,v in table.iteritems():
curveNumber[((lu == k[0]) & (sc == k[1]))] = v
cnArray = arcpy.NumPyArrayToRaster(curveNumber,pnt,cellSize,cellSize)
return cnArray
#*******************************************************************************************************
# Calculate curve number for poor soil condition using RAGAN or SSURGO soil data using MD/DE tables
#*******************************************************************************************************
def mddelookupfair(inLU,inSoil):
lu = arcpy.RasterToNumPyArray(inLU)
sc = arcpy.RasterToNumPyArray(inSoil)
desc = arcpy.Describe(inLU)
cellSize = desc.meanCellHeight
extent = desc.Extent
pnt = arcpy.Point(extent.XMin,extent.YMin)
table = {(10,1):68,(10,2):80,(10,3):86,(10,4):89,
(11,1):61,(11,2):76,(11,3):84,(11,4):88,
(12,1):68,(12,2):80,(12,3):86,(12,4):89,
(13,1):81,(13,2):88,(13,3):91,(13,4):93,
(14,1):91,(14,2):94,(14,3):95,(14,4):96,
(15,1):81,(15,2):88,(15,3):91,(15,4):93,
(16,1):74,(16,2):84,(16,3):89,(16,4):91,
(17,1):77,(17,2):86,(17,3):91,(17,4):94,
(18,1):49,(18,2):69,(18,3):79,(18,4):84,
(20,1):70,(20,2):80,(20,3):87,(20,4):90,
(21,1):70,(21,2):80,(21,3):87,(21,4):90,
(22,1):49,(22,2):69,(22,3):79,(22,4):84,
(23,1):43,(23,2):65,(23,3):76,(23,4):82,
(24,1):59,(24,2):74,(24,3):82,(24,4):86,
(25,1):70,(25,2):80,(25,3):87,(25,4):90,
(40,1):36,(40,2):60,(40,3):73,(40,4):79,
(41,1):36,(41,2):60,(41,3):73,(41,4):79,
(42,1):36,(42,2):60,(42,3):73,(42,4):79,
(43,1):36,(43,2):60,(43,3):73,(43,4):79,
(44,1):35,(44,2):56,(44,3):70,(44,4):77,
(50,1):100,(50,2):100,(50,3):100,(50,4):100,
(60,1):100,(60,2):100,(60,3):100,(60,4):100,
(70,1):77,(70,2):86,(70,3):91,(70,4):94,
(71,1):77,(71,2):86,(71,3):91,(71,4):94,
(72,1):77,(72,2):86,(72,3):91,(72,4):94,
(73,1):77,(73,2):86,(73,3):91,(73,4):94,
(80,1):98,(80,2):98,(80,3):98,(80,4):98,
(191,1):70,(191,2):80,(191,3):87,(191,4):90,
(192,1):36,(192,2):60,(192,3):73,(192,4):79,
(241,1):59,(241,2):74,(241,3):82,(241,4):86,
(242,1):59,(242,2):74,(242,3):82,(242,4):86,
(1110,1):68,(1110,2):80,(1110,3):86,(1110,4):89,
(1120,1):81,(1120,2):88,(1120,3):91,(1120,4):93,
(1140,1):81,(1140,2):88,(1140,3):91,(1140,4):93,
(1200,1):91,(1200,2):94,(1200,3):95,(1200,4):96,
(1210,1):91,(1210,2):94,(1210,3):95,(1210,4):96,
(1220,1):91,(1220,2):94,(1220,3):95,(1220,4):96,
(1230,1):91,(1230,2):94,(1230,3):95,(1230,4):96,
(1250,1):91,(1250,2):94,(1250,3):95,(1250,4):96,
(1290,1):91,(1290,2):94,(1290,3):95,(1290,4):96,
(1300,1):81,(1300,2):88,(1300,3):91,(1300,4):93,
(1400,1):98,(1400,2):98,(1400,3):98,(1400,4):98,
(1410,1):98,(1410,2):98,(1410,3):98,(1410,4):98,
(1420,1):91,(1420,2):94,(1420,3):95,(1420,4):96,
(1430,1):98,(1430,2):98,(1430,3):98,(1430,4):98,
(1440,1):98,(1440,2):98,(1440,3):98,(1440,4):98,
(1450,1):49,(1450,2):69,(1450,3):79,(1450,4):84,
(1460,1):98,(1460,2):98,(1460,3):98,(1460,4):98,
(1490,1):98,(1490,2):98,(1490,3):98,(1490,4):98,
(1500,1):49,(1500,2):69,(1500,3):79,(1500,4):84,
(1600,1):74,(1600,2):84,(1600,3):89,(1600,4):91,
(1700,1):74,(1700,2):84,(1700,3):89,(1700,4):91,
(1800,1):74,(1800,2):84,(1800,3):89,(1800,4):91,
(1900,1):49,(1900,2):69,(1900,3):79,(1900,4):84,
(2100,1):70,(2100,2):80,(2100,3):87,(2100,4):90,
(2110,1):70,(2110,2):80,(2110,3):87,(2110,4):90,
(2120,1):49,(2120,2):69,(2120,3):79,(2120,4):84,
(2130,1):70,(2130,2):80,(2130,3):87,(2130,4):90,
(2140,1):70,(2140,2):80,(2140,3):87,(2140,4):90,
(2150,1):70,(2150,2):80,(2150,3):87,(2150,4):90,
(2200,1):43,(2200,2):65,(2200,3):76,(2200,4):82,
(2300,1):59,(2300,2):74,(2300,3):82,(2300,4):86,
(2400,1):59,(2400,2):74,(2400,3):82,(2400,4):86,
(2900,1):70,(2900,2):80,(2900,3):87,(2900,4):90,
(3100,1):36,(3100,2):60,(3100,3):73,(3100,4):79,
(3200,1):36,(3200,2):60,(3200,3):73,(3200,4):79,
(3300,1):36,(3300,2):60,(3300,3):73,(3300,4):79,
(4100,1):36,(4100,2):60,(4100,3):73,(4100,4):79,
(4200,1):36,(4200,2):60,(4200,3):73,(4200,4):79,
(4300,1):36,(4300,2):60,(4300,3):73,(4300,4):79,
(4400,1):36,(4400,2):60,(4400,3):73,(4400,4):79,
(5100,1):100,(5100,2):100,(5100,3):100,(5100,4):100,
(5200,1):100,(5200,2):100,(5200,3):100,(5200,4):100,
(5300,1):100,(5300,2):100,(5300,3):100,(5300,4):100,
(5400,1):100,(5400,2):100,(5400,3):100,(5400,4):100,
(5900,1):100,(5900,2):100,(5900,3):100,(5900,4):100,
(6000,1):100,(6000,2):100,(6000,3):100,(6000,4):100,
(7200,1):77,(7200,2):86,(7200,3):91,(7200,4):94,
(7300,1):77,(7300,2):86,(7300,3):91,(7300,4):94,
(7500,1):77,(7500,2):86,(7500,3):91,(7500,4):94,
(7600,1):77,(7600,2):86,(7600,3):91,(7600,4):94}
curveNumber = np.ones(lu.shape, dtype = np.int)
for k,v in table.iteritems():
curveNumber[((lu == k[0]) & (sc == k[1]))] = v
cnArray = arcpy.NumPyArrayToRaster(curveNumber,pnt,cellSize,cellSize)
return cnArray
def mddelookupgood(inLU,inSoil):
lu = arcpy.RasterToNumPyArray(inLU)
sc = arcpy.RasterToNumPyArray(inSoil)
desc = arcpy.Describe(inLU)
cellSize = desc.meanCellHeight
extent = desc.Extent
pnt = arcpy.Point(extent.XMin,extent.YMin)
table = {(10,1):61,(10,2):75,(10,3):83,(10,4):87,
(11,1):54,(11,2):70,(11,3):80,(11,4):85,
(12,1):61,(12,2):75,(12,3):83,(12,4):87,
(13,1):77,(13,2):85,(13,3):90,(13,4):92,
(14,1):89,(14,2):92,(14,3):94,(14,4):95,
(15,1):81,(15,2):88,(15,3):91,(15,4):93,
(16,1):69,(16,2):80,(16,3):86,(16,4):89,
(17,1):77,(17,2):86,(17,3):91,(17,4):94,
(18,1):39,(18,2):61,(18,3):74,(18,4):80,
(20,1):67,(20,2):78,(20,3):85,(20,4):89,
(21,1):67,(21,2):78,(21,3):85,(21,4):89,
(22,1):39,(22,2):61,(22,3):74,(22,4):80,
(23,1):32,(23,2):58,(23,3):72,(23,4):79,
(24,1):59,(24,2):74,(24,3):82,(24,4):86,
(25,1):67,(25,2):78,(25,3):85,(25,4):89,
(40,1):30,(40,2):55,(40,3):70,(40,4):77,
(41,1):30,(41,2):55,(41,3):70,(41,4):77,
(42,1):30,(42,2):55,(42,3):70,(42,4):77,
(43,1):30,(43,2):55,(43,3):70,(43,4):77,
(44,1):30,(44,2):48,(44,3):65,(44,4):73,
(50,1):100,(50,2):100,(50,3):100,(50,4):100,
(60,1):100,(60,2):100,(60,3):100,(60,4):100,
(70,1):77,(70,2):86,(70,3):91,(70,4):94,
(71,1):77,(71,2):86,(71,3):91,(71,4):94,
(72,1):77,(72,2):86,(72,3):91,(72,4):94,
(73,1):77,(73,2):86,(73,3):91,(73,4):94,
(80,1):83,(80,2):89,(80,3):92,(80,4):94,
(191,1):67,(191,2):78,(191,3):85,(191,4):89,
(192,1):30,(192,2):55,(192,3):70,(192,4):77,
(241,1):59,(241,2):74,(241,3):82,(241,4):86,
(242,1):59,(242,2):74,(242,3):82,(242,4):86,
(1110,1):61,(1110,2):75,(1110,3):83,(1110,4):87,
(1120,1):77,(1120,2):85,(1120,3):90,(1120,4):92,
(1140,1):77,(1140,2):85,(1140,3):90,(1140,4):92,
(1200,1):89,(1200,2):92,(1200,3):94,(1200,4):95,
(1210,1):89,(1210,2):92,(1210,3):94,(1210,4):95,
(1220,1):89,(1220,2):92,(1220,3):94,(1220,4):95,
(1230,1):89,(1230,2):92,(1230,3):94,(1230,4):95,
(1250,1):89,(1250,2):92,(1250,3):94,(1250,4):95,
(1290,1):89,(1290,2):92,(1290,3):94,(1290,4):95,
(1300,1):81,(1300,2):88,(1300,3):91,(1300,4):93,
(1400,1):83,(1400,2):89,(1400,3):92,(1400,4):94,
(1410,1):83,(1410,2):89,(1410,3):92,(1410,4):94,
(1420,1):89,(1420,2):92,(1420,3):94,(1420,4):95,
(1430,1):83,(1430,2):89,(1430,3):92,(1430,4):94,
(1440,1):83,(1440,2):89,(1440,3):92,(1440,4):94,
(1450,1):39,(1450,2):61,(1450,3):74,(1450,4):80,
(1460,1):83,(1460,2):89,(1460,3):92,(1460,4):94,
(1490,1):83,(1490,2):89,(1490,3):92,(1490,4):94,
(1500,1):39,(1500,2):61,(1500,3):74,(1500,4):80,
(1600,1):69,(1600,2):80,(1600,3):86,(1600,4):89,
(1700,1):69,(1700,2):80,(1700,3):86,(1700,4):89,
(1800,1):69,(1800,2):80,(1800,3):86,(1800,4):89,
(1900,1):39,(1900,2):61,(1900,3):74,(1900,4):80,
(2100,1):67,(2100,2):78,(2100,3):85,(2100,4):89,
(2110,1):67,(2110,2):78,(2110,3):85,(2110,4):89,
(2120,1):39,(2120,2):61,(2120,3):74,(2120,4):80,
(2130,1):67,(2130,2):78,(2130,3):85,(2130,4):89,
(2140,1):67,(2140,2):78,(2140,3):85,(2140,4):89,
(2150,1):67,(2150,2):78,(2150,3):85,(2150,4):89,
(2200,1):32,(2200,2):58,(2200,3):72,(2200,4):79,
(2300,1):59,(2300,2):74,(2300,3):82,(2300,4):86,
(2400,1):59,(2400,2):74,(2400,3):82,(2400,4):86,
(2900,1):67,(2900,2):78,(2900,3):85,(2900,4):89,
(3100,1):30,(3100,2):55,(3100,3):70,(3100,4):77,
(3200,1):30,(3200,2):55,(3200,3):70,(3200,4):77,
(3300,1):30,(3300,2):55,(3300,3):70,(3300,4):77,
(4100,1):30,(4100,2):55,(4100,3):70,(4100,4):77,
(4200,1):30,(4200,2):55,(4200,3):70,(4200,4):77,
(4300,1):30,(4300,2):55,(4300,3):70,(4300,4):77,
(4400,1):30,(4400,2):55,(4400,3):70,(4400,4):77,
(5100,1):100,(5100,2):100,(5100,3):100,(5100,4):100,
(5200,1):100,(5200,2):100,(5200,3):100,(5200,4):100,
(5300,1):100,(5300,2):100,(5300,3):100,(5300,4):100,
(5400,1):100,(5400,2):100,(5400,3):100,(5400,4):100,
(5900,1):100,(5900,2):100,(5900,3):100,(5900,4):100,
(6000,1):100,(6000,2):100,(6000,3):100,(6000,4):100,
(7200,1):77,(7200,2):86,(7200,3):91,(7200,4):94,
(7300,1):77,(7300,2):86,(7300,3):91,(7300,4):94,
(7500,1):77,(7500,2):86,(7500,3):91,(7500,4):94,
(7600,1):77,(7600,2):86,(7600,3):91,(7600,4):94}
curveNumber = np.ones(lu.shape, dtype = np.int)
for k,v in table.iteritems():
curveNumber[((lu == k[0]) & (sc == k[1]))] = v
cnArray = arcpy.NumPyArrayToRaster(curveNumber,pnt,cellSize,cellSize)
return cnArray
def mddelookuppoor(inLU,inSoil):
lu = arcpy.RasterToNumPyArray(inLU)
sc = arcpy.RasterToNumPyArray(inSoil)
desc = arcpy.Describe(inLU)
cellSize = desc.meanCellHeight
extent = desc.Extent
pnt = arcpy.Point(extent.XMin,extent.YMin)
table = {(10,1):79,(10,2):86,(10,3):91,(10,4):92,
(11,1):76,(11,2):84,(11,3):89,(11,4):91,
(12,1):79,(12,2):86,(12,3):91,(12,4):92,
(13,1):88,(13,2):91,(13,3):94,(13,4):95,
(14,1):94,(14,2):95,(14,3):96,(14,4):97,
(15,1):90,(15,2):93,(15,3):95,(15,4):95,
(16,1):83,(16,2):89,(16,3):92,(16,4):94,
(17,1):77,(17,2):86,(17,3):91,(17,4):94,
(18,1):68,(18,2):79,(18,3):86,(18,4):89,
(20,1):72,(20,2):81,(20,3):88,(20,4):91,
(21,1):72,(21,2):81,(21,3):88,(21,4):91,
(22,1):68,(22,2):79,(22,3):86,(22,4):89,
(23,1):57,(23,2):73,(23,3):82,(23,4):86,
(24,1):59,(24,2):74,(24,3):82,(24,4):86,
(25,1):72,(25,2):81,(25,3):88,(25,4):91,
(40,1):45,(40,2):66,(40,3):77,(40,4):83,
(41,1):45,(41,2):66,(41,3):77,(41,4):83,
(42,1):45,(42,2):66,(42,3):77,(42,4):83,
(43,1):45,(43,2):66,(43,3):77,(43,4):83,
(44,1):48,(44,2):67,(44,3):77,(44,4):83,
(50,1):100,(50,2):100,(50,3):100,(50,4):100,
(60,1):100,(60,2):100,(60,3):100,(60,4):100,
(70,1):77,(70,2):86,(70,3):91,(70,4):94,
(71,1):77,(71,2):86,(71,3):91,(71,4):94,
(72,1):77,(72,2):86,(72,3):91,(72,4):94,
(73,1):77,(73,2):86,(73,3):91,(73,4):94,
(80,1):98,(80,2):98,(80,3):98,(80,4):98,
(191,1):72,(191,2):81,(191,3):88,(191,4):91,
(192,1):45,(192,2):66,(192,3):77,(192,4):83,
(241,1):59,(241,2):74,(241,3):82,(241,4):86,
(242,1):59,(242,2):74,(242,3):82,(242,4):86,
(1110,1):79,(1110,2):86,(1110,3):91,(1110,4):92,
(1120,1):88,(1120,2):91,(1120,3):94,(1120,4):95,
(1140,1):88,(1140,2):91,(1140,3):94,(1140,4):95,
(1200,1):94,(1200,2):95,(1200,3):96,(1200,4):97,
(1210,1):94,(1210,2):95,(1210,3):96,(1210,4):97,
(1220,1):94,(1220,2):95,(1220,3):96,(1220,4):97,
(1230,1):94,(1230,2):95,(1230,3):96,(1230,4):97,
(1250,1):94,(1250,2):95,(1250,3):96,(1250,4):97,
(1290,1):94,(1290,2):95,(1290,3):96,(1290,4):97,
(1300,1):90,(1300,2):93,(1300,3):95,(1300,4):95,
(1400,1):98,(1400,2):98,(1400,3):98,(1400,4):98,
(1410,1):98,(1410,2):98,(1410,3):98,(1410,4):98,
(1420,1):94,(1420,2):95,(1420,3):96,(1420,4):97,
(1430,1):98,(1430,2):98,(1430,3):98,(1430,4):98,
(1440,1):98,(1440,2):98,(1440,3):98,(1440,4):98,
(1450,1):68,(1450,2):79,(1450,3):86,(1450,4):89,
(1460,1):98,(1460,2):98,(1460,3):98,(1460,4):98,
(1490,1):98,(1490,2):98,(1490,3):98,(1490,4):98,
(1500,1):68,(1500,2):79,(1500,3):86,(1500,4):89,
(1600,1):83,(1600,2):89,(1600,3):92,(1600,4):94,
(1700,1):83,(1700,2):89,(1700,3):92,(1700,4):94,
(1800,1):83,(1800,2):89,(1800,3):92,(1800,4):94,
(1900,1):68,(1900,2):79,(1900,3):86,(1900,4):89,
(2100,1):72,(2100,2):81,(2100,3):88,(2100,4):91,
(2110,1):72,(2110,2):81,(2110,3):88,(2110,4):91,
(2120,1):68,(2120,2):79,(2120,3):86,(2120,4):89,
(2130,1):72,(2130,2):81,(2130,3):88,(2130,4):91,
(2140,1):72,(2140,2):81,(2140,3):88,(2140,4):91,
(2150,1):72,(2150,2):81,(2150,3):88,(2150,4):91,
(2200,1):57,(2200,2):73,(2200,3):82,(2200,4):86,
(2300,1):59,(2300,2):74,(2300,3):82,(2300,4):86,
(2400,1):59,(2400,2):74,(2400,3):82,(2400,4):86,
(2900,1):72,(2900,2):81,(2900,3):88,(2900,4):91,
(3100,1):45,(3100,2):66,(3100,3):77,(3100,4):83,
(3200,1):45,(3200,2):66,(3200,3):77,(3200,4):83,
(3300,1):45,(3300,2):66,(3300,3):77,(3300,4):83,
(4100,1):45,(4100,2):66,(4100,3):77,(4100,4):83,
(4200,1):45,(4200,2):66,(4200,3):77,(4200,4):83,
(4300,1):45,(4300,2):66,(4300,3):77,(4300,4):83,
(4400,1):45,(4400,2):66,(4400,3):77,(4400,4):83,
(5100,1):100,(5100,2):100,(5100,3):100,(5100,4):100,
(5200,1):100,(5200,2):100,(5200,3):100,(5200,4):100,
(5300,1):100,(5300,2):100,(5300,3):100,(5300,4):100,
(5400,1):100,(5400,2):100,(5400,3):100,(5400,4):100,
(5900,1):100,(5900,2):100,(5900,3):100,(5900,4):100,
(6000,1):100,(6000,2):100,(6000,3):100,(6000,4):100,
(7200,1):77,(7200,2):86,(7200,3):91,(7200,4):94,
(7300,1):77,(7300,2):86,(7300,3):91,(7300,4):94,
(7500,1):77,(7500,2):86,(7500,3):91,(7500,4):94,
(7600,1):77,(7600,2):86,(7600,3):91,(7600,4):94}
curveNumber = np.ones(lu.shape, dtype = np.int)
for k,v in table.iteritems():
curveNumber[((lu == k[0]) & (sc == k[1]))] = v
cnArray = arcpy.NumPyArrayToRaster(curveNumber,pnt,cellSize,cellSize)
return cnArray
#*******************************************************************************************************
# Calculate curve number for poor soil condition using RAGAN or SSURGO soil data using zoning tables
#*******************************************************************************************************
def zoninglookupfair(inLU,inSoil):
lu = arcpy.RasterToNumPyArray(inLU)
sc = arcpy.RasterToNumPyArray(inSoil)
desc = arcpy.Describe(inLU)
cellSize = desc.meanCellHeight
extent = desc.Extent
pnt = arcpy.Point(extent.XMin,extent.YMin)
table = {(10,1):68,(10,2):80,(10,3):86,(10,4):89,
(11,1):61,(11,2):76,(11,3):84,(11,4):88,
(12,1):68,(12,2):80,(12,3):86,(12,4):89,
(13,1):81,(13,2):88,(13,3):91,(13,4):93,
(14,1):91,(14,2):94,(14,3):95,(14,4):96,
(15,1):81,(15,2):88,(15,3):91,(15,4):93,
(16,1):74,(16,2):84,(16,3):89,(16,4):91,
(17,1):77,(17,2):86,(17,3):91,(17,4):94,
(18,1):49,(18,2):69,(18,3):79,(18,4):84,
(20,1):70,(20,2):80,(20,3):87,(20,4):90,
(21,1):70,(21,2):80,(21,3):87,(21,4):90,
(22,1):49,(22,2):69,(22,3):79,(22,4):84,
(23,1):43,(23,2):65,(23,3):76,(23,4):82,
(24,1):59,(24,2):74,(24,3):82,(24,4):86,
(25,1):70,(25,2):80,(25,3):87,(25,4):90,
(40,1):36,(40,2):60,(40,3):73,(40,4):79,
(41,1):36,(41,2):60,(41,3):73,(41,4):79,
(42,1):36,(42,2):60,(42,3):73,(42,4):79,
(43,1):36,(43,2):60,(43,3):73,(43,4):79,
(44,1):35,(44,2):56,(44,3):70,(44,4):77,
(50,1):100,(50,2):100,(50,3):100,(50,4):100,
(60,1):100,(60,2):100,(60,3):100,(60,4):100,
(70,1):77,(70,2):86,(70,3):91,(70,4):94,
(71,1):77,(71,2):86,(71,3):91,(71,4):94,
(72,1):77,(72,2):86,(72,3):91,(72,4):94,
(73,1):77,(73,2):86,(73,3):91,(73,4):94,
(80,1):98,(80,2):98,(80,3):98,(80,4):98,
(111,1):55,(111,2):72,(111,3):81,(111,4):86,
(112,1):59,(112,2):75,(112,3):83,(112,4):87,
(113,1):61,(113,2):76,(113,3):84,(113,4):88,
(114,1):64,(114,2):78,(114,3):85,(114,4):88,
(115,1):68,(115,2):80,(115,3):86,(115,4):89,
(116,1):81,(116,2):88,(116,3):91,(116,4):93,
(191,1):70,(191,2):80,(191,3):87,(191,4):90,
(192,1):36,(192,2):60,(192,3):73,(192,4):79,
(241,1):59,(241,2):74,(241,3):82,(241,4):86,
(242,1):59,(242,2):74,(242,3):82,(242,4):86,
(1110,1):68,(1110,2):80,(1110,3):86,(1110,4):89,
(1120,1):81,(1120,2):88,(1120,3):91,(1120,4):93,
(1140,1):81,(1140,2):88,(1140,3):91,(1140,4):93,
(1200,1):91,(1200,2):94,(1200,3):95,(1200,4):96,
(1210,1):91,(1210,2):94,(1210,3):95,(1210,4):96,
(1220,1):91,(1220,2):94,(1220,3):95,(1220,4):96,
(1230,1):91,(1230,2):94,(1230,3):95,(1230,4):96,
(1250,1):91,(1250,2):94,(1250,3):95,(1250,4):96,
(1290,1):91,(1290,2):94,(1290,3):95,(1290,4):96,
(1300,1):81,(1300,2):88,(1300,3):91,(1300,4):93,
(1400,1):98,(1400,2):98,(1400,3):98,(1400,4):98,
(1410,1):98,(1410,2):98,(1410,3):98,(1410,4):98,
(1420,1):91,(1420,2):94,(1420,3):95,(1420,4):96,
(1430,1):98,(1430,2):98,(1430,3):98,(1430,4):98,
(1440,1):98,(1440,2):98,(1440,3):98,(1440,4):98,
(1450,1):49,(1450,2):69,(1450,3):79,(1450,4):84,
(1460,1):98,(1460,2):98,(1460,3):98,(1460,4):98,
(1490,1):98,(1490,2):98,(1490,3):98,(1490,4):98,
(1500,1):49,(1500,2):69,(1500,3):79,(1500,4):84,
(1600,1):74,(1600,2):84,(1600,3):89,(1600,4):91,
(1700,1):74,(1700,2):84,(1700,3):89,(1700,4):91,
(1800,1):74,(1800,2):84,(1800,3):89,(1800,4):91,
(1900,1):49,(1900,2):69,(1900,3):79,(1900,4):84,
(2100,1):70,(2100,2):80,(2100,3):87,(2100,4):90,
(2110,1):70,(2110,2):80,(2110,3):87,(2110,4):90,
(2120,1):49,(2120,2):69,(2120,3):79,(2120,4):84,
(2130,1):70,(2130,2):80,(2130,3):87,(2130,4):90,
(2140,1):70,(2140,2):80,(2140,3):87,(2140,4):90,
(2150,1):70,(2150,2):80,(2150,3):87,(2150,4):90,
(2200,1):43,(2200,2):65,(2200,3):76,(2200,4):82,
(2300,1):59,(2300,2):74,(2300,3):82,(2300,4):86,
(2400,1):59,(2400,2):74,(2400,3):82,(2400,4):86,
(2900,1):70,(2900,2):80,(2900,3):87,(2900,4):90,
(3100,1):36,(3100,2):60,(3100,3):73,(3100,4):79,
(3200,1):36,(3200,2):60,(3200,3):73,(3200,4):79,
(3300,1):36,(3300,2):60,(3300,3):73,(3300,4):79,
(4100,1):36,(4100,2):60,(4100,3):73,(4100,4):79,
(4200,1):36,(4200,2):60,(4200,3):73,(4200,4):79,
(4300,1):36,(4300,2):60,(4300,3):73,(4300,4):79,
(4400,1):36,(4400,2):60,(4400,3):73,(4400,4):79,
(5100,1):100,(5100,2):100,(5100,3):100,(5100,4):100,
(5200,1):100,(5200,2):100,(5200,3):100,(5200,4):100,
(5300,1):100,(5300,2):100,(5300,3):100,(5300,4):100,
(5400,1):100,(5400,2):100,(5400,3):100,(5400,4):100,
(5900,1):100,(5900,2):100,(5900,3):100,(5900,4):100,
(6000,1):100,(6000,2):100,(6000,3):100,(6000,4):100,
(7200,1):77,(7200,2):86,(7200,3):91,(7200,4):94,
(7300,1):77,(7300,2):86,(7300,3):91,(7300,4):94,
(7500,1):77,(7500,2):86,(7500,3):91,(7500,4):94,
(7600,1):77,(7600,2):86,(7600,3):91,(7600,4):94}
curveNumber = np.ones(lu.shape, dtype = np.int)
for k,v in table.iteritems():
curveNumber[((lu == k[0]) & (sc == k[1]))] = v
cnArray = arcpy.NumPyArrayToRaster(curveNumber,pnt,cellSize,cellSize)
return cnArray
def zoninglookupgood(inLU,inSoil):
lu = arcpy.RasterToNumPyArray(inLU)
sc = arcpy.RasterToNumPyArray(inSoil)
desc = arcpy.Describe(inLU)
cellSize = desc.meanCellHeight
extent = desc.Extent
pnt = arcpy.Point(extent.XMin,extent.YMin)
table = {(10,1):61,(10,2):75,(10,3):83,(10,4):87,
(11,1):54,(11,2):70,(11,3):80,(11,4):85,
(12,1):61,(12,2):75,(12,3):83,(12,4):87,
(13,1):77,(13,2):85,(13,3):90,(13,4):92,
(14,1):89,(14,2):92,(14,3):94,(14,4):95,
(15,1):81,(15,2):88,(15,3):91,(15,4):93,
(16,1):69,(16,2):80,(16,3):86,(16,4):89,
(17,1):77,(17,2):86,(17,3):91,(17,4):94,
(18,1):39,(18,2):61,(18,3):74,(18,4):80,
(20,1):67,(20,2):78,(20,3):85,(20,4):89,
(21,1):67,(21,2):78,(21,3):85,(21,4):89,
(22,1):39,(22,2):61,(22,3):74,(22,4):80,
(23,1):32,(23,2):58,(23,3):72,(23,4):79,
(24,1):59,(24,2):74,(24,3):82,(24,4):86,
(25,1):67,(25,2):78,(25,3):85,(25,4):89,
(40,1):30,(40,2):55,(40,3):70,(40,4):77,
(41,1):30,(41,2):55,(41,3):70,(41,4):77,
(42,1):30,(42,2):55,(42,3):70,(42,4):77,
(43,1):30,(43,2):55,(43,3):70,(43,4):77,
(44,1):30,(44,2):48,(44,3):65,(44,4):73,
(50,1):100,(50,2):100,(50,3):100,(50,4):100,
(60,1):100,(60,2):100,(60,3):100,(60,4):100,
(70,1):77,(70,2):86,(70,3):91,(70,4):94,
(71,1):77,(71,2):86,(71,3):91,(71,4):94,
(72,1):77,(72,2):86,(72,3):91,(72,4):94,
(73,1):77,(73,2):86,(73,3):91,(73,4):94,
(80,1):98,(80,2):98,(80,3):98,(80,4):98,
(111,1):46,(111,2):65,(111,3):77,(111,4):82,
(112,1):51,(112,2):68,(112,3):79,(112,4):84,
(113,1):54,(113,2):70,(113,3):80,(113,4):85,
(114,1):57,(114,2):72,(114,3):81,(114,4):86,
(115,1):61,(115,2):75,(115,3):83,(115,4):87,
(116,1):77,(116,2):85,(116,3):90,(116,4):92,
(191,1):67,(191,2):78,(191,3):85,(191,4):89,
(192,1):30,(192,2):55,(192,3):70,(192,4):77,
(241,1):59,(241,2):74,(241,3):82,(241,4):86,
(242,1):59,(242,2):74,(242,3):82,(242,4):86,
(1110,1):61,(1110,2):75,(1110,3):83,(1110,4):87,
(1120,1):77,(1120,2):85,(1120,3):90,(1120,4):92,
(1140,1):77,(1140,2):85,(1140,3):90,(1140,4):92,
(1200,1):89,(1200,2):92,(1200,3):94,(1200,4):95,
(1210,1):89,(1210,2):92,(1210,3):94,(1210,4):95,
(1220,1):89,(1220,2):92,(1220,3):94,(1220,4):95,
(1230,1):89,(1230,2):92,(1230,3):94,(1230,4):95,
(1250,1):89,(1250,2):92,(1250,3):94,(1250,4):95,
(1290,1):89,(1290,2):92,(1290,3):94,(1290,4):95,
(1300,1):81,(1300,2):88,(1300,3):91,(1300,4):93,
(1400,1):83,(1400,2):89,(1400,3):92,(1400,4):94,
(1410,1):83,(1410,2):89,(1410,3):92,(1410,4):94,
(1420,1):89,(1420,2):92,(1420,3):94,(1420,4):95,
(1430,1):83,(1430,2):89,(1430,3):92,(1430,4):94,
(1440,1):83,(1440,2):89,(1440,3):92,(1440,4):94,
(1450,1):39,(1450,2):61,(1450,3):74,(1450,4):80,
(1460,1):83,(1460,2):89,(1460,3):92,(1460,4):94,
(1490,1):83,(1490,2):89,(1490,3):92,(1490,4):94,
(1500,1):39,(1500,2):61,(1500,3):74,(1500,4):80,
(1600,1):69,(1600,2):80,(1600,3):86,(1600,4):89,
(1700,1):69,(1700,2):80,(1700,3):86,(1700,4):89,
(1800,1):69,(1800,2):80,(1800,3):86,(1800,4):89,
(1900,1):39,(1900,2):61,(1900,3):74,(1900,4):80,
(2100,1):67,(2100,2):78,(2100,3):85,(2100,4):89,
(2110,1):67,(2110,2):78,(2110,3):85,(2110,4):89,
(2120,1):39,(2120,2):61,(2120,3):74,(2120,4):80,
(2130,1):67,(2130,2):78,(2130,3):85,(2130,4):89,
(2140,1):67,(2140,2):78,(2140,3):85,(2140,4):89,
(2150,1):67,(2150,2):78,(2150,3):85,(2150,4):89,
(2200,1):32,(2200,2):58,(2200,3):72,(2200,4):79,
(2300,1):59,(2300,2):74,(2300,3):82,(2300,4):86,
(2400,1):59,(2400,2):74,(2400,3):82,(2400,4):86,
(2900,1):67,(2900,2):78,(2900,3):85,(2900,4):89,
(3100,1):30,(3100,2):55,(3100,3):70,(3100,4):77,
(3200,1):30,(3200,2):55,(3200,3):70,(3200,4):77,
(3300,1):30,(3300,2):55,(3300,3):70,(3300,4):77,
(4100,1):30,(4100,2):55,(4100,3):70,(4100,4):77,
(4200,1):30,(4200,2):55,(4200,3):70,(4200,4):77,
(4300,1):30,(4300,2):55,(4300,3):70,(4300,4):77,
(4400,1):30,(4400,2):55,(4400,3):70,(4400,4):77,
(5100,1):100,(5100,2):100,(5100,3):100,(5100,4):100,
(5200,1):100,(5200,2):100,(5200,3):100,(5200,4):100,
(5300,1):100,(5300,2):100,(5300,3):100,(5300,4):100,
(5400,1):100,(5400,2):100,(5400,3):100,(5400,4):100,
(5900,1):100,(5900,2):100,(5900,3):100,(5900,4):100,
(6000,1):100,(6000,2):100,(6000,3):100,(6000,4):100,
(7200,1):77,(7200,2):86,(7200,3):91,(7200,4):94,
(7300,1):77,(7300,2):86,(7300,3):91,(7300,4):94,
(7500,1):77,(7500,2):86,(7500,3):91,(7500,4):94,
(7600,1):77,(7600,2):86,(7600,3):91,(7600,4):94}
curveNumber = np.ones(lu.shape, dtype = np.int)
for k,v in table.iteritems():
curveNumber[((lu == k[0]) & (sc == k[1]))] = v
cnArray = arcpy.NumPyArrayToRaster(curveNumber,pnt,cellSize,cellSize)
return cnArray
def zoninglookuppoor(inLU,inSoil):
lu = arcpy.RasterToNumPyArray(inLU)
sc = arcpy.RasterToNumPyArray(inSoil)
desc = arcpy.Describe(inLU)
cellSize = desc.meanCellHeight
extent = desc.Extent
pnt = arcpy.Point(extent.XMin,extent.YMin)
table = {(10,1):79,(10,2):86,(10,3):91,(10,4):92,
(11,1):76,(11,2):84,(11,3):89,(11,4):91,
(12,1):79,(12,2):86,(12,3):91,(12,4):92,
(13,1):88,(13,2):91,(13,3):94,(13,4):95,
(14,1):94,(14,2):95,(14,3):96,(14,4):97,
(15,1):90,(15,2):93,(15,3):95,(15,4):95,
(16,1):83,(16,2):89,(16,3):92,(16,4):94,
(17,1):77,(17,2):86,(17,3):91,(17,4):94,
(18,1):68,(18,2):79,(18,3):86,(18,4):89,
(20,1):72,(20,2):81,(20,3):88,(20,4):91,
(21,1):72,(21,2):81,(21,3):88,(21,4):91,
(22,1):68,(22,2):79,(22,3):86,(22,4):89,
(23,1):57,(23,2):73,(23,3):82,(23,4):86,
(24,1):59,(24,2):74,(24,3):82,(24,4):86,
(25,1):72,(25,2):81,(25,3):88,(25,4):91,
(40,1):45,(40,2):66,(40,3):77,(40,4):83,
(41,1):45,(41,2):66,(41,3):77,(41,4):83,
(42,1):45,(42,2):66,(42,3):77,(42,4):83,
(43,1):45,(43,2):66,(43,3):77,(43,4):83,
(44,1):48,(44,2):67,(44,3):77,(44,4):83,
(50,1):100,(50,2):100,(50,3):100,(50,4):100,
(60,1):100,(60,2):100,(60,3):100,(60,4):100,
(70,1):77,(70,2):86,(70,3):91,(70,4):94,
(71,1):77,(71,2):86,(71,3):91,(71,4):94,
(72,1):77,(72,2):86,(72,3):91,(72,4):94,
(73,1):77,(73,2):86,(73,3):91,(73,4):94,
(80,1):98,(80,2):98,(80,3):98,(80,4):98,
(111,1):72,(111,2):81,(111,3):87,(111,4):90,
(112,1):74,(112,2):83,(112,3):88,(112,4):91,
(113,1):76,(113,2):84,(113,3):89,(113,4):91,
(114,1):77,(114,2):85,(114,3):90,(114,4):92,
(115,1):79,(115,2):86,(115,3):91,(115,4):92,
(116,1):88,(116,2):91,(116,3):94,(116,4):95,
(191,1):72,(191,2):81,(191,3):88,(191,4):91,
(192,1):45,(192,2):66,(192,3):77,(192,4):83,
(241,1):59,(241,2):74,(241,3):82,(241,4):86,
(242,1):59,(242,2):74,(242,3):82,(242,4):86,
(1110,1):79,(1110,2):86,(1110,3):91,(1110,4):92,
(1120,1):88,(1120,2):91,(1120,3):94,(1120,4):95,
(1140,1):88,(1140,2):91,(1140,3):94,(1140,4):95,
(1200,1):94,(1200,2):95,(1200,3):96,(1200,4):97,
(1210,1):94,(1210,2):95,(1210,3):96,(1210,4):97,
(1220,1):94,(1220,2):95,(1220,3):96,(1220,4):97,
(1230,1):94,(1230,2):95,(1230,3):96,(1230,4):97,
(1250,1):94,(1250,2):95,(1250,3):96,(1250,4):97,
(1290,1):94,(1290,2):95,(1290,3):96,(1290,4):97,
(1300,1):90,(1300,2):93,(1300,3):95,(1300,4):95,
(1400,1):98,(1400,2):98,(1400,3):98,(1400,4):98,
(1410,1):98,(1410,2):98,(1410,3):98,(1410,4):98,
(1420,1):94,(1420,2):95,(1420,3):96,(1420,4):97,
(1430,1):98,(1430,2):98,(1430,3):98,(1430,4):98,
(1440,1):98,(1440,2):98,(1440,3):98,(1440,4):98,
(1450,1):68,(1450,2):79,(1450,3):86,(1450,4):89,
(1460,1):98,(1460,2):98,(1460,3):98,(1460,4):98,
(1490,1):98,(1490,2):98,(1490,3):98,(1490,4):98,
(1500,1):68,(1500,2):79,(1500,3):86,(1500,4):89,
(1600,1):83,(1600,2):89,(1600,3):92,(1600,4):94,
(1700,1):83,(1700,2):89,(1700,3):92,(1700,4):94,
(1800,1):83,(1800,2):89,(1800,3):92,(1800,4):94,
(1900,1):68,(1900,2):79,(1900,3):86,(1900,4):89,
(2100,1):72,(2100,2):81,(2100,3):88,(2100,4):91,
(2110,1):72,(2110,2):81,(2110,3):88,(2110,4):91,
(2120,1):68,(2120,2):79,(2120,3):86,(2120,4):89,
(2130,1):72,(2130,2):81,(2130,3):88,(2130,4):91,
(2140,1):72,(2140,2):81,(2140,3):88,(2140,4):91,
(2150,1):72,(2150,2):81,(2150,3):88,(2150,4):91,
(2200,1):57,(2200,2):73,(2200,3):82,(2200,4):86,
(2300,1):59,(2300,2):74,(2300,3):82,(2300,4):86,
(2400,1):59,(2400,2):74,(2400,3):82,(2400,4):86,
(2900,1):72,(2900,2):81,(2900,3):88,(2900,4):91,
(3100,1):45,(3100,2):66,(3100,3):77,(3100,4):83,
(3200,1):45,(3200,2):66,(3200,3):77,(3200,4):83,
(3300,1):45,(3300,2):66,(3300,3):77,(3300,4):83,
(4100,1):45,(4100,2):66,(4100,3):77,(4100,4):83,
(4200,1):45,(4200,2):66,(4200,3):77,(4200,4):83,
(4300,1):45,(4300,2):66,(4300,3):77,(4300,4):83,
(4400,1):45,(4400,2):66,(4400,3):77,(4400,4):83,
(5100,1):100,(5100,2):100,(5100,3):100,(5100,4):100,
(5200,1):100,(5200,2):100,(5200,3):100,(5200,4):100,
(5300,1):100,(5300,2):100,(5300,3):100,(5300,4):100,
(5400,1):100,(5400,2):100,(5400,3):100,(5400,4):100,
(5900,1):100,(5900,2):100,(5900,3):100,(5900,4):100,
(6000,1):100,(6000,2):100,(6000,3):100,(6000,4):100,
(7200,1):77,(7200,2):86,(7200,3):91,(7200,4):94,
(7300,1):77,(7300,2):86,(7300,3):91,(7300,4):94,
(7500,1):77,(7500,2):86,(7500,3):91,(7500,4):94,
(7600,1):77,(7600,2):86,(7600,3):91,(7600,4):94}
curveNumber = np.ones(lu.shape, dtype = np.int)
for k,v in table.iteritems():
curveNumber[((lu == k[0]) & (sc == k[1]))] = v
cnArray = arcpy.NumPyArrayToRaster(curveNumber,pnt,cellSize,cellSize)
return cnArray
#*******************************************************************************************************
# Calculate curve number for fair soil condition using RAGAN or SSURGO soil data using MRLC tables
#*******************************************************************************************************
def mrlclookupfair(inLU,inSoil):
lu = arcpy.RasterToNumPyArray(inLU)
sc = arcpy.RasterToNumPyArray(inSoil)
desc = arcpy.Describe(inLU)
cellSize = desc.meanCellHeight
extent = desc.Extent
pnt = arcpy.Point(extent.XMin,extent.YMin)
table = {(11,1):100,(11,2):100,(11,3):100,(11,4):100,
(21,1):61,(21,2):76,(21,3):84,(21,4):88,
(22,1):81,(22,2):88,(22,3):91,(22,4):93,
(23,1):91,(23,2):94,(23,3):95,(23,4):96,
(24,1):91,(24,2):94,(24,3):95,(24,4):96,
(31,1):77,(31,2):86,(31,3):91,(31,4):94,
(32,1):77,(32,2):86,(32,3):91,(32,4):94,
(33,1):77,(33,2):86,(33,3):91,(33,4):94,
(41,1):36,(41,2):60,(41,3):73,(41,4):79,
(42,1):36,(42,2):60,(42,3):73,(42,4):79,
(43,1):36,(43,2):60,(43,3):73,(43,4):79,
(81,1):49,(81,2):69,(81,3):79,(81,4):84,
(82,1):70,(82,2):80,(82,3):87,(82,4):90,
(85,1):49,(85,2):69,(85,3):79,(85,4):84,
(90,1):100,(90,2):100,(90,3):100,(90,4):100,
(91,1):100,(91,2):100,(91,3):100,(91,4):100,
(92,1):100,(92,2):100,(92,3):100,(92,4):100,
(95,1):100,(95,2):100,(95,3):100,(95,4):100}
curveNumber = np.ones(lu.shape, dtype = np.int)
for k,v in table.iteritems():
curveNumber[((lu == k[0]) & (sc == k[1]))] = v
cnArray = arcpy.NumPyArrayToRaster(curveNumber,pnt,cellSize,cellSize)
return cnArray
def mrlclookupgood(inLU,inSoil):
lu = arcpy.RasterToNumPyArray(inLU)
sc = arcpy.RasterToNumPyArray(inSoil)
desc = arcpy.Describe(inLU)
cellSize = desc.meanCellHeight
extent = desc.Extent
pnt = arcpy.Point(extent.XMin,extent.YMin)
table = {(11,1):100,(11,2):100,(11,3):100,(11,4):100,
(21,1):54,(21,2):70,(21,3):80,(21,4):85,
(22,1):77,(22,2):85,(22,3):90,(22,4):92,
(23,1):89,(23,2):92,(23,3):94,(23,4):95,
(24,1):89,(24,2):92,(24,3):94,(24,4):95,
(31,1):77,(31,2):86,(31,3):91,(31,4):94,
(32,1):77,(32,2):86,(32,3):91,(32,4):94,
(33,1):77,(33,2):86,(33,3):91,(33,4):94,
(41,1):30,(41,2):55,(41,3):70,(41,4):77,
(42,1):30,(42,2):55,(42,3):70,(42,4):77,
(43,1):30,(43,2):55,(43,3):70,(43,4):77,
(81,1):39,(81,2):61,(81,3):74,(81,4):80,
(82,1):67,(82,2):78,(82,3):85,(82,4):89,
(85,1):39,(85,2):61,(85,3):74,(85,4):80,
(90,1):100,(90,2):100,(90,3):100,(90,4):100,
(91,1):100,(91,2):100,(91,3):100,(91,4):100,
(92,1):100,(92,2):100,(92,3):100,(92,4):100,
(95,1):100,(95,2):100,(95,3):100,(95,4):100}
curveNumber = np.ones(lu.shape, dtype = np.int)
for k,v in table.iteritems():
curveNumber[((lu == k[0]) & (sc == k[1]))] = v
cnArray = arcpy.NumPyArrayToRaster(curveNumber,pnt,cellSize,cellSize)
return cnArray
def mrlclookuppoor(inLU,inSoil):
lu = arcpy.RasterToNumPyArray(inLU)
sc = arcpy.RasterToNumPyArray(inSoil)
desc = arcpy.Describe(inLU)
cellSize = desc.meanCellHeight
extent = desc.Extent
pnt = arcpy.Point(extent.XMin,extent.YMin)
table = {(11,1):100,(11,2):100,(11,3):100,(11,4):100,
(21,1):76,(21,2):84,(21,3):89,(21,4):91,
(22,1):88,(22,2):91,(22,3):94,(22,4):95,
(23,1):89,(23,2):92,(23,3):94,(23,4):95,
(24,1):94,(24,2):95,(24,3):96,(24,4):97,
(31,1):77,(31,2):86,(31,3):91,(31,4):94,
(32,1):77,(32,2):86,(32,3):91,(32,4):94,
(33,1):77,(33,2):86,(33,3):91,(33,4):94,
(41,1):45,(41,2):66,(41,3):77,(41,4):83,
(42,1):45,(42,2):66,(42,3):77,(42,4):83,
(43,1):45,(43,2):66,(43,3):77,(43,4):83,
(81,1):68,(81,2):79,(81,3):86,(81,4):89,
(82,1):72,(82,2):81,(82,3):88,(82,4):91,
(85,1):68,(85,2):79,(85,3):86,(85,4):89,
(90,1):100,(90,2):100,(90,3):100,(90,4):100,
(91,1):100,(91,2):100,(91,3):100,(91,4):100,
(92,1):100,(92,2):100,(92,3):100,(92,4):100,
(95,1):100,(95,2):100,(95,3):100,(95,4):100}
curveNumber = np.ones(lu.shape, dtype = np.int)
for k,v in table.iteritems():
curveNumber[((lu == k[0]) & (sc == k[1]))] = v
cnArray = arcpy.NumPyArrayToRaster(curveNumber,pnt,cellSize,cellSize)
return cnArray
#*******************************************************************************************************
# Calculate curve number for poor soil condition using RAGAN or SSURGO soil data using USGS tables
#*******************************************************************************************************
def usgslookupfair(inLU,inSoil):
lu = arcpy.RasterToNumPyArray(inLU)
sc = arcpy.RasterToNumPyArray(inSoil)
desc = arcpy.Describe(inLU)
cellSize = desc.meanCellHeight
extent = desc.Extent
pnt = arcpy.Point(extent.XMin,extent.YMin)
table = {(11,1):68,(11,2):80,(11,3):86,(11,4):89,
(12,1):91,(12,2):94,(12,3):95,(12,4):96,
(13,1):81,(13,2):88,(13,3):91,(13,4):93,
(14,1):98,(14,2):98,(14,3):98,(14,4):98,
(15,1):91,(15,2):94,(15,3):95,(15,4):96,
(16,1):68,(16,2):80,(16,3):86,(16,4):89,
(17,1):54,(17,2):72,(17,3):81,(17,4):86,
(21,1):70,(21,2):80,(21,3):87,(21,4):90,
(22,1):43,(22,2):65,(22,3):76,(22,4):82,
(23,1):59,(23,2):74,(23,3):82,(23,4):86,
(24,1):70,(24,2):80,(24,3):87,(24,4):90,
(41,1):36,(41,2):60,(41,3):73,(41,4):79,
(42,1):36,(42,2):60,(42,3):73,(42,4):79,
(43,1):36,(43,2):60,(43,3):73,(43,4):79,
(51,1):100,(51,2):100,(51,3):100,(51,4):100,
(52,1):100,(52,2):100,(52,3):100,(52,4):100,
(53,1):100,(53,2):100,(53,3):100,(53,4):100,
(54,1):100,(54,2):100,(54,3):100,(54,4):100,
(61,1):100,(61,2):100,(61,3):100,(61,4):100,
(62,1):100,(62,2):100,(62,3):100,(62,4):100,
(70,1):77,(70,2):86,(70,3):91,(70,4):94,
(71,1):100,(71,2):100,(71,3):100,(71,4):100,
(72,1):77,(72,2):86,(72,3):91,(72,4):94,
(73,1):77,(73,2):86,(73,3):91,(73,4):94,
(74,1):100,(74,2):100,(74,3):100,(74,4):100,
(75,1):77,(75,2):86,(75,3):91,(75,4):94,
(76,1):77,(76,2):86,(76,3):91,(76,4):94,
(77,1):77,(77,2):86,(77,3):91,(76,4):94}
curveNumber = np.ones(lu.shape, dtype = np.int)
for k,v in table.iteritems():
curveNumber[((lu == k[0]) & (sc == k[1]))] = v
cnArray = arcpy.NumPyArrayToRaster(curveNumber,pnt,cellSize,cellSize)
return cnArray
def usgslookupgood(inLU,inSoil):
lu = arcpy.RasterToNumPyArray(inLU)
sc = arcpy.RasterToNumPyArray(inSoil)
desc = arcpy.Describe(inLU)
cellSize = desc.meanCellHeight
extent = desc.Extent
pnt = arcpy.Point(extent.XMin,extent.YMin)
table = {(11,1):61,(11,2):75,(11,3):83,(11,4):87,
(12,1):89,(12,2):92,(12,3):94,(12,4):95,
(13,1):81,(13,2):88,(13,3):91,(13,4):93,
(14,1):98,(14,2):98,(14,3):98,(14,4):98,
(15,1):89,(15,2):92,(15,3):94,(15,4):95,
(16,1):61,(16,2):75,(16,3):83,(16,4):87,
(17,1):45,(17,2):65,(17,3):77,(17,4):82,
(21,1):65,(21,2):78,(21,3):85,(21,4):89,
(22,1):32,(22,2):58,(22,3):72,(22,4):79,
(23,1):59,(23,2):74,(23,3):82,(23,4):86,
(24,1):67,(24,2):78,(24,3):85,(24,4):89,
(41,1):30,(41,2):55,(41,3):70,(41,4):77,
(42,1):30,(42,2):55,(42,3):70,(42,4):77,
(43,1):30,(43,2):55,(43,3):70,(43,4):77,
(51,1):100,(51,2):100,(51,3):100,(51,4):100,
(52,1):100,(52,2):100,(52,3):100,(52,4):100,
(53,1):100,(53,2):100,(53,3):100,(53,4):100,
(54,1):100,(54,2):100,(54,3):100,(54,4):100,
(61,1):100,(61,2):100,(61,3):100,(61,4):100,
(62,1):100,(62,2):100,(62,3):100,(62,4):100,
(70,1):77,(70,2):86,(70,3):91,(70,4):94,
(71,1):100,(71,2):100,(71,3):100,(71,4):100,
(72,1):77,(72,2):86,(72,3):91,(72,4):94,
(73,1):77,(73,2):86,(73,3):91,(73,4):94,
(74,1):77,(74,2):86,(74,3):91,(74,4):94,
(75,1):77,(75,2):86,(75,3):91,(75,4):94,
(76,1):77,(76,2):86,(76,3):91,(76,4):94,
(77,1):77,(77,2):86,(77,3):91,(76,4):94}
curveNumber = np.ones(lu.shape, dtype = np.int)
for k,v in table.iteritems():
curveNumber[((lu == k[0]) & (sc == k[1]))] = v
cnArray = arcpy.NumPyArrayToRaster(curveNumber,pnt,cellSize,cellSize)
return cnArray
def usgslookuppoor(inLU,inSoil):
lu = arcpy.RasterToNumPyArray(inLU)
sc = arcpy.RasterToNumPyArray(inSoil)
desc = arcpy.Describe(inLU)
cellSize = desc.meanCellHeight
extent = desc.Extent
pnt = arcpy.Point(extent.XMin,extent.YMin)
table = {(11,1):79,(11,2):86,(11,3):91,(11,4):92,
(12,1):94,(12,2):95,(12,3):96,(12,4):97,
(13,1):90,(13,2):93,(13,3):95,(13,4):95,
(14,1):98,(14,2):98,(14,3):98,(14,4):98,
(15,1):94,(15,2):95,(15,3):96,(15,4):97,
(16,1):79,(16,2):86,(16,3):91,(16,4):92,
(17,1):71,(17,2):81,(17,3):87,(17,4):90,
(21,1):72,(21,2):81,(21,3):88,(21,4):91,
(22,1):57,(22,2):73,(22,3):82,(22,4):86,
(23,1):59,(23,2):74,(23,3):82,(23,4):86,
(24,1):72,(24,2):81,(24,3):88,(24,4):91,
(41,1):45,(41,2):66,(41,3):77,(41,4):83,
(42,1):45,(42,2):66,(42,3):77,(42,4):83,
(43,1):45,(43,2):66,(43,3):77,(43,4):83,
(51,1):100,(51,2):100,(51,3):100,(51,4):100,
(52,1):100,(52,2):100,(52,3):100,(52,4):100,
(53,1):100,(53,2):100,(53,3):100,(53,4):100,
(54,1):100,(54,2):100,(54,3):100,(54,4):100,
(61,1):100,(61,2):100,(61,3):100,(61,4):100,
(62,1):100,(62,2):100,(62,3):100,(62,4):100,
(70,1):77,(70,2):86,(70,3):91,(70,4):94,
(71,1):100,(71,2):100,(71,3):100,(71,4):100,
(72,1):77,(72,2):86,(72,3):91,(72,4):94,
(73,1):77,(73,2):86,(73,3):91,(73,4):94,
(74,1):100,(74,2):100,(74,3):100,(74,4):100,
(75,1):77,(75,2):86,(75,3):91,(75,4):94,
(76,1):77,(76,2):86,(76,3):91,(76,4):94,
(77,1):77,(77,2):86,(77,3):91,(76,4):94}
curveNumber = np.ones(lu.shape, dtype = np.int)
for k,v in table.iteritems():
curveNumber[((lu == k[0]) & (sc == k[1]))] = v
cnArray = arcpy.NumPyArrayToRaster(curveNumber,pnt,cellSize,cellSize)
return cnArray
#*******************************************************************************************************
# Calculate curve number for fair, good, and poor soil conditions using STATSGO soil data
# for NLCD Look Up table
#*******************************************************************************************************
def statsgonlcdlookupfair(inLU,inSoil,opath):
arcpy.env.workspace = opath
arcpy.env.scratchWorkspace = opath + "/aux_folder"
# read statsgo shapefile
arcpy.FeatureToRaster_conversion(inSoil, "PCT_A", opath + "/soilG_a", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_B", opath + "/soilG_b", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_C", opath + "/soilG_c", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_D", opath + "/soilG_d", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_WATER", opath + "/water_grid", 30)
# read landuse data
lusG_a = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,100],[21,49],[22,61],[23,68],[24,81],
[31,77],[41,36],[42,36],[43,36],[52,35],
[71,49],[81,49],[82,70],[90,100],
[95,100]]),"DATA")
lusG_b = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,100],[21,69],[22,76],[23,80],[24,88],
[31,86],[41,60],[42,60],[43,60],[52,56],
[71,69],[81,69],[82,80],[90,100],
[95,100]]),"DATA")
lusG_c = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,100],[21,79],[22,84],[23,86],[24,91],
[31,91],[41,73],[42,73],[43,73],[52,70],
[71,79],[81,79],[82,87],[90,100],
[95,100]]),"DATA")
lusG_d = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,100],[21,84],[22,88],[23,89],[24,93],
[31,94],[41,79],[42,79],[43,79],[52,77],
[71,84],[81,84],[82,90],[90,100],
[95,100]]),"DATA")
# Integer grid
numer = ((arcpy.sa.Times(opath + "/soilG_a",lusG_a))
+(arcpy.sa.Times(opath + "/soilG_b",lusG_b))
+(arcpy.sa.Times(opath + "/soilG_c",lusG_c))
+(arcpy.sa.Times(opath + "/soilG_d",lusG_d))
+(arcpy.sa.Times(opath + "/water_grid",100)))
denom = (arcpy.sa.Plus(opath + "/soilG_a",opath + "/soilG_b")
+(arcpy.sa.Plus(opath + "/soilG_c",opath + "/soilG_d"))
+(opath + "/water_grid"))
cnGrid = numer / denom
return cnGrid
def statsgonlcdlookupgood(inLU,inSoil,opath):
arcpy.env.workspace = opath
arcpy.env.scratchWorkspace = opath + "/aux_folder"
# read statsgo shapefile
arcpy.FeatureToRaster_conversion(inSoil, "PCT_A", opath + "/soilG_a", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_B", opath + "/soilG_b", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_C", opath + "/soilG_c", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_D", opath + "/soilG_d", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_WATER", opath + "/water_grid", 30)
# read landuse data
lusG_a = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,100],[21,39],[22,54],[23,61],[24,77],
[31,77],[41,30],[42,30],[43,30],[52,30],
[71,39],[81,39],[82,67],[90,100],
[95,100]]),"DATA")
lusG_b = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,100],[21,61],[22,70],[23,75],[24,85],
[31,86],[41,55],[42,55],[43,55],[52,48],
[71,61],[81,61],[82,78],[90,100],
[95,100]]),"DATA")
lusG_c = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,100],[21,74],[22,80],[23,83],[24,90],
[31,91],[41,70],[42,70],[43,70],[52,65],
[71,74],[81,74],[82,85],[90,100],
[95,100]]),"DATA")
lusG_d = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,100],[21,80],[22,85],[23,87],[24,92],
[31,94],[41,77],[42,77],[43,77],[52,77],
[71,80],[81,80],[82,89],[90,100],
[95,100]]),"DATA")
# Integer grid
numer = ((arcpy.sa.Times(opath + "/soilG_a",lusG_a))
+(arcpy.sa.Times(opath + "/soilG_b",lusG_b))
+(arcpy.sa.Times(opath + "/soilG_c",lusG_c))
+(arcpy.sa.Times(opath + "/soilG_d",lusG_d))
+(arcpy.sa.Times(opath + "/water_grid",100)))
denom = (arcpy.sa.Plus(opath + "/soilG_a",opath + "/soilG_b")
+(arcpy.sa.Plus(opath + "/soilG_c",opath + "/soilG_d"))
+(opath + "/water_grid"))
cnGrid = numer / denom
return cnGrid
def statsgonlcdlookuppoor(inLU,inSoil,opath):
arcpy.env.workspace = opath
arcpy.env.scratchWorkspace = opath + "/aux_folder"
# read statsgo shapefile
arcpy.FeatureToRaster_conversion(inSoil, "PCT_A", opath + "/soilG_a", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_B", opath + "/soilG_b", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_C", opath + "/soilG_c", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_D", opath + "/soilG_d", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_WATER", opath + "/water_grid", 30)
# read landuse data
lusG_a = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,100],[21,68],[22,76],[23,79],[24,88],
[31,77],[41,45],[42,45],[43,45],[52,48],
[71,68],[81,68],[82,72],[90,100],
[95,100]]),"DATA")
lusG_b = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,100],[21,79],[22,84],[23,86],[24,91],
[31,86],[41,66],[42,66],[43,66],[52,67],
[71,79],[81,79],[82,81],[90,100],
[95,100]]),"DATA")
lusG_c = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,100],[21,86],[22,89],[23,91],[24,94],
[31,91],[41,77],[42,77],[43,77],[52,77],
[71,86],[81,86],[82,86],[90,100],
[95,100]]),"DATA")
lusG_d = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,100],[21,89],[22,91],[23,92],[24,95],
[31,94],[41,83],[42,83],[43,83],[52,83],
[71,89],[81,89],[82,91],[90,100],
[95,100]]),"DATA")
# Integer grid
numer = ((arcpy.sa.Times(opath + "/soilG_a",lusG_a))
+(arcpy.sa.Times(opath + "/soilG_b",lusG_b))
+(arcpy.sa.Times(opath + "/soilG_c",lusG_c))
+(arcpy.sa.Times(opath + "/soilG_d",lusG_d))
+(arcpy.sa.Times(opath + "/water_grid",100)))
denom = (arcpy.sa.Plus(opath + "/soilG_a",opath + "/soilG_b")
+(arcpy.sa.Plus(opath + "/soilG_c",opath + "/soilG_d"))
+(opath + "/water_grid"))
cnGrid = numer / denom
return cnGrid
#*******************************************************************************************************
# Calculate curve number for fair, good, and poor soil conditions using STATSGO soil data
# for Anderson Look Up table
#*******************************************************************************************************
def statsgoandlookupfair(inLU,inSoil,opath):
arcpy.env.workspace = opath
arcpy.env.scratchWorkspace = opath + "/aux_folder"
# read statsgo shapefile
arcpy.FeatureToRaster_conversion(inSoil, "PCT_A", opath + "/soilG_a", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_B", opath + "/soilG_b", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_C", opath + "/soilG_c", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_D", opath + "/soilG_d", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_WATER", opath + "/water_grid", 30)
# read landuse data
lusG_a = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,68],[11,61],[12,68],[13,81],[14,91],
[15,81],[16,74],[17,77],[18,49],[20,70],
[21,70],[22,49],[23,43],[24,59],[25,70],
[40,36],[41,36],[42,36],[43,36],[44,35],
[50,100],[60,100],[70,77],[71,77],[72,77],
[73,77],[80,98],[191,70],[192,36],[241,59],
[242,59]]),"DATA")
lusG_b = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,80],[11,76],[12,80],[13,88],[14,94],
[15,88],[16,84],[17,86],[18,69],[20,80],
[21,80],[22,69],[23,65],[24,74],[25,80],
[40,60],[41,60],[42,60],[43,60],[44,56],
[50,100],[60,100],[70,86],[71,86],[72,86],
[73,86],[80,98],[191,80],[192,60],[241,74],
[242,74]]),"DATA")
lusG_c = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,86],[11,84],[12,86],[13,91],[14,95],
[15,91],[16,89],[17,91],[18,79],[20,87],
[21,87],[22,79],[23,76],[24,82],[25,87],
[40,73],[41,73],[42,73],[43,73],[44,70],
[50,100],[60,100],[70,91],[71,91],[72,91],
[73,91],[80,98],[191,87],[192,73],[241,82],
[242,82]]),"DATA")
lusG_d = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,89],[11,88],[12,89],[13,93],[14,96],
[15,93],[16,91],[17,94],[18,84],[20,90],
[21,90],[22,84],[23,82],[24,86],[25,90],
[40,79],[41,79],[42,79],[43,79],[44,77],
[50,100],[60,100],[70,94],[71,94],[72,94],
[73,94],[80,98],[191,90],[192,79],[241,86],
[242,86]]),"DATA")
# Integer grid
numer = ((arcpy.sa.Times(opath + "/soilG_a",lusG_a))
+(arcpy.sa.Times(opath + "/soilG_b",lusG_b))
+(arcpy.sa.Times(opath + "/soilG_c",lusG_c))
+(arcpy.sa.Times(opath + "/soilG_d",lusG_d))
+(arcpy.sa.Times(opath + "/water_grid",100)))
denom = (arcpy.sa.Plus(opath + "/soilG_a",opath + "/soilG_b")
+(arcpy.sa.Plus(opath + "/soilG_c",opath + "/soilG_d"))
+(opath + "/water_grid"))
cnGrid = numer / denom
return cnGrid
def statsgoandlookupgood(inLU,inSoil,opath):
arcpy.env.workspace = opath
arcpy.env.scratchWorkspace = opath + "/aux_folder"
# read statsgo shapefile
arcpy.FeatureToRaster_conversion(inSoil, "PCT_A", opath + "/soilG_a", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_B", opath + "/soilG_b", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_C", opath + "/soilG_c", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_D", opath + "/soilG_d", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_WATER", opath + "/water_grid", 30)
# read landuse data
lusG_a = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,61],[11,54],[12,61],[13,77],[14,89],
[15,81],[16,69],[17,77],[18,39],[20,67],
[21,67],[22,39],[23,32],[24,59],[25,67],
[40,30],[41,30],[42,30],[43,30],[44,30],
[50,100],[60,100],[70,77],[71,77],[72,77],
[73,77],[80,83],[191,67],[192,30],[241,59],
[242,59]]),"DATA")
lusG_b = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,75],[11,70],[12,75],[13,85],[14,92],
[15,88],[16,80],[17,86],[18,61],[20,78],
[21,78],[22,61],[23,58],[24,74],[25,78],
[40,55],[41,55],[42,55],[43,55],[44,48],
[50,100],[60,100],[70,86],[71,86],[72,86],
[73,86],[80,89],[191,78],[192,55],[241,74],
[242,74]]),"DATA")
lusG_c = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,83],[11,80],[12,83],[13,90],[14,94],
[15,91],[16,86],[17,91],[18,74],[20,85],
[21,85],[22,74],[23,72],[24,82],[25,85],
[40,70],[41,70],[42,70],[43,70],[44,65],
[50,100],[60,100],[70,91],[71,91],[72,91],
[73,91],[80,92],[191,85],[192,70],[241,82],
[242,82]]),"DATA")
lusG_d = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,87],[11,85],[12,87],[13,92],[14,95],
[15,93],[16,89],[17,94],[18,80],[20,89],
[21,89],[22,80],[23,79],[24,86],[25,89],
[40,77],[41,77],[42,77],[43,77],[44,73],
[50,100],[60,100],[70,94],[71,94],[72,94],
[73,94],[80,94],[191,89],[192,77],[241,86],
[242,86]]),"DATA")
# Integer grid
numer = ((arcpy.sa.Times(opath + "/soilG_a",lusG_a))
+(arcpy.sa.Times(opath + "/soilG_b",lusG_b))
+(arcpy.sa.Times(opath + "/soilG_c",lusG_c))
+(arcpy.sa.Times(opath + "/soilG_d",lusG_d))
+(arcpy.sa.Times(opath + "/water_grid",100)))
denom = (arcpy.sa.Plus(opath + "/soilG_a",opath + "/soilG_b")
+(arcpy.sa.Plus(opath + "/soilG_c",opath + "/soilG_d"))
+(opath + "/water_grid"))
cnGrid = numer / denom
return cnGrid
def statsgoandlookuppoor(inLU,inSoil,opath):
arcpy.env.workspace = opath
arcpy.env.scratchWorkspace = opath + "/aux_folder"
# read statsgo shapefile
arcpy.FeatureToRaster_conversion(inSoil, "PCT_A", opath + "/soilG_a", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_B", opath + "/soilG_b", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_C", opath + "/soilG_c", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_D", opath + "/soilG_d", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_WATER", opath + "/water_grid", 30)
# read landuse data
lusG_a = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,79],[11,76],[12,79],[13,88],[14,94],
[15,90],[16,83],[17,77],[18,68],[20,72],
[21,72],[22,68],[23,57],[24,59],[25,72],
[40,45],[41,45],[42,45],[43,45],[44,48],
[50,100],[60,100],[70,77],[71,77],[72,77],
[73,77],[80,98],[191,72],[192,45],[241,59],
[242,59]]),"DATA")
lusG_b = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,86],[11,84],[12,86],[13,91],[14,95],
[15,93],[16,89],[17,86],[18,79],[20,81],
[21,81],[22,79],[23,73],[24,74],[25,81],
[40,66],[41,66],[42,66],[43,66],[44,67],
[50,100],[60,100],[70,86],[71,86],[72,86],
[73,86],[80,98],[191,81],[192,66],[241,74],
[242,74]]),"DATA")
lusG_c = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,91],[11,89],[12,91],[13,94],[14,96],
[15,95],[16,92],[17,91],[18,86],[20,88],
[21,88],[22,86],[23,82],[24,82],[25,88],
[40,77],[41,77],[42,77],[43,77],[44,77],
[50,100],[60,100],[70,91],[71,91],[72,91],
[73,91],[80,98],[191,88],[192,77],[241,82],
[242,82]]),"DATA")
lusG_d = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,92],[11,91],[12,92],[13,95],[14,97],
[15,95],[16,94],[17,94],[18,89],[20,91],
[21,91],[22,89],[23,86],[24,86],[25,91],
[40,83],[41,83],[42,83],[43,83],[44,83],
[50,100],[60,100],[70,94],[71,94],[72,94],
[73,94],[80,98],[191,91],[192,83],[241,86],
[242,86]]),"DATA")
# Integer grid
numer = ((arcpy.sa.Times(opath + "/soilG_a",lusG_a))
+(arcpy.sa.Times(opath + "/soilG_b",lusG_b))
+(arcpy.sa.Times(opath + "/soilG_c",lusG_c))
+(arcpy.sa.Times(opath + "/soilG_d",lusG_d))
+(arcpy.sa.Times(opath + "/water_grid",100)))
denom = (arcpy.sa.Plus(opath + "/soilG_a",opath + "/soilG_b")
+(arcpy.sa.Plus(opath + "/soilG_c",opath + "/soilG_d"))
+(opath + "/water_grid"))
cnGrid = numer / denom
return cnGrid
#*******************************************************************************************************
# Calculate curve number for fair, good, and poor soil conditions using STATSGO soil data
# for MD/DE Look Up table
#*******************************************************************************************************
def statsgomddelookupfair(inLU,inSoil,opath):
arcpy.env.workspace = opath
arcpy.env.scratchWorkspace = opath + "/aux_folder"
# read statsgo shapefile
arcpy.FeatureToRaster_conversion(inSoil, "PCT_A", opath + "/soilG_a", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_B", opath + "/soilG_b", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_C", opath + "/soilG_c", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_D", opath + "/soilG_d", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_WATER", opath + "/water_grid", 30)
# read landuse data
lusG_a = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,68],[11,61],[12,68],[13,81],[14,91],
[15,81],[16,74],[17,77],[18,49],[20,70],
[21,70],[22,49],[23,43],[24,59],[25,70],
[40,36],[41,36],[42,36],[43,36],[44,35],
[50,100],[60,100],[70,77],[71,77],[72,77],
[73,77],[80,98],[191,70],[192,36],[241,59],
[242,59],[1110,68],[1120,81],[1140,81],[1200,91],
[1210,91],[1220,91],[1230,91],[1250,91],[1290,91],
[1300,81],[1400,98],[1410,98],[1420,91],[1430,98],
[1440,98],[1450,49],[1460,98],[1490,98],[1500,49],
[1600,74],[1700,74],[1800,74],[1900,49],[2100,70],
[2110,70],[2120,49],[2130,70],[2140,70],[2150,70],
[2200,43],[2300,59],[2400,59],[2900,70],[3100,36],
[3200,36],[3300,36],[4100,36],[4200,36],[4300,36],
[4400,36],[5100,100],[5200,100],[5300,100],[5400,100],
[5900,100],[6000,100],[7200,77],[7300,77],[7500,77],
[7600,77]]),"DATA")
lusG_b = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,80],[11,76],[12,80],[13,88],[14,94],
[15,88],[16,84],[17,86],[18,69],[20,80],
[21,80],[22,69],[23,65],[24,74],[25,80],
[40,60],[41,60],[42,60],[43,60],[44,56],
[50,100],[60,100],[70,86],[71,86],[72,86],
[73,86],[80,98],[191,80],[192,60],[241,74],
[242,74],[1110,80],[1120,88],[1140,88],[1200,94],
[1210,94],[1220,94],[1230,94],[1250,94],[1290,94],
[1300,88],[1400,98],[1410,98],[1420,94],[1430,98],
[1440,98],[1450,69],[1460,98],[1490,98],[1500,69],
[1600,84],[1700,84],[1800,84],[1900,69],[2100,80],
[2110,80],[2120,69],[2130,80],[2140,80],[2150,80],
[2200,65],[2300,74],[2400,74],[2900,80],[3100,60],
[3200,60],[3300,60],[4100,60],[4200,60],[4300,60],
[4400,60],[5100,100],[5200,100],[5300,100],[5400,100],
[5900,100],[6000,100],[7200,86],[7300,86],[7500,86],
[7600,86]]),"DATA")
lusG_c = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,86],[11,84],[12,86],[13,91],[14,95],
[15,91],[16,89],[17,91],[18,79],[20,87],
[21,87],[22,79],[23,76],[24,82],[25,87],
[40,73],[41,73],[42,73],[43,73],[44,70],
[50,100],[60,100],[70,91],[71,91],[72,91],
[73,91],[80,98],[191,87],[192,73],[241,82],
[242,82],[1110,86],[1120,91],[1140,91],[1200,95],
[1210,95],[1220,95],[1230,95],[1250,95],[1290,95],
[1300,91],[1400,98],[1410,98],[1420,95],[1430,98],
[1440,98],[1450,79],[1460,98],[1490,98],[1500,79],
[1600,89],[1700,89],[1800,89],[1900,79],[2100,87],
[2110,87],[2120,79],[2130,87],[2140,87],[2150,87],
[2200,76],[2300,82],[2400,82],[2900,87],[3100,73],
[3200,73],[3300,73],[4100,73],[4200,73],[4300,73],
[4400,73],[5100,100],[5200,100],[5300,100],[5400,100],
[5900,100],[6000,100],[7200,91],[7300,91],[7500,91],
[7600,91]]),"DATA")
lusG_d = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,89],[11,88],[12,89],[13,93],[14,96],
[15,93],[16,91],[17,94],[18,84],[20,90],
[21,90],[22,84],[23,82],[24,86],[25,90],
[40,79],[41,79],[42,79],[43,79],[44,77],
[50,100],[60,100],[70,94],[71,94],[72,94],
[73,94],[80,98],[191,90],[192,79],[241,86],
[242,86],[1110,89],[1120,93],[1140,93],[1200,96],
[1210,96],[1220,96],[1230,96],[1250,96],[1290,96],
[1300,93],[1400,98],[1410,98],[1420,96],[1430,98],
[1440,98],[1450,84],[1460,98],[1490,98],[1500,84],
[1600,91],[1700,91],[1800,91],[1900,84],[2100,90],
[2110,90],[2120,84],[2130,90],[2140,90],[2150,90],
[2200,82],[2300,86],[2400,86],[2900,90],[3100,79],
[3200,79],[3300,79],[4100,79],[4200,79],[4300,79],
[4400,79],[5100,100],[5200,100],[5300,100],[5400,100],
[5900,100],[6000,100],[7200,94],[7300,94],[7500,94],
[7600,94]]),"DATA")
# Integer grid
numer = ((arcpy.sa.Times(opath + "/soilG_a",lusG_a))
+(arcpy.sa.Times(opath + "/soilG_b",lusG_b))
+(arcpy.sa.Times(opath + "/soilG_c",lusG_c))
+(arcpy.sa.Times(opath + "/soilG_d",lusG_d))
+(arcpy.sa.Times(opath + "/water_grid",100)))
denom = (arcpy.sa.Plus(opath + "/soilG_a",opath + "/soilG_b")
+(arcpy.sa.Plus(opath + "/soilG_c",opath + "/soilG_d"))
+(opath + "/water_grid"))
cnGrid = numer / denom
return cnGrid
def statsgomddelookupgood(inLU,inSoil,opath):
arcpy.env.workspace = opath
arcpy.env.scratchWorkspace = opath + "/aux_folder"
# read statsgo shapefile
arcpy.FeatureToRaster_conversion(inSoil, "PCT_A", opath + "/soilG_a", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_B", opath + "/soilG_b", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_C", opath + "/soilG_c", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_D", opath + "/soilG_d", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_WATER", opath + "/water_grid", 30)
# read landuse data
lusG_a = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,61],[11,54],[12,61],[13,77],[14,89],
[15,81],[16,69],[17,77],[18,39],[20,67],
[21,67],[22,39],[23,32],[24,59],[25,67],
[40,30],[41,30],[42,30],[43,30],[44,30],
[50,100],[60,100],[70,77],[71,77],[72,77],
[73,77],[80,83],[191,67],[192,30],[241,59],
[242,59],[1110,61],[1120,77],[1140,77],[1200,89],
[1210,89],[1220,89],[1230,89],[1250,89],[1290,89],
[1300,81],[1400,83],[1410,83],[1420,89],[1430,83],
[1440,83],[1450,39],[1460,83],[1490,83],[1500,39],
[1600,69],[1700,69],[1800,69],[1900,39],[2100,67],
[2110,67],[2120,39],[2130,67],[2140,67],[2150,67],
[2200,32],[2300,59],[2400,59],[2900,67],[3100,30],
[3200,30],[3300,30],[4100,30],[4200,30],[4300,30],
[4400,30],[5100,100],[5200,100],[5300,100],[5400,100],
[5900,100],[6000,100],[7200,77],[7300,77],[7500,77],
[7600,77]]),"DATA")
lusG_b = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,75],[11,70],[12,75],[13,85],[14,92],
[15,88],[16,80],[17,86],[18,61],[20,78],
[21,78],[22,61],[23,58],[24,74],[25,78],
[40,55],[41,55],[42,55],[43,55],[44,48],
[50,100],[60,100],[70,86],[71,86],[72,86],
[73,86],[80,89],[191,78],[192,55],[241,74],
[242,74],[1110,75],[1120,85],[1140,85],[1200,92],
[1210,92],[1220,92],[1230,92],[1250,92],[1290,92],
[1300,88],[1400,89],[1410,89],[1420,92],[1430,89],
[1440,89],[1450,61],[1460,89],[1490,89],[1500,61],
[1600,80],[1700,80],[1800,80],[1900,61],[2100,78],
[2110,78],[2120,61],[2130,78],[2140,78],[2150,78],
[2200,58],[2300,74],[2400,74],[2900,78],[3100,55],
[3200,55],[3300,55],[4100,55],[4200,55],[4300,55],
[4400,55],[5100,100],[5200,100],[5300,100],[5400,100],
[5900,100],[6000,100],[7200,86],[7300,86],[7500,86],
[7600,86]]),"DATA")
lusG_c = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,83],[11,80],[12,83],[13,90],[14,94],
[15,91],[16,86],[17,91],[18,74],[20,85],
[21,85],[22,74],[23,72],[24,82],[25,85],
[40,70],[41,70],[42,70],[43,70],[44,65],
[50,100],[60,100],[70,91],[71,91],[72,91],
[73,91],[80,92],[191,85],[192,70],[241,82],
[242,82],[1110,83],[1120,90],[1140,90],[1200,94],
[1210,94],[1220,94],[1230,94],[1250,94],[1290,94],
[1300,91],[1400,92],[1410,92],[1420,94],[1430,92],
[1440,92],[1450,74],[1460,92],[1490,92],[1500,74],
[1600,86],[1700,86],[1800,86],[1900,74],[2100,85],
[2110,85],[2120,74],[2130,85],[2140,85],[2150,85],
[2200,72],[2300,82],[2400,82],[2900,85],[3100,70],
[3200,70],[3300,70],[4100,70],[4200,70],[4300,70],
[4400,70],[5100,100],[5200,100],[5300,100],[5400,100],
[5900,100],[6000,100],[7200,91],[7300,91],[7500,91],
[7600,91]]),"DATA")
lusG_d = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,87],[11,85],[12,87],[13,92],[14,95],
[15,93],[16,89],[17,94],[18,80],[20,89],
[21,89],[22,80],[23,79],[24,86],[25,89],
[40,77],[41,77],[42,77],[43,77],[44,73],
[50,100],[60,100],[70,94],[71,94],[72,94],
[73,94],[80,94],[191,89],[192,77],[241,86],
[242,86],[1110,87],[1120,92],[1140,92],[1200,95],
[1210,95],[1220,95],[1230,95],[1250,95],[1290,95],
[1300,93],[1400,94],[1410,94],[1420,95],[1430,94],
[1440,94],[1450,80],[1460,94],[1490,94],[1500,80],
[1600,89],[1700,89],[1800,89],[1900,80],[2100,89],
[2110,89],[2120,80],[2130,89],[2140,89],[2150,89],
[2200,79],[2300,86],[2400,86],[2900,89],[3100,77],
[3200,77],[3300,77],[4100,77],[4200,77],[4300,77],
[4400,77],[5100,100],[5200,100],[5300,100],[5400,100],
[5900,100],[6000,100],[7200,94],[7300,94],[7500,94],
[7600,94]]),"DATA")
# Integer grid
numer = ((arcpy.sa.Times(opath + "/soilG_a",lusG_a))
+(arcpy.sa.Times(opath + "/soilG_b",lusG_b))
+(arcpy.sa.Times(opath + "/soilG_c",lusG_c))
+(arcpy.sa.Times(opath + "/soilG_d",lusG_d))
+(arcpy.sa.Times(opath + "/water_grid",100)))
denom = (arcpy.sa.Plus(opath + "/soilG_a",opath + "/soilG_b")
+(arcpy.sa.Plus(opath + "/soilG_c",opath + "/soilG_d"))
+(opath + "/water_grid"))
cnGrid = numer / denom
return cnGrid
def statsgomddelookuppoor(inLU,inSoil,opath):
arcpy.env.workspace = opath
arcpy.env.scratchWorkspace = opath + "/aux_folder"
# read statsgo shapefile
arcpy.FeatureToRaster_conversion(inSoil, "PCT_A", opath + "/soilG_a", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_B", opath + "/soilG_b", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_C", opath + "/soilG_c", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_D", opath + "/soilG_d", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_WATER", opath + "/water_grid", 30)
# read landuse data
lusG_a = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,79],[11,76],[12,79],[13,88],[14,94],
[15,90],[16,83],[17,77],[18,68],[20,72],
[21,72],[22,68],[23,57],[24,59],[25,72],
[40,45],[41,45],[42,45],[43,45],[44,48],
[50,100],[60,100],[70,77],[71,77],[72,77],
[73,77],[80,98],[191,72],[192,45],[241,59],
[242,59],[1110,79],[1120,88],[1140,88],[1200,94],
[1210,94],[1220,94],[1230,94],[1250,94],[1290,94],
[1300,90],[1400,98],[1410,98],[1420,94],[1430,98],
[1440,98],[1450,68],[1460,98],[1490,98],[1500,68],
[1600,83],[1700,83],[1800,83],[1900,68],[2100,72],
[2110,72],[2120,68],[2130,72],[2140,72],[2150,72],
[2200,57],[2300,59],[2400,59],[2900,72],[3100,45],
[3200,45],[3300,45],[4100,45],[4200,45],[4300,45],
[4400,45],[5100,100],[5200,100],[5300,100],[5400,100],
[5900,100],[6000,100],[7200,77],[7300,77],[7500,77],
[7600,77]]),"DATA")
lusG_b = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,86],[11,84],[12,86],[13,91],[14,95],
[15,93],[16,89],[17,86],[18,79],[20,81],
[21,81],[22,79],[23,73],[24,74],[25,81],
[40,66],[41,66],[42,66],[43,66],[44,67],
[50,100],[60,100],[70,86],[71,86],[72,86],
[73,86],[80,98],[191,81],[192,66],[241,74],
[242,74],[1110,86],[1120,91],[1140,91],[1200,95],
[1210,95],[1220,95],[1230,95],[1250,95],[1290,95],
[1300,93],[1400,98],[1410,98],[1420,95],[1430,98],
[1440,98],[1450,79],[1460,98],[1490,98],[1500,79],
[1600,89],[1700,89],[1800,89],[1900,79],[2100,81],
[2110,81],[2120,79],[2130,81],[2140,81],[2150,81],
[2200,73],[2300,74],[2400,74],[2900,81],[3100,66],
[3200,66],[3300,66],[4100,66],[4200,66],[4300,66],
[4400,66],[5100,100],[5200,100],[5300,100],[5400,100],
[5900,100],[6000,100],[7200,86],[7300,86],[7500,86],
[7600,86]]),"DATA")
lusG_c = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,91],[11,89],[12,91],[13,94],[14,96],
[15,95],[16,92],[17,91],[18,86],[20,88],
[21,88],[22,86],[23,82],[24,82],[25,88],
[40,77],[41,77],[42,77],[43,77],[44,77],
[50,100],[60,100],[70,91],[71,91],[72,91],
[73,91],[80,98],[191,88],[192,77],[241,82],
[242,82],[1110,91],[1120,94],[1140,94],[1200,96],
[1210,96],[1220,96],[1230,96],[1250,96],[1290,96],
[1300,95],[1400,98],[1410,98],[1420,96],[1430,98],
[1440,98],[1450,86],[1460,98],[1490,98],[1500,86],
[1600,92],[1700,92],[1800,92],[1900,86],[2100,88],
[2110,88],[2120,86],[2130,88],[2140,88],[2150,88],
[2200,82],[2300,82],[2400,82],[2900,88],[3100,77],
[3200,77],[3300,77],[4100,77],[4200,77],[4300,77],
[4400,77],[5100,100],[5200,100],[5300,100],[5400,100],
[5900,100],[6000,100],[7200,91],[7300,91],[7500,91],
[7600,91]]),"DATA")
lusG_d = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,92],[11,91],[12,92],[13,95],[14,97],
[15,95],[16,94],[17,94],[18,89],[20,91],
[21,91],[22,89],[23,86],[24,86],[25,91],
[40,83],[41,83],[42,83],[43,83],[44,83],
[50,100],[60,100],[70,94],[71,94],[72,94],
[73,94],[80,98],[191,91],[192,83],[241,86],
[242,86],[1110,91],[1120,95],[1140,95],[1200,97],
[1210,97],[1220,97],[1230,97],[1250,97],[1290,97],
[1300,95],[1400,98],[1410,98],[1420,97],[1430,98],
[1440,98],[1450,89],[1460,98],[1490,98],[1500,89],
[1600,94],[1700,94],[1800,94],[1900,89],[2100,91],
[2110,91],[2120,89],[2130,91],[2140,91],[2150,91],
[2200,86],[2300,86],[2400,86],[2900,91],[3100,83],
[3200,83],[3300,83],[4100,83],[4200,83],[4300,83],
[4400,83],[5100,100],[5200,100],[5300,100],[5400,100],
[5900,100],[6000,100],[7200,94],[7300,94],[7500,94],
[7600,94]]),"DATA")
# Integer grid
numer = ((arcpy.sa.Times(opath + "/soilG_a",lusG_a))
+(arcpy.sa.Times(opath + "/soilG_b",lusG_b))
+(arcpy.sa.Times(opath + "/soilG_c",lusG_c))
+(arcpy.sa.Times(opath + "/soilG_d",lusG_d))
+(arcpy.sa.Times(opath + "/water_grid",100)))
denom = (arcpy.sa.Plus(opath + "/soilG_a",opath + "/soilG_b")
+(arcpy.sa.Plus(opath + "/soilG_c",opath + "/soilG_d"))
+(opath + "/water_grid"))
cnGrid = numer / denom
return cnGrid
#*******************************************************************************************************
# Calculate curve number for fair, good, and poor soil conditions using STATSGO soil data
# for MRLC Look Up table
#*******************************************************************************************************
def statsgomrlclookupfair(inLU,inSoil,opath):
arcpy.env.workspace = opath
arcpy.env.scratchWorkspace = opath + "/aux_folder"
# read statsgo shapefile
arcpy.FeatureToRaster_conversion(inSoil, "PCT_A", opath + "/soilG_a", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_B", opath + "/soilG_b", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_C", opath + "/soilG_c", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_D", opath + "/soilG_d", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_WATER", opath + "/water_grid", 30)
# read landuse data
lusG_a = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,100],[21,61],[22,81],[23,91],[24,91],
[31,77],[32,77],[33,77],[41,36],[42,36],
[43,36],[81,49],[82,70],[85,49],[90,100],
[91,100],[92,100],[95,100]]),"DATA")
lusG_b = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,100],[21,76],[22,88],[23,94],[24,94],
[31,86],[32,86],[33,86],[41,60],[42,60],
[43,60],[81,69],[82,80],[85,69],[90,100],
[91,100],[92,100],[95,100]]),"DATA")
lusG_c = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,100],[21,84],[22,91],[23,95],[24,95],
[31,91],[32,91],[33,91],[41,73],[42,73],
[43,73],[81,79],[82,87],[85,79],[90,100],
[91,100],[92,100],[95,100]]),"DATA")
lusG_d = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,100],[21,88],[22,93],[23,96],[24,96],
[31,94],[32,94],[33,94],[41,79],[42,79],
[43,79],[81,84],[82,90],[85,84],[90,100],
[91,100],[92,100],[95,100]]),"DATA")
# Integer grid
numer = ((arcpy.sa.Times(opath + "/soilG_a",lusG_a))
+(arcpy.sa.Times(opath + "/soilG_b",lusG_b))
+(arcpy.sa.Times(opath + "/soilG_c",lusG_c))
+(arcpy.sa.Times(opath + "/soilG_d",lusG_d))
+(arcpy.sa.Times(opath + "/water_grid",100)))
denom = (arcpy.sa.Plus(opath + "/soilG_a",opath + "/soilG_b")
+(arcpy.sa.Plus(opath + "/soilG_c",opath + "/soilG_d"))
+(opath + "/water_grid"))
cnGrid = numer / denom
return cnGrid
def statsgomrlclookupgood(inLU,inSoil,opath):
arcpy.env.workspace = opath
arcpy.env.scratchWorkspace = opath + "/aux_folder"
# read statsgo shapefile
arcpy.FeatureToRaster_conversion(inSoil, "PCT_A", opath + "/soilG_a", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_B", opath + "/soilG_b", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_C", opath + "/soilG_c", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_D", opath + "/soilG_d", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_WATER", opath + "/water_grid", 30)
# read landuse data
lusG_a = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,100],[21,54],[22,77],[23,89],[24,89],
[31,77],[32,77],[33,77],[41,30],[42,30],
[43,30],[81,39],[82,67],[85,39],[90,100],
[91,100],[92,100],[95,100]]),"DATA")
lusG_b = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,100],[21,70],[22,85],[23,92],[24,92],
[31,86],[32,86],[33,86],[41,55],[42,55],
[43,55],[81,61],[82,78],[85,61],[90,100],
[91,100],[92,100],[95,100]]),"DATA")
lusG_c = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,100],[21,80],[22,90],[23,94],[24,94],
[31,91],[32,91],[33,91],[41,70],[42,70],
[43,70],[81,74],[82,85],[85,74],[90,100],
[91,100],[92,100],[95,100]]),"DATA")
lusG_d = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,100],[21,85],[22,92],[23,95],[24,95],
[31,94],[32,94],[33,94],[41,77],[42,77],
[43,77],[81,80],[82,89],[85,80],[90,100],
[91,100],[92,100],[95,100]]),"DATA")
# Integer grid
numer = ((arcpy.sa.Times(opath + "/soilG_a",lusG_a))
+(arcpy.sa.Times(opath + "/soilG_b",lusG_b))
+(arcpy.sa.Times(opath + "/soilG_c",lusG_c))
+(arcpy.sa.Times(opath + "/soilG_d",lusG_d))
+(arcpy.sa.Times(opath + "/water_grid",100)))
denom = (arcpy.sa.Plus(opath + "/soilG_a",opath + "/soilG_b")
+(arcpy.sa.Plus(opath + "/soilG_c",opath + "/soilG_d"))
+(opath + "/water_grid"))
cnGrid = numer / denom
return cnGrid
def statsgomrlclookuppoor(inLU,inSoil,opath):
arcpy.env.workspace = opath
arcpy.env.scratchWorkspace = opath + "/aux_folder"
# read statsgo shapefile
arcpy.FeatureToRaster_conversion(inSoil, "PCT_A", opath + "/soilG_a", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_B", opath + "/soilG_b", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_C", opath + "/soilG_c", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_D", opath + "/soilG_d", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_WATER", opath + "/water_grid", 30)
# read landuse data
lusG_a = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,100],[21,76],[22,88],[23,89],[24,94],
[31,77],[32,77],[33,77],[41,45],[42,45],
[43,45],[81,68],[82,72],[85,68],[90,100],
[91,100],[92,100],[95,100]]),"DATA")
lusG_b = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,100],[21,84],[22,91],[23,92],[24,95],
[31,86],[32,86],[33,86],[41,66],[42,66],
[43,66],[81,79],[82,81],[85,79],[90,100],
[91,100],[92,100],[95,100]]),"DATA")
lusG_c = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,100],[21,89],[22,94],[23,94],[24,96],
[31,91],[32,91],[33,91],[41,77],[42,77],
[43,77],[81,86],[82,88],[85,86],[90,100],
[91,100],[92,100],[95,100]]),"DATA")
lusG_d = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,100],[21,91],[22,95],[23,95],[24,97],
[31,94],[32,94],[33,94],[41,83],[42,83],
[43,83],[81,89],[82,91],[85,89],[90,100],
[91,100],[92,100],[95,100]]),"DATA")
# Integer grid
numer = ((arcpy.sa.Times(opath + "/soilG_a",lusG_a))
+(arcpy.sa.Times(opath + "/soilG_b",lusG_b))
+(arcpy.sa.Times(opath + "/soilG_c",lusG_c))
+(arcpy.sa.Times(opath + "/soilG_d",lusG_d))
+(arcpy.sa.Times(opath + "/water_grid",100)))
denom = (arcpy.sa.Plus(opath + "/soilG_a",opath + "/soilG_b")
+(arcpy.sa.Plus(opath + "/soilG_c",opath + "/soilG_d"))
+(opath + "/water_grid"))
cnGrid = numer / denom
return cnGrid
#*******************************************************************************************************
# Calculate curve number for fair, good, and poor soil conditions using STATSGO soil data
# for USGS Look Up table
#*******************************************************************************************************
def statsgousgslookupfair(inLU,inSoil,opath):
arcpy.env.workspace = opath
arcpy.env.scratchWorkspace = opath + "/aux_folder"
# read statsgo shapefile
arcpy.FeatureToRaster_conversion(inSoil, "PCT_A", opath + "/soilG_a", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_B", opath + "/soilG_b", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_C", opath + "/soilG_c", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_D", opath + "/soilG_d", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_WATER", opath + "/water_grid", 30)
# read landuse data
lusG_a = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,68],[12,91],[13,81],[14,98],[15,91],
[16,68],[17,54],[21,70],[22,43],[23,59],
[24,70],[41,36],[42,36],[43,36],[51,100],
[52,100],[53,100],[54,100],[61,100],[62,100],
[70,77],[71,100],[72,77],[73,77],[74,100],
[75,77],[76,77],[77,77]]),"DATA")
lusG_b = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,80],[12,94],[13,88],[14,98],[15,94],
[16,80],[17,72],[21,80],[22,65],[23,74],
[24,80],[41,60],[42,60],[43,60],[51,100],
[52,100],[53,100],[54,100],[61,100],[62,100],
[70,86],[71,100],[72,86],[73,86],[74,100],
[75,86],[76,86],[77,86]]),"DATA")
lusG_c = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,86],[12,95],[13,91],[14,98],[15,95],
[16,86],[17,81],[21,87],[22,76],[23,82],
[24,87],[41,73],[42,73],[43,73],[51,100],
[52,100],[53,100],[54,100],[61,100],[62,100],
[70,91],[71,100],[72,91],[73,91],[74,100],
[75,91],[76,91],[77,91]]),"DATA")
lusG_d = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,89],[12,96],[13,93],[14,98],[15,96],
[16,89],[17,86],[21,90],[22,82],[23,86],
[24,90],[41,79],[42,79],[43,79],[51,100],
[52,100],[53,100],[54,100],[61,100],[62,100],
[70,94],[71,100],[72,94],[73,94],[74,100],
[75,94],[76,94],[77,94]]),"DATA")
# Integer grid
numer = ((arcpy.sa.Times(opath + "/soilG_a",lusG_a))
+(arcpy.sa.Times(opath + "/soilG_b",lusG_b))
+(arcpy.sa.Times(opath + "/soilG_c",lusG_c))
+(arcpy.sa.Times(opath + "/soilG_d",lusG_d))
+(arcpy.sa.Times(opath + "/water_grid",100)))
denom = (arcpy.sa.Plus(opath + "/soilG_a",opath + "/soilG_b")
+(arcpy.sa.Plus(opath + "/soilG_c",opath + "/soilG_d"))
+(opath + "/water_grid"))
cnGrid = numer / denom
return cnGrid
def statsgousgslookupgood(inLU,inSoil,opath):
arcpy.env.workspace = opath
arcpy.env.scratchWorkspace = opath + "/aux_folder"
# read statsgo shapefile
arcpy.FeatureToRaster_conversion(inSoil, "PCT_A", opath + "/soilG_a", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_B", opath + "/soilG_b", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_C", opath + "/soilG_c", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_D", opath + "/soilG_d", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_WATER", opath + "/water_grid", 30)
# read landuse data
lusG_a = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,61],[12,89],[13,81],[14,98],[15,89],
[16,61],[17,45],[21,67],[22,32],[23,59],
[24,67],[41,30],[42,30],[43,30],[51,100],
[52,100],[53,100],[54,100],[61,100],[62,100],
[70,77],[71,100],[72,77],[73,77],[74,77],
[75,77],[76,77],[77,77]]),"DATA")
lusG_b = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,75],[12,92],[13,88],[14,98],[15,92],
[16,75],[17,65],[21,78],[22,58],[23,74],
[24,78],[41,55],[42,55],[43,55],[51,100],
[52,100],[53,100],[54,100],[61,100],[62,100],
[70,86],[71,100],[72,86],[73,86],[74,86],
[75,86],[76,86],[77,86]]),"DATA")
lusG_c = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,83],[12,94],[13,91],[14,98],[15,94],
[16,83],[17,77],[21,85],[22,72],[23,82],
[24,85],[41,70],[42,70],[43,70],[51,100],
[52,100],[53,100],[54,100],[61,100],[62,100],
[70,91],[71,100],[72,91],[73,91],[74,91],
[75,91],[76,91],[77,91]]),"DATA")
lusG_d = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,87],[12,95],[13,93],[14,98],[15,95],
[16,87],[17,82],[21,89],[22,79],[23,86],
[24,89],[41,77],[42,77],[43,77],[51,100],
[52,100],[53,100],[54,100],[61,100],[62,100],
[70,94],[71,100],[72,94],[73,94],[74,94],
[75,94],[76,94],[77,94]]),"DATA")
# Integer grid
numer = ((arcpy.sa.Times(opath + "/soilG_a",lusG_a))
+(arcpy.sa.Times(opath + "/soilG_b",lusG_b))
+(arcpy.sa.Times(opath + "/soilG_c",lusG_c))
+(arcpy.sa.Times(opath + "/soilG_d",lusG_d))
+(arcpy.sa.Times(opath + "/water_grid",100)))
denom = (arcpy.sa.Plus(opath + "/soilG_a",opath + "/soilG_b")
+(arcpy.sa.Plus(opath + "/soilG_c",opath + "/soilG_d"))
+(opath + "/water_grid"))
cnGrid = numer / denom
return cnGrid
def statsgousgslookuppoor(inLU,inSoil,opath):
arcpy.env.workspace = opath
arcpy.env.scratchWorkspace = opath + "/aux_folder"
# read statsgo shapefile
arcpy.FeatureToRaster_conversion(inSoil, "PCT_A", opath + "/soilG_a", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_B", opath + "/soilG_b", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_C", opath + "/soilG_c", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_D", opath + "/soilG_d", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_WATER", opath + "/water_grid", 30)
# read landuse data
lusG_a = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,79],[12,94],[13,90],[14,98],[15,94],
[16,79],[17,71],[21,72],[22,57],[23,59],
[24,72],[41,45],[42,45],[43,45],[51,100],
[52,100],[53,100],[54,100],[61,100],[62,100],
[70,77],[71,100],[72,77],[73,77],[74,100],
[75,77],[76,77],[77,77]]),"DATA")
lusG_b = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,86],[12,95],[13,93],[14,98],[15,95],
[16,86],[17,81],[21,81],[22,73],[23,74],
[24,81],[41,66],[42,66],[43,66],[51,100],
[52,100],[53,100],[54,100],[61,100],[62,100],
[70,86],[71,100],[72,86],[73,86],[74,100],
[75,86],[76,86],[77,86]]),"DATA")
lusG_c = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,91],[12,96],[13,95],[14,98],[15,96],
[16,91],[17,87],[21,88],[22,82],[23,82],
[24,88],[41,77],[42,77],[43,77],[51,100],
[52,100],[53,100],[54,100],[61,100],[62,100],
[70,91],[71,100],[72,91],[73,91],[74,100],
[75,91],[76,91],[77,91]]),"DATA")
lusG_d = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,92],[12,97],[13,95],[14,98],[15,97],
[16,92],[17,90],[21,91],[22,86],[23,86],
[24,91],[41,83],[42,83],[43,83],[51,100],
[52,100],[53,100],[54,100],[61,100],[62,100],
[70,94],[71,100],[72,94],[73,94],[74,100],
[75,94],[76,94],[77,94]]),"DATA")
# Integer grid
numer = ((arcpy.sa.Times(opath + "/soilG_a",lusG_a))
+(arcpy.sa.Times(opath + "/soilG_b",lusG_b))
+(arcpy.sa.Times(opath + "/soilG_c",lusG_c))
+(arcpy.sa.Times(opath + "/soilG_d",lusG_d))
+(arcpy.sa.Times(opath + "/water_grid",100)))
denom = (arcpy.sa.Plus(opath + "/soilG_a",opath + "/soilG_b")
+(arcpy.sa.Plus(opath + "/soilG_c",opath + "/soilG_d"))
+(opath + "/water_grid"))
cnGrid = numer / denom
return cnGrid
#*******************************************************************************************************
# Calculate curve number for fair, good, and poor soil conditions using STATSGO soil data
# for Zoning Look Up table
#*******************************************************************************************************
def statsgozoninglookupfair(inLU,inSoil,opath):
arcpy.env.workspace = opath
arcpy.env.scratchWorkspace = opath + "/aux_folder"
# read statsgo shapefile
arcpy.FeatureToRaster_conversion(inSoil, "PCT_A", opath + "/soilG_a", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_B", opath + "/soilG_b", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_C", opath + "/soilG_c", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_D", opath + "/soilG_d", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_WATER", opath + "/water_grid", 30)
# read landuse data
lusG_a = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,68],[11,61],[12,68],[13,81],[14,91],
[15,81],[16,74],[17,77],[18,49],[20,70],
[21,70],[22,49],[23,43],[24,59],[25,70],
[40,36],[41,36],[42,36],[43,36],[44,35],
[50,100],[60,100],[70,77],[71,77],[72,77],
[73,77],[80,98],[111,55],[112,59],[113,61],
[114,64],[115,68],[116,81],[191,70],[192,36],
[241,59],[242,59],[1110,68],[1120,81],[1140,81],
[1200,91],[1210,91],[1220,91],[1230,91],[1250,91],
[1290,91],[1300,81],[1400,98],[1410,98],[1420,91],
[1430,98],[1440,98],[1450,49],[1460,98],[1490,98],
[1500,49],[1600,74],[1700,74],[1800,74],[1900,49],
[2100,70],[2110,70],[2120,49],[2130,70],[2140,70],
[2150,70],[2200,43],[2300,59],[2400,59],[2900,70],
[3100,36],[3200,36],[3300,36],[4100,36],[4200,36],
[4300,36],[4400,36],[5100,100],[5200,100],[5300,100],
[5400,100],[5900,100],[6000,100],[7200,77],[7300,77],
[7500,77],[7600,77]]),"DATA")
lusG_b = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,80],[11,76],[12,80],[13,88],[14,94],
[15,88],[16,84],[17,86],[18,69],[20,80],
[21,80],[22,69],[23,65],[24,74],[25,80],
[40,60],[41,60],[42,60],[43,60],[44,56],
[50,100],[60,100],[70,86],[71,86],[72,86],
[73,86],[80,98],[111,72],[112,75],[113,76],
[114,78],[115,80],[116,88],[191,80],[192,60],
[241,74],[242,74],[1110,80],[1120,88],[1140,88],
[1200,94],[1210,94],[1220,94],[1230,94],[1250,94],
[1290,94],[1300,88],[1400,98],[1410,98],[1420,94],
[1430,98],[1440,98],[1450,69],[1460,98],[1490,98],
[1500,69],[1600,84],[1700,84],[1800,84],[1900,69],
[2100,80],[2110,80],[2120,69],[2130,80],[2140,80],
[2150,80],[2200,65],[2300,74],[2400,74],[2900,80],
[3100,60],[3200,60],[3300,60],[4100,60],[4200,60],
[4300,60],[4400,60],[5100,100],[5200,100],[5300,100],
[5400,100],[5900,100],[6000,100],[7200,86],[7300,86],
[7500,86],[7600,86]]),"DATA")
lusG_c = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,86],[11,84],[12,86],[13,91],[14,95],
[15,91],[16,89],[17,91],[18,79],[20,87],
[21,87],[22,79],[23,76],[24,82],[25,87],
[40,73],[41,73],[42,73],[43,73],[44,70],
[50,100],[60,100],[70,91],[71,91],[72,91],
[73,91],[80,98],[111,81],[112,83],[113,84],
[114,85],[115,86],[116,91],[191,87],[192,73],
[241,82],[242,82],[1110,86],[1120,91],[1140,91],
[1200,95],[1210,95],[1220,95],[1230,95],[1250,95],
[1290,95],[1300,91],[1400,98],[1410,98],[1420,95],
[1430,98],[1440,98],[1450,79],[1460,98],[1490,98],
[1500,79],[1600,89],[1700,89],[1800,89],[1900,79],
[2100,87],[2110,87],[2120,79],[2130,87],[2140,87],
[2150,87],[2200,76],[2300,82],[2400,82],[2900,87],
[3100,73],[3200,73],[3300,73],[4100,73],[4200,73],
[4300,73],[4400,73],[5100,100],[5200,100],[5300,100],
[5400,100],[5900,100],[6000,100],[7200,91],[7300,91],
[7500,91],[7600,91]]),"DATA")
lusG_d = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,89],[11,88],[12,89],[13,93],[14,96],
[15,93],[16,91],[17,94],[18,84],[20,90],
[21,90],[22,84],[23,82],[24,86],[25,90],
[40,79],[41,79],[42,79],[43,79],[44,77],
[50,100],[60,100],[70,94],[71,94],[72,94],
[73,94],[80,98],[111,86],[112,87],[113,88],
[114,88],[115,89],[116,93],[191,90],[192,79],
[241,86],[242,86],[1110,89],[1120,93],[1140,93],
[1200,96],[1210,96],[1220,96],[1230,96],[1250,96],
[1290,96],[1300,93],[1400,98],[1410,98],[1420,96],
[1430,98],[1440,98],[1450,84],[1460,98],[1490,98],
[1500,84],[1600,91],[1700,91],[1800,91],[1900,84],
[2100,90],[2110,90],[2120,84],[2130,90],[2140,90],
[2150,90],[2200,82],[2300,86],[2400,86],[2900,90],
[3100,79],[3200,79],[3300,79],[4100,79],[4200,79],
[4300,79],[4400,79],[5100,100],[5200,100],[5300,100],
[5400,100],[5900,100],[6000,100],[7200,94],[7300,94],
[7500,94],[7600,94]]),"DATA")
# Integer grid
numer = ((arcpy.sa.Times(opath + "/soilG_a",lusG_a))
+(arcpy.sa.Times(opath + "/soilG_b",lusG_b))
+(arcpy.sa.Times(opath + "/soilG_c",lusG_c))
+(arcpy.sa.Times(opath + "/soilG_d",lusG_d))
+(arcpy.sa.Times(opath + "/water_grid",100)))
denom = (arcpy.sa.Plus(opath + "/soilG_a",opath + "/soilG_b")
+(arcpy.sa.Plus(opath + "/soilG_c",opath + "/soilG_d"))
+(opath + "/water_grid"))
cnGrid = numer / denom
return cnGrid
def statsgozoninglookupgood(inLU,inSoil,opath):
arcpy.env.workspace = opath
arcpy.env.scratchWorkspace = opath + "/aux_folder"
# read statsgo shapefile
arcpy.FeatureToRaster_conversion(inSoil, "PCT_A", opath + "/soilG_a", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_B", opath + "/soilG_b", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_C", opath + "/soilG_c", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_D", opath + "/soilG_d", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_WATER", opath + "/water_grid", 30)
# read landuse data
lusG_a = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,61],[11,54],[12,61],[13,77],[14,89],
[15,81],[16,69],[17,77],[18,39],[20,67],
[21,67],[22,39],[23,32],[24,59],[25,67],
[40,30],[41,30],[42,30],[43,30],[44,30],
[50,100],[60,100],[70,77],[71,77],[72,77],
[73,77],[80,98],[111,46],[112,51],[113,54],
[114,57],[115,61],[116,77],[191,67],[192,30],
[241,59],[242,59],[1110,61],[1120,77],[1140,77],
[1200,89],[1210,89],[1220,89],[1230,89],[1250,89],
[1290,89],[1300,81],[1400,83],[1410,83],[1420,89],
[1430,83],[1440,83],[1450,39],[1460,83],[1490,83],
[1500,39],[1600,69],[1700,69],[1800,69],[1900,39],
[2100,67],[2110,67],[2120,39],[2130,67],[2140,67],
[2150,67],[2200,32],[2300,59],[2400,59],[2900,67],
[3100,30],[3200,30],[3300,30],[4100,30],[4200,30],
[4300,30],[4400,30],[5100,100],[5200,100],[5300,100],
[5400,100],[5900,100],[6000,100],[7200,77],[7300,77],
[7500,77],[7600,77]]),"DATA")
lusG_b = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,75],[11,70],[12,75],[13,85],[14,92],
[15,88],[16,80],[17,86],[18,61],[20,78],
[21,78],[22,61],[23,58],[24,74],[25,78],
[40,55],[41,55],[42,55],[43,55],[44,48],
[50,100],[60,100],[70,86],[71,86],[72,86],
[73,86],[80,98],[111,65],[112,68],[113,70],
[114,72],[115,75],[116,85],[191,78],[192,55],
[241,74],[242,74],[1110,75],[1120,85],[1140,85],
[1200,92],[1210,92],[1220,92],[1230,92],[1250,92],
[1290,92],[1300,88],[1400,99],[1410,89],[1420,92],
[1430,89],[1440,89],[1450,61],[1460,89],[1490,89],
[1500,61],[1600,80],[1700,80],[1800,80],[1900,61],
[2100,78],[2110,78],[2120,61],[2130,78],[2140,78],
[2150,78],[2200,58],[2300,74],[2400,74],[2900,78],
[3100,55],[3200,55],[3300,55],[4100,55],[4200,55],
[4300,55],[4400,55],[5100,100],[5200,100],[5300,100],
[5400,100],[5900,100],[6000,100],[7200,86],[7300,86],
[7500,86],[7600,86]]),"DATA")
lusG_c = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,83],[11,80],[12,83],[13,90],[14,94],
[15,91],[16,86],[17,91],[18,74],[20,85],
[21,85],[22,74],[23,72],[24,82],[25,85],
[40,70],[41,70],[42,70],[43,70],[44,65],
[50,100],[60,100],[70,91],[71,91],[72,91],
[73,91],[80,98],[111,77],[112,79],[113,80],
[114,81],[115,83],[116,90],[191,85],[192,70],
[241,82],[242,82],[1110,83],[1120,90],[1140,90],
[1200,94],[1210,94],[1220,94],[1230,94],[1250,94],
[1290,94],[1300,91],[1400,92],[1410,92],[1420,94],
[1430,92],[1440,92],[1450,74],[1460,92],[1490,92],
[1500,74],[1600,86],[1700,86],[1800,86],[1900,74],
[2100,85],[2110,85],[2120,74],[2130,85],[2140,85],
[2150,85],[2200,72],[2300,82],[2400,82],[2900,85],
[3100,70],[3200,70],[3300,70],[4100,70],[4200,70],
[4300,70],[4400,70],[5100,100],[5200,100],[5300,100],
[5400,100],[5900,100],[6000,100],[7200,91],[7300,91],
[7500,91],[7600,91]]),"DATA")
lusG_d = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,87],[11,85],[12,87],[13,92],[14,95],
[15,93],[16,89],[17,94],[18,80],[20,89],
[21,89],[22,80],[23,79],[24,86],[25,89],
[40,77],[41,77],[42,77],[43,77],[44,73],
[50,100],[60,100],[70,94],[71,94],[72,94],
[73,94],[80,98],[111,82],[112,84],[113,85],
[114,86],[115,87],[116,92],[191,89],[192,77],
[241,86],[242,86],[1110,87],[1120,92],[1140,92],
[1200,95],[1210,95],[1220,95],[1230,95],[1250,95],
[1290,95],[1300,93],[1400,94],[1410,94],[1420,95],
[1430,94],[1440,94],[1450,80],[1460,94],[1490,94],
[1500,80],[1600,89],[1700,89],[1800,89],[1900,80],
[2100,89],[2110,89],[2120,80],[2130,89],[2140,89],
[2150,89],[2200,79],[2300,86],[2400,86],[2900,89],
[3100,77],[3200,77],[3300,77],[4100,77],[4200,77],
[4300,77],[4400,77],[5100,100],[5200,100],[5300,100],
[5400,100],[5900,100],[6000,100],[7200,94],[7300,94],
[7500,94],[7600,94]]),"DATA")
# Integer grid
numer = ((arcpy.sa.Times(opath + "/soilG_a",lusG_a))
+(arcpy.sa.Times(opath + "/soilG_b",lusG_b))
+(arcpy.sa.Times(opath + "/soilG_c",lusG_c))
+(arcpy.sa.Times(opath + "/soilG_d",lusG_d))
+(arcpy.sa.Times(opath + "/water_grid",100)))
denom = (arcpy.sa.Plus(opath + "/soilG_a",opath + "/soilG_b")
+(arcpy.sa.Plus(opath + "/soilG_c",opath + "/soilG_d"))
+(opath + "/water_grid"))
cnGrid = numer / denom
return cnGrid
def statsgozoninglookuppoor(inLU,inSoil,opath):
arcpy.env.workspace = opath
arcpy.env.scratchWorkspace = opath + "/aux_folder"
# read statsgo shapefile
arcpy.FeatureToRaster_conversion(inSoil, "PCT_A", opath + "/soilG_a", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_B", opath + "/soilG_b", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_C", opath + "/soilG_c", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_D", opath + "/soilG_d", 30)
arcpy.FeatureToRaster_conversion(inSoil, "PCT_WATER", opath + "/water_grid", 30)
# read landuse data
lusG_a = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,79],[11,76],[12,79],[13,88],[14,94],
[15,90],[16,83],[17,77],[18,68],[20,72],
[21,72],[22,68],[23,57],[24,59],[25,72],
[40,45],[41,45],[42,45],[43,45],[44,48],
[50,100],[60,100],[70,77],[71,77],[72,77],
[73,77],[80,98],[111,72],[112,74],[113,76],
[114,77],[115,79],[116,88],[191,72],[192,45],
[241,59],[242,59],[1110,79],[1120,88],[1140,88],
[1200,94],[1210,94],[1220,94],[1230,94],[1250,94],
[1290,94],[1300,90],[1400,98],[1410,98],[1420,94],
[1430,98],[1440,98],[1450,68],[1460,98],[1490,98],
[1500,68],[1600,83],[1700,83],[1800,83],[1900,68],
[2100,72],[2110,72],[2120,68],[2130,72],[2140,72],
[2150,72],[2200,57],[2300,59],[2400,59],[2900,72],
[3100,45],[3200,45],[3300,45],[4100,45],[4200,45],
[4300,45],[4400,45],[5100,100],[5200,100],[5300,100],
[5400,100],[5900,100],[6000,100],[7200,77],[7300,77],
[7500,77],[7600,77]]),"DATA")
lusG_b = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,86],[11,84],[12,86],[13,91],[14,95],
[15,93],[16,89],[17,86],[18,79],[20,81],
[21,81],[22,79],[23,73],[24,74],[25,81],
[40,66],[41,66],[42,66],[43,66],[44,67],
[50,100],[60,100],[70,86],[71,86],[72,86],
[73,86],[80,98],[111,81],[112,83],[113,84],
[114,85],[115,86],[116,91],[191,81],[192,66],
[241,74],[242,74],[1110,86],[1120,91],[1140,91],
[1200,95],[1210,95],[1220,95],[1230,95],[1250,95],
[1290,95],[1300,93],[1400,98],[1410,98],[1420,95],
[1430,98],[1440,98],[1450,79],[1460,98],[1490,98],
[1500,79],[1600,89],[1700,89],[1800,89],[1900,79],
[2100,81],[2110,81],[2120,79],[2130,81],[2140,81],
[2150,81],[2200,73],[2300,74],[2400,74],[2900,81],
[3100,66],[3200,66],[3300,66],[4100,66],[4200,66],
[4300,66],[4400,66],[5100,100],[5200,100],[5300,100],
[5400,100],[5900,100],[6000,100],[7200,86],[7300,86],
[7500,86],[7600,86]]),"DATA")
lusG_c = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,91],[11,89],[12,91],[13,94],[14,96],
[15,95],[16,92],[17,91],[18,86],[20,88],
[21,88],[22,86],[23,82],[24,82],[25,88],
[40,77],[41,77],[42,77],[43,77],[44,77],
[50,100],[60,100],[70,91],[71,91],[72,91],
[73,91],[80,98],[111,87],[112,88],[113,89],
[114,90],[115,91],[116,94],[191,88],[192,77],
[241,82],[242,82],[1110,91],[1120,94],[1140,94],
[1200,96],[1210,96],[1220,96],[1230,96],[1250,96],
[1290,96],[1300,95],[1400,98],[1410,98],[1420,96],
[1430,98],[1440,98],[1450,86],[1460,98],[1490,98],
[1500,86],[1600,92],[1700,92],[1800,92],[1900,86],
[2100,86],[2110,88],[2120,86],[2130,88],[2140,88],
[2150,88],[2200,82],[2300,82],[2400,82],[2900,88],
[3100,77],[3200,77],[3300,77],[4100,77],[4200,77],
[4300,77],[4400,77],[5100,100],[5200,100],[5300,100],
[5400,100],[5900,100],[6000,100],[7200,91],[7300,91],
[7500,91],[7600,91]]),"DATA")
lusG_d = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,92],[11,91],[12,92],[13,95],[14,97],
[15,95],[16,94],[17,94],[18,89],[20,91],
[21,91],[22,89],[23,86],[24,86],[25,91],
[40,83],[41,83],[42,83],[43,83],[44,83],
[50,100],[60,100],[70,94],[71,94],[72,94],
[73,94],[80,98],[111,90],[112,91],[113,91],
[114,92],[115,92],[116,95],[191,91],[192,83],
[241,86],[242,86],[1110,92],[1120,95],[1140,95],
[1200,97],[1210,97],[1220,97],[1230,97],[1250,97],
[1290,97],[1300,95],[1400,98],[1410,98],[1420,97],
[1430,98],[1440,98],[1450,98],[1460,98],[1490,98],
[1500,89],[1600,94],[1700,94],[1800,94],[1900,89],
[2100,91],[2110,91],[2120,89],[2130,91],[2140,91],
[2150,91],[2200,86],[2300,86],[2400,86],[2900,91],
[3100,83],[3200,83],[3300,83],[4100,83],[4200,83],
[4300,83],[4400,83],[5100,100],[5200,100],[5300,100],
[5400,100],[5900,100],[6000,100],[7200,94],[7300,94],
[7500,94],[7600,94]]),"DATA")
# Integer grid
numer = ((arcpy.sa.Times(opath + "/soilG_a",lusG_a))
+(arcpy.sa.Times(opath + "/soilG_b",lusG_b))
+(arcpy.sa.Times(opath + "/soilG_c",lusG_c))
+(arcpy.sa.Times(opath + "/soilG_d",lusG_d))
+(arcpy.sa.Times(opath + "/water_grid",100)))
denom = (arcpy.sa.Plus(opath + "/soilG_a",opath + "/soilG_b")
+(arcpy.sa.Plus(opath + "/soilG_c",opath + "/soilG_d"))
+(opath + "/water_grid"))
cnGrid = numer / denom
return cnGrid
#*******************************************************************************************************
# Calculate soil percentages using selected soil data for soil percentages
#*******************************************************************************************************
##pctAsoilR = 0
##pctBsoilR = 0
##pctCsoilR = 0
##pctDsoilR = 0
##pctWsoilR = 0
def SoilPct(basinarea,inSoil):
global pctAsoilR
global pctBsoilR
global pctCsoilR
global pctDsoilR
global pctWsoilR
pctAsoilR = 0
pctBsoilR = 0
pctCsoilR = 0
pctDsoilR = 0
pctWsoilR = 0
temptab = arcpy.SearchCursor(inSoil,"","","Value;Count","")
for row in temptab:
thecount = float(row.getValue("Count"))
if row.getValue("Value") == 1:
## pythonaddins.MessageBox(type(pctAsoilR),"pctAsoilR",0)
pctAsoilR = pctAsoilR + thecount
elif row.getValue("Value") == 2:
pctBsoilR = pctBsoilR + thecount
elif row.getValue("Value") == 3:
pctCsoilR = pctCsoilR + thecount
elif row.getValue("Value") == 4:
pctDsoilR = pctDsoilR + thecount
elif row.getValue("Value") == -1:
pctWsoilR = pctWsoilR + thecount
del row
return pctAsoilR, pctBsoilR, pctCsoilR, pctDsoilR, pctWsoilR
#*******************************************************************************************************
# Calculate soil percentages using SSURGO for use in regression equations
#*******************************************************************************************************
# soil types are initialized with zero values before function arguments to avoid following error:
# UnboundLocalError: local variable referenced before assignment
#
def SSURGOPct(basinarea,ssurgo):
global pctAsoil
global pctBsoil
global pctCsoil
global pctDsoil
global pctWsoil
pctAsoil = 0
pctBsoil = 0
pctCsoil = 0
pctDsoil = 0
pctWsoil = 0
temptab = arcpy.SearchCursor(ssurgo,"","","Value;Count","")
for row in temptab:
count = float(row.getValue("Count"))
if row.getValue("Value") == 1:
pctAsoil = pctAsoil + float((count/basinarea)* 100)
elif row.getValue("Value") == 2:
pctBsoil = pctBsoil + float((count/basinarea)* 100)
elif row.getValue("Value") == 3:
pctCsoil = pctCsoil + float((count/basinarea)* 100)
elif row.getValue("Value") == 4:
pctDsoil = pctDsoil + float((count/basinarea)* 100)
elif row.getValue("Value") == -1:
pctWsoil = pctWsoil + float((count/basinarea)* 100)
del row
return pctAsoil, pctBsoil, pctCsoil, pctDsoil, pctWsoil
#*******************************************************************************************************
# Calculate LU cell count for delineated watershed -- LU count will vary based on landuse type
#*******************************************************************************************************
def GetLUCountNLCD(basinarea,inLU):
global Urb
global nil
global FC
global ST
Urb = 0
nil = 0
FC = 0
ST = 0
# urban = 1, n = 2, forest = 3, storage = 4
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,4],[21,2],[22,1],[23,1],[24,1],
[31,2],[41,3],[42,3],[43,3],[52,3],
[71,2],[81,2],[82,2],[90,4],
[95,4]]),"DATA")
#rows = arcpy.da.SearchCursor(LU,"Value;Count") # input arg order changed with v10.1
# BUT its better to use simple search cursor
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
count = float(row.getValue("Count"))
if row.getValue("Value") == 1:
Urb = Urb + float((count/basinarea)* 100)
elif row.getValue("Value") == 2:
nil = nil + float((count/basinarea)* 100)
elif row.getValue("Value") == 3:
FC = FC + float((count/basinarea)* 100)
elif row.getValue("Value") == 4:
ST = ST + float((count/basinarea)* 100)
del row
return Urb, nil, FC, ST # map these values and index out as required
def GetLUCountAnderson(basinarea,inLU):
global Urb
global nil
global FC
global ST
Urb = 0
nil = 0
FC = 0
ST = 0
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,1],[11,1],[12,1],[13,1],[14,1],
[15,1],[16,2],[17,2],[18,2],[20,2],
[21,2],[22,2],[23,2],[24,2],[25,2],
[40,3],[41,3],[42,3],[43,3],[44,3],
[50,4],[60,4],[70,2],[71,2],[72,2],
[73,2],[80,2],[191,2],[192,3],[241,2],
[242,2]]),"DATA")
#rows = arcpy.da.SearchCursor(LU,"Value;Count") # input arg order changed with v10.1
# BUT its better to use simple search cursor
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
count = float(row.getValue("Count"))
if row.getValue("Value") == 1:
Urb = Urb + float((count/basinarea)* 100)
elif row.getValue("Value") == 2:
nil = nil + float((count/basinarea)* 100)
elif row.getValue("Value") == 3:
FC = FC + float((count/basinarea)* 100)
elif row.getValue("Value") == 4:
ST = ST + float((count/basinarea)* 100)
del row
return Urb, nil, FC, ST # map these values and index out as required
def GetLUCountMDDE(basinarea,inLU):
global Urb
global nil
global FC
global ST
Urb = 0
nil = 0
FC = 0
ST = 0
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,1],[11,1],[12,1],[13,1],[14,1],
[15,1],[16,2],[17,2],[18,2],[20,2],
[21,2],[22,2],[23,2],[24,2],[25,2],
[40,3],[41,3],[42,3],[43,3],[44,3],
[50,4],[60,4],[70,2],[71,2],[72,2],
[73,2],[80,2],[191,2],[192,3],[241,2],
[242,2],[1110,1],[1120,1],[1140,1],
[1200,1],[1210,1],[1220,1],[1230,1],
[1250,1],[1290,1],[1300,1],[1400,2],
[1410,2],[1420,1],[1430,2],[1440,2],
[1450,2],[1460,2],[1490,2],[1500,2],
[1600,2],[1700,2],[1800,2],[1900,2],
[2100,2],[2110,2],[2120,2],[2130,2],
[2140,2],[2150,2],[2200,2],[2300,2],
[2400,2],[2900,2],[3100,3],[3200,3],
[3300,3],[4100,3],[4200,3],[4300,3],
[4400,3],[5100,4],[5200,4],[5300,4],
[5400,4],[5900,4],[6000,4],[7200,2],
[7300,2],[7500,2],[7600,2]]),"DATA")
#rows = arcpy.da.SearchCursor(LU,"Value;Count") # input arg order changed with v10.1
# BUT its better to use simple search cursor
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
count = float(row.getValue("Count"))
if row.getValue("Value") == 1:
Urb = Urb + float((count/basinarea)* 100)
elif row.getValue("Value") == 2:
nil = nil + float((count/basinarea)* 100)
elif row.getValue("Value") == 3:
FC = FC + float((count/basinarea)* 100)
elif row.getValue("Value") == 4:
ST = ST + float((count/basinarea)* 100)
del row
return Urb, nil, FC, ST # map these values and index out as required
def GetLUCountUltimate(basinarea,inLU):
global Urb
global nil
global FC
global ST
Urb = 0
nil = 0
FC = 0
ST = 0
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,1],[11,1],[12,1],[13,1],[14,1],
[15,1],[16,2],[17,2],[18,2],[20,2],
[21,2],[22,2],[23,2],[24,2],[25,2],
[40,3],[41,3],[42,3],[43,3],[44,3],
[50,4],[60,4],[70,2],[71,2],[72,2],
[73,2],[80,2],[111,2],[112,2],[113,1],
[114,1],[115,1],[116,1],[191,2],[192,3],
[241,2],[242,2],[1110,1],[1120,1],[1140,1],
[1200,1],[1210,1],[1220,1],[1230,1],
[1250,1],[1290,1],[1300,1],[1400,2],
[1410,2],[1420,1],[1430,2],[1440,2],
[1450,2],[1460,2],[1490,2],[1500,2],
[1600,2],[1700,2],[1800,2],[1900,2],
[2100,2],[2110,2],[2120,2],[2130,2],
[2140,2],[2150,2],[2200,2],[2300,2],
[2400,2],[2900,2],[3100,3],[3200,3],
[3300,3],[4100,3],[4200,3],[4300,3],
[4400,3],[5100,4],[5200,4],[5300,4],
[5400,4],[5900,4],[6000,4],[7200,2],
[7300,2],[7500,2],[7600,2]]),"DATA")
#rows = arcpy.da.SearchCursor(LU,"Value;Count") # input arg order changed with v10.1
# BUT its better to use simple search cursor
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
count = float(row.getValue("Count"))
if row.getValue("Value") == 1:
Urb = Urb + float((count/basinarea)* 100)
elif row.getValue("Value") == 2:
nil = nil + float((count/basinarea)* 100)
elif row.getValue("Value") == 3:
FC = FC + float((count/basinarea)* 100)
elif row.getValue("Value") == 4:
ST = ST + float((count/basinarea)* 100)
del row
return Urb, nil, FC, ST # map these values and index out as required
def GetLUCountMRLC(basinarea,inLU):
global Urb
global nil
global FC
global ST
Urb = 0
nil = 0
FC = 0
ST = 0
# urban = 1, n = 2, forest = 3, storage = 4
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,4],[21,1],[22,1],[23,1],[24,1],
[31,2],[32,2],[33,2],[41,3],[42,3],
[43,3],[81,2],[82,2],[85,2],[90,4],
[91,4],[92,4],[95,4]]),"DATA")
#rows = arcpy.da.SearchCursor(LU,"Value;Count") # input arg order changed with v10.1
# BUT its better to use simple search cursor
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
count = float(row.getValue("Count"))
if row.getValue("Value") == 1:
Urb = Urb + float((count/basinarea)* 100)
elif row.getValue("Value") == 2:
nil = nil + float((count/basinarea)* 100)
elif row.getValue("Value") == 3:
FC = FC + float((count/basinarea)* 100)
elif row.getValue("Value") == 4:
ST = ST + float((count/basinarea)* 100)
del row
return Urb, nil, FC, ST # map these values and index out as required
def GetLUCountUSGS(basinarea,inLU):
global Urb
global nil
global FC
global ST
Urb = 0
nil = 0
FC = 0
ST = 0
# urban = 1, n = 2, forest = 3, storage = 4
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,1],[12,1],[13,1],[14,1],[15,1],
[16,2],[17,2],[21,2],[22,2],[23,2],
[24,2],[41,3],[42,3],[43,3],[51,4],
[52,4],[53,4],[54,4],[61,4],[62,4],
[70,2],[71,2],[72,2],[73,2],[74,2],
[75,2],[76,2],[77,2]]),"DATA")
#rows = arcpy.da.SearchCursor(LU,"Value;Count") # input arg order changed with v10.1
# BUT its better to use simple search cursor
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
count = float(row.getValue("Count"))
if row.getValue("Value") == 1:
Urb = Urb + float((count/basinarea)* 100)
elif row.getValue("Value") == 2:
nil = nil + float((count/basinarea)* 100)
elif row.getValue("Value") == 3:
FC = FC + float((count/basinarea)* 100)
elif row.getValue("Value") == 4:
ST = ST + float((count/basinarea)* 100)
del row
return Urb, nil, FC, ST # map these values and index out as required
#*******************************************************************************************************
# Reclassify LU data to corresponding Imp values from NLCD table for Impervious area calculation
#*******************************************************************************************************
##Impcount = 0
def GetImpCountNLCDFair(inLU):
global Impcount
Impcount = 0
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,0],[21,11],[22,25],[23,38],[24,65],
[31,50],[41,0],[42,0],[43,0],[52,0],
[71,0],[81,0],[82,0],[90,0],
[95,0]]),"DATA")
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
val = float(row.getValue("Value"))
value = float(val/100)
count = row.getValue("Count")
# Positive Impervious count
if value > 0:
Impcount = Impcount + float(value* count)
return Impcount
##Impcount = 0
def GetImpCountNLCDGood(inLU):
global Impcount
Impcount = 0
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,0],[21,11],[22,25],[23,38],[24,65],
[31,50],[41,0],[42,0],[43,0],[52,0],
[71,0],[81,0],[82,0],[90,0],
[95,0]]),"DATA")
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
val = float(row.getValue("Value"))
value = float(val/100)
count = row.getValue("Count")
# Positive Impervious count
if value > 0:
Impcount = Impcount + float(value* count)
return Impcount
def GetImpCountNLCDPoor(inLU):
global Impcount
Impcount = 0
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,0],[21,11],[22,25],[23,38],[24,65],
[31,50],[41,0],[42,0],[43,0],[52,0],
[71,0],[81,0],[82,0],[90,0],
[95,0]]),"DATA")
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
val = float(row.getValue("Value"))
value = float(val/100)
count = row.getValue("Count")
# Positive Impervious count
if value > 0:
Impcount = Impcount + float(value* count)
return Impcount
#*******************************************************************************************************
# Reclassify LU data to corresponding Imp values from Anderson table for Impervious area calculation
#*******************************************************************************************************
##Impcount = 0
def GetImpCountAndersonFair(inLU):
global Impcount
Impcount = 0
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,38],[11,25],[12,38],[13,65],[14,85],[15,72],
[16,50],[17,11],[18,11],[20,0],[21,0],[22,0],
[23,0],[24,0],[25,0],[40,0],[41,0],[42,0],
[43,0],[44,0],[50,0],[60,0],[70,50],[71,0],
[72,100],[73,50],[80,100],[191,15],[192,15],
[241,10],[242,10]]),"DATA")
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
val = float(row.getValue("Value"))
value = float(val/100)
count = row.getValue("Count")
# Positive Impervious count
if value > 0:
Impcount = Impcount + float(value* count)
return Impcount
##Impcount = 0
def GetImpCountAndersonGood(inLU):
global Impcount
Impcount = 0
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,38],[11,25],[12,38],[13,65],[14,85],[15,72],
[16,50],[17,11],[18,11],[20,0],[21,0],[22,0],
[23,0],[24,0],[25,0],[40,0],[41,0],[42,0],
[43,0],[44,0],[50,0],[60,0],[70,50],[71,0],
[72,100],[73,50],[80,75],[191,15],[192,15],
[241,10],[242,10]]),"DATA")
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
val = float(row.getValue("Value"))
value = float(val/100)
count = row.getValue("Count")
# Positive Impervious count
if value > 0:
Impcount = Impcount + float(value* count)
return Impcount
##Impcount = 0
def GetImpCountAndersonPoor(inLU):
global Impcount
Impcount = 0
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,38],[11,25],[12,38],[13,65],[14,85],[15,72],
[16,50],[17,11],[18,11],[20,0],[21,0],[22,0],
[23,0],[24,0],[25,0],[40,0],[41,0],[42,0],
[43,0],[44,0],[50,0],[60,0],[70,50],[71,0],
[72,100],[73,50],[80,100],[191,15],[192,15],
[241,10],[242,10]]),"DATA")
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
val = float(row.getValue("Value"))
value = float(val/100)
count = row.getValue("Count")
# Positive Impervious count
if value > 0:
Impcount = Impcount + float(value* count)
return Impcount
#*******************************************************************************************************
# Reclassify LU data to corresponding Imp values from MDDE table for Impervious area calculation
#*******************************************************************************************************
##Impcount = 0
def GetImpCountMDDEFair(inLU):
global Impcount
Impcount = 0
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,38],[11,25],[12,38],[13,65],[14,85],[15,72],
[16,50],[17,11],[18,11],[20,0],[21,0],[22,0],
[23,0],[24,0],[25,0],[40,0],[41,0],[42,0],
[43,0],[44,0],[50,0],[60,0],[70,50],[71,0],
[72,100],[73,50],[80,100],[191,15],[192,15],
[241,10],[242,10],[1110,38],[1120,65],[1140,65],
[1200,85],[1210,85],[1220,85],[1230,85],
[1250,85],[1290,85],[1300,72],[1400,100],
[1410,100],[1420,85],[1430,100],[1440,100],
[1450,11],[1460,100],[1490,100],[1500,11],
[1600,50],[1700,50],[1800,50],[1900,11],
[2100,0],[2110,0],[2120,0],[2130,0],[2140,0],
[2150,0],[2200,0],[2300,10],[2400,10],[2900,0],
[3100,0],[3200,0],[3300,0],[4100,0],[4200,0],
[4300,0],[4400,0],[5100,0],[5200,0],[5300,0],
[5400,0],[5900,0],[6000,0],[7200,0],[7300,0],
[7500,11],[7600,50]]),"DATA")
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
val = float(row.getValue("Value"))
value = float(val/100)
count = row.getValue("Count")
# Positive Impervious count
if value > 0:
Impcount = Impcount + float(value* count)
return Impcount
def GetImpCountMDDEGood(inLU):
global Impcount
Impcount = 0
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,38],[11,25],[12,38],[13,65],[14,85],[15,72],
[16,50],[17,11],[18,11],[20,0],[21,0],[22,0],
[23,0],[24,0],[25,0],[40,0],[41,0],[42,0],
[43,0],[44,0],[50,0],[60,0],[70,50],[71,0],
[72,100],[73,50],[80,75],[191,15],[192,15],
[241,10],[242,10],[1110,38],[1120,65],[1140,65],
[1200,85],[1210,85],[1220,85],[1230,85],
[1250,85],[1290,85],[1300,72],[1400,75],
[1410,75],[1420,85],[1430,75],[1440,75],
[1450,11],[1460,75],[1490,75],[1500,11],
[1600,50],[1700,50],[1800,50],[1900,11],
[2100,0],[2110,0],[2120,0],[2130,0],[2140,0],
[2150,0],[2200,0],[2300,10],[2400,10],[2900,0],
[3100,0],[3200,0],[3300,0],[4100,0],[4200,0],
[4300,0],[4400,0],[5100,0],[5200,0],[5300,0],
[5400,0],[5900,0],[6000,0],[7200,0],[7300,0],
[7500,11],[7600,50]]),"DATA")
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
val = float(row.getValue("Value"))
value = float(val/100)
count = row.getValue("Count")
# Positive Impervious count
if value > 0:
Impcount = Impcount + float(value* count)
return Impcount
def GetImpCountMDDEPoor(inLU):
global Impcount
Impcount = 0
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,38],[11,25],[12,38],[13,65],[14,85],[15,72],
[16,50],[17,11],[18,11],[20,0],[21,0],[22,0],
[23,0],[24,0],[25,0],[40,0],[41,0],[42,0],
[43,0],[44,0],[50,0],[60,0],[70,50],[71,0],
[72,100],[73,50],[80,100],[191,15],[192,15],
[241,10],[242,10],[1110,38],[1120,65],[1140,65],
[1200,85],[1210,85],[1220,85],[1230,85],
[1250,85],[1290,85],[1300,72],[1400,100],
[1410,100],[1420,85],[1430,100],[1440,100],
[1450,11],[1460,100],[1490,100],[1500,11],
[1600,50],[1700,50],[1800,50],[1900,11],
[2100,0],[2110,0],[2120,0],[2130,0],[2140,0],
[2150,0],[2200,0],[2300,10],[2400,10],[2900,0],
[3100,0],[3200,0],[3300,0],[4100,0],[4200,0],
[4300,0],[4400,0],[5100,0],[5200,0],[5300,0],
[5400,0],[5900,0],[6000,0],[7200,0],[7300,0],
[7500,11],[7600,50]]),"DATA")
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
val = float(row.getValue("Value"))
value = float(val/100)
count = row.getValue("Count")
# Positive Impervious count
if value > 0:
Impcount = Impcount + float(value* count)
return Impcount
#*******************************************************************************************************
# Reclassify LU data to corresponding Imp values from Ultimate table for Impervious area calculation
#*******************************************************************************************************
##Impcount = 0
def GetImpCountUltimateFair(inLU):
global Impcount
Impcount = 0
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,38],[11,25],[12,38],[13,65],[14,85],[15,72],
[16,50],[17,11],[18,11],[20,0],[21,0],[22,0],
[23,0],[24,0],[25,0],[40,0],[41,0],[42,0],
[43,0],[44,0],[50,0],[60,0],[70,50],[71,0],
[72,100],[73,50],[80,100],[111,12],[112,20],
[113,25],[114,30],[115,38],[116,65],[191,15],
[192,15],[241,10],[242,10],[1110,38],[1120,65],
[1140,65],[1200,85],[1210,85],[1220,85],[1230,85],
[1250,85],[1290,85],[1300,72],[1400,100],
[1410,100],[1420,85],[1430,100],[1440,100],
[1450,11],[1460,100],[1490,100],[1500,11],
[1600,50],[1700,50],[1800,50],[1900,11],
[2100,0],[2110,0],[2120,0],[2130,0],[2140,0],
[2150,0],[2200,0],[2300,10],[2400,10],[2900,0],
[3100,0],[3200,0],[3300,0],[4100,0],[4200,0],
[4300,0],[4400,0],[5100,0],[5200,0],[5300,0],
[5400,0],[5900,0],[6000,0],[7200,0],[7300,0],
[7500,11],[7600,50]]),"DATA")
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
val = float(row.getValue("Value"))
value = float(val/100)
count = row.getValue("Count")
# Positive Impervious count
if value > 0:
Impcount = Impcount + float(value* count)
return Impcount
def GetImpCountUltimateGood(inLU):
global Impcount
Impcount = 0
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,38],[11,25],[12,38],[13,65],[14,85],[15,72],
[16,50],[17,11],[18,11],[20,0],[21,0],[22,0],
[23,0],[24,0],[25,0],[40,0],[41,0],[42,0],
[43,0],[44,0],[50,0],[60,0],[70,50],[71,0],
[72,100],[73,50],[80,100],[111,12],[112,20],
[113,25],[114,30],[115,38],[116,65],[191,15],
[192,15],[241,10],[242,10],[1110,38],[1120,65],
[1140,65],[1200,85],[1210,85],[1220,85],[1230,85],
[1250,85],[1290,85],[1300,72],[1400,75],
[1410,75],[1420,85],[1430,75],[1440,75],
[1450,11],[1460,75],[1490,75],[1500,11],
[1600,50],[1700,50],[1800,50],[1900,11],
[2100,0],[2110,0],[2120,0],[2130,0],[2140,0],
[2150,0],[2200,0],[2300,10],[2400,10],[2900,0],
[3100,0],[3200,0],[3300,0],[4100,0],[4200,0],
[4300,0],[4400,0],[5100,0],[5200,0],[5300,0],
[5400,0],[5900,0],[6000,0],[7200,0],[7300,0],
[7500,11],[7600,50]]),"DATA")
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
val = float(row.getValue("Value"))
value = float(val/100)
count = row.getValue("Count")
# Positive Impervious count
if value > 0:
Impcount = Impcount + float(value* count)
return Impcount
def GetImpCountUltimatePoor(inLU):
global Impcount
Impcount = 0
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[10,38],[11,25],[12,38],[13,65],[14,85],[15,72],
[16,50],[17,11],[18,11],[20,0],[21,0],[22,0],
[23,0],[24,0],[25,0],[40,0],[41,0],[42,0],
[43,0],[44,0],[50,0],[60,0],[70,50],[71,0],
[72,100],[73,50],[80,100],[111,12],[112,20],
[113,25],[114,30],[115,38],[116,65],[191,15],
[192,15],[241,10],[242,10],[1110,38],[1120,65],
[1140,65],[1200,85],[1210,85],[1220,85],[1230,85],
[1250,85],[1290,85],[1300,72],[1400,100],
[1410,100],[1420,85],[1430,100],[1440,100],
[1450,11],[1460,100],[1490,100],[1500,11],
[1600,50],[1700,50],[1800,50],[1900,11],
[2100,0],[2110,0],[2120,0],[2130,0],[2140,0],
[2150,0],[2200,0],[2300,10],[2400,10],[2900,0],
[3100,0],[3200,0],[3300,0],[4100,0],[4200,0],
[4300,0],[4400,0],[5100,0],[5200,0],[5300,0],
[5400,0],[5900,0],[6000,0],[7200,0],[7300,0],
[7500,11],[7600,50]]),"DATA")
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
val = float(row.getValue("Value"))
value = float(val/100)
count = row.getValue("Count")
# Positive Impervious count
if value > 0:
Impcount = Impcount + float(value* count)
return Impcount
#*******************************************************************************************************
# Reclassify LU data to corresponding Imp values from MRLC table for Impervious area calculation
#*******************************************************************************************************
##Impcount = 0
def GetImpCountMRLCFair(inLU):
global Impcount
Impcount = 0
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,0],[21,25],[22,65],[23,85],[24,85],
[31,100],[32,11],[33,50],[41,0],[42,0],
[43,0],[81,0],[82,0],[85,11],[90,0],[91,0],
[92,0],[95,0]]),"DATA")
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
val = float(row.getValue("Value"))
value = float(val/100)
count = row.getValue("Count")
# Positive Impervious count
if value > 0:
Impcount = Impcount + float(value* count)
return Impcount
def GetImpCountMRLCGood(inLU):
global Impcount
Impcount = 0
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,0],[21,25],[22,65],[23,85],[24,85],
[31,100],[32,11],[33,50],[41,0],[42,0],
[43,0],[81,0],[82,0],[85,11],[90,0],[91,0],
[92,0],[95,0]]),"DATA")
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
val = float(row.getValue("Value"))
value = float(val/100)
count = row.getValue("Count")
# Positive Impervious count
if value > 0:
Impcount = Impcount + float(value* count)
return Impcount
def GetImpCountMRLCPoor(inLU):
global Impcount
Impcount = 0
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,0],[21,25],[22,65],[23,82],[24,85],
[31,100],[32,11],[33,50],[41,0],[42,0],
[43,0],[81,0],[82,0],[85,11],[90,0],[91,0],
[92,0],[95,0]]),"DATA")
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
val = float(row.getValue("Value"))
value = float(val/100)
count = row.getValue("Count")
# Positive Impervious count
if value > 0:
Impcount = Impcount + float(value* count)
return Impcount
#*******************************************************************************************************
# Reclassify LU data to corresponding Imp values from USGS table for Impervious area calculation
#*******************************************************************************************************
##Impcount = 0
def GetImpCountUSGSFair(inLU):
global Impcount
Impcount = 0
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,38],[12,85],[13,72],[14,100],[15,85],
[16,38],[17,11],[21,0],[22,0],
[23,0],[24,0],[41,0],[42,0],
[43,0],[51,0],[52,0],[53,0],[54,0],[61,0],
[62,0],[70,50],[71,100],[72,0],[73,0],[74,100],
[75,11],[76,50],[77,50]]),"DATA")
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
val = float(row.getValue("Value"))
value = float(val/100)
count = row.getValue("Count")
# Positive Impervious count
if value > 0:
Impcount = Impcount + float(value* count)
return Impcount
def GetImpCountUSGSGood(inLU):
global Impcount
Impcount = 0
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,38],[12,85],[13,72],[14,100],[15,85],
[16,38],[17,11],[21,0],[22,0],
[23,0],[24,0],[41,0],[42,0],
[43,0],[51,0],[52,0],[53,0],[54,0],[61,0],
[62,0],[70,50],[71,100],[72,0],[73,0],[74,100],
[75,11],[76,50],[77,50]]),"DATA")
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
val = float(row.getValue("Value"))
value = float(val/100)
count = row.getValue("Count")
# Positive Impervious count
if value > 0:
Impcount = Impcount + float(value* count)
return Impcount
def GetImpCountUSGSPoor(inLU):
global Impcount
Impcount = 0
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,38],[12,85],[13,72],[14,100],[15,85],
[16,38],[17,11],[21,0],[22,0],
[23,0],[24,0],[41,0],[42,0],
[43,0],[51,0],[52,0],[53,0],[54,0],[61,0],
[62,0],[70,50],[71,100],[72,0],[73,0],[74,100],
[75,11],[76,50],[77,50]]),"DATA")
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
val = float(row.getValue("Value"))
value = float(val/100)
count = row.getValue("Count")
# Positive Impervious count
if value > 0:
Impcount = Impcount + float(value* count)
return Impcount
#*******************************************************************************************************
# Get FC for subwatersheds based on landuse data type
#*******************************************************************************************************
def TcFC(inLU):
rows = arcpy.SearchCursor(inLU,"","","Value;Count","")
FCcnt = 0
for row in rows:
if (row.getValue("Value") >= 40) & (row.getValue("Value") < 44):
thecount = row.getValue("Count")
FCcnt = FCcnt + thecount
return FCcnt
#*******************************************************************************************************
# Get ST for subwatersheds based on landuse data type
#*******************************************************************************************************
def TcST_MRLC(inLU):
rows = arcpy.SearchCursor(inLU,"","","Value;Count","")
STcnt = 0
for row in rows:
if (row.getValue("Value") == 11) or (row.getValue("Value") == 91) or (row.getValue("Value") == 92):
thecount = row.getValue("Count")
STcnt = STcnt + thecount
return STcnt
def TcST(inLU):
rows = arcpy.SearchCursor(inLU,"","","Value;Count","")
STcnt = 0
for row in rows:
if (row.getValue("Value") >= 50) & (row.getValue("Value") <= 62):
thecount = row.getValue("Count")
STcnt = STcnt + thecount
return STcnt
#*******************************************************************************************************
# Reclassify LU data to corresponding Imp values to compute Tc -- repeat it for all different landuse
# types
#*******************************************************************************************************
def TcImpAnderson(inLU):
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,0.25],[12,0.3],[13,0.65],[14,0.82],
[15,0.70],[16,0.5],[17,0.11],[18,0.11],
[70,0.5],[72,1],[73,0.5],[80,1],[191,0.15],
[192,0.15],[241,0.1],[242,0.1]]))
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
Impcnt = float(row.getValue("Count"))
return Impcnt
def TcImpMRLC(inLU):
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[21,0.25],[22,0.65],[23,0.82],[31,1.0],
[32,0.11],[33,0.5],[85,0.11]]))
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
Impcnt = float(row.getValue("Count"))
return Impcnt
def TcImpUltimate(inLU):
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,0.25],[12,0.30],[13,0.65],[14,0.82],
[15,0.7],[16,0.5],[17,0.11],[18,0.11],
[70,0.5],[72,0.1],[73,0.5],[80,1],[191,0.15],
[192,0.15],[241,0.1],[242,0.1],[111,0.65],
[112,0.38],[113,0.3],[114,0.25],[115,0.2],[116,0.12]]))
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
Impcnt = float(row.getValue("Count"))
return Impcnt
def TcImpUSGS(inLU):
LU = arcpy.sa.Reclassify(inLU,"Value",
RemapValue([[11,0.3],[12,0.82],[13,0.7],[14,1],
[15,0.76],[16,0.11],[17,0.11],[70,0.5],
[71,1],[74,1],[75,0.11],[76,0.5],[77,0.5]]))
rows = arcpy.SearchCursor(LU,"","","Value;Count","")
for row in rows:
Impcnt = float(row.getValue("Count"))
return Impcnt
#*******************************************************************************************************
# Check gages in the vicinity of delineated watershed and return gagelist and gagelistthomas
#*******************************************************************************************************
##def CheckGages(usgsgages,mdgages,outletpoint,outletbuffer,gagefield): # outletbuffer and gagefiled are only
## arcpy.MakeFeatureLayer_management(mdgages, "mdlayer") # path name to save files
## usgsfield = "GAGE_ID"
## mdgageslayer = "GAGE_ID"
## arcpy.AddJoin_management("mdlayer",mdgageslayer,usgsgages,usgsfield)
## #arcpy.Buffer_analysis(outletpoint, outletbuffer, "300","FULL","FLAT","ALL","FID") # ArcInfo License
## arcpy.Buffer_analysis(outletpoint, outletbuffer, "300","#","#","ALL","FID") # ArcView License
## arcpy.Intersect_analysis(["mdlayer",outletbuffer],gagefield,"ALL","#","INPUT")
## gagelist = []
## gagevalue = arcpy.SearchCursor(gagefield,"","","usgsgag_14;usgsgag_16","")
## for isgood in gagevalue:
## if isgood.getValue("usgsgag_16") == "y":
## gage = isgood.getValue("usgsgag_14")
## gagelist.append(gage)
def CheckGages(usgsgages,mdgages,outletpoint,mask,opath):
arcpy.env.workspace = opath
arcpy.env.scratchWorkspace = opath + "/aux_folder"
# added on 10-13-2017: Only using intersect tool to identify gauges which overlap watershed and mdgagugedstreams2016 file
arcpy.Intersect_analysis([mdgages,mask], opath + "/mask_ints.shp","ALL","#","INPUT")
arcpy.JoinField_management(opath + "/mask_ints.shp","GAGE_ID",usgsgages,"GAGE_ID","GAGE_ID")
arcpy.Intersect_analysis([opath + "/mask_ints.shp",outletpoint], opath + "/gauge_outlet.shp","ALL","#","INPUT")
gagelist = []
gagevalue = arcpy.SearchCursor(opath + "/gauge_outlet.shp","","","GAGE_ID","")
for g in gagevalue:
gid = g.getValue("GAGE_ID")
gagelist.append(gid)
return gagelist
#*******************************************************************************************************
# Generate tasker string and return percentile flow values
#*******************************************************************************************************
def TaskerString(period,qlist):
v1 = qlist[0]
v1 = int(v1)
v2 = qlist[1]
v2 = int(v2)
v3 = qlist[2]
v3 = int(v3)
v4 = qlist[3]
v4 = int(v4)
v5 = qlist[4]
v5 = int(v5)
v6 = qlist[5]
v6 = int(v6)
v7 = qlist[6]
v7 = int(v7)
v8 = qlist[7]
v8 = int(v8)
datastring = str(period).rjust(7)
datastring = datastring + str(v1).rjust(7)
datastring = datastring + str(v2).rjust(9)
datastring = datastring + str(v3).rjust(9)
datastring = datastring + str(v4).rjust(9)
datastring = datastring + str(v5).rjust(9)
datastring = datastring + str(v6).rjust(9)
datastring = datastring + str(v7).rjust(9)
datastring = datastring + str(v8).rjust(9)
return datastring
#*******************************************************************************************************
# Rating Table -- width, ymax, area, perimeter, and discharge calculations
#*******************************************************************************************************
def RatingTable(ymax, x, y, n, slope):
minElev = min(y)
imin = y.index(minElev)
nrow = len(y)
ibeg = 0
iend = nrow - 1
for i in range(imin, 1, -1):
if y[i] < ymax and y[i-1] >= ymax:
ibeg = i - 1
x[ibeg] = (x[i] + ((ymax - y[i])/(y[i-1] - y[i]) * (x[i-1] - x[i])))
y[ibeg] = ymax
for i in range(imin, (nrow-2)):
if y[i] < ymax and y[i+1] >= ymax:
iend = i + 1
x[iend] = (x[i] + ((ymax - y[i])/(y[i+1] - y[i]) * (x[i+1] - x[i])))
y[iend] = ymax
x1 = []
y1 = []
for i in range(ibeg,iend):
x1.append(x)
y1.append(y)
nrow2 = len(y)
sum_list = 0
sum2_list = 0
if nrow2 > 4:
for i in range(1, (nrow2 - 4)):
dx = x[i+1] - x[i]
dh = y[i+1] - y[i]
havg = 2*ymax - dh*0.5
sum_list = sum_list + (dx*havg) # area
sum2_list = sum2_list + (math.sqrt(pow(dx,2) + pow(dh,2))) # perimeter
if nrow2 > 3:
dx1 = x[1] - x[0]
dx2 = x[nrow2 - 1] - x[nrow2 - 2]
dh1 = ymax - y[1]
dh2 = ymax - y[nrow2 - 2]
sum_list = sum_list + ((dx1*dh1)/2) + ((dx2*dh2)/2)
sum2_list = sum2_list + (math.sqrt(pow(dx1,2)*pow(dh1,2))) + (math.sqrt(pow(dx2,2)*pow(dh2,2)))
else:
dx1 = x[2] - x[1]
dx2 = x[1] - x[0]
dh = ymax - y[1]
sum_list = sum_list + ((dx1*dh)/2) + ((dx2*dh)/2)
sum2_list = sum2_list + (math.sqrt(pow(dx1,2)*pow(dh,2))) + (math.sqrt(pow(dx2,2)*pow(dh,2)))
area = sum_list
width = math.fabs(x[ibeg] - x[iend])
depth = ymax
perimeter = sum2_list
vel = (1.49/n)* (math.sqrt(slope))* (pow((area/perimeter),2/3))
discharge = area*vel
return width, ymax, area, perimeter, discharge
#*******************************************************************************************************
# Get average precipitation list with all values set to -999 except last four set to 0
#*******************************************************************************************************
def GetAvgPrecList(yearlist,durlist):
avgpreclist = []
for i in yearlist:
for j in durlist:
avgpreclist.append(-999)
for k in avgpreclist[32:36]:
avgpreclist.remove(k)
avgpreclist.append(0)
return avgpreclist
#*******************************************************************************************************
# Get average precipitation list with all values set to -999 except last four set to 0
#*******************************************************************************************************
def GetPrecipFrequency(pr_path,durlist,theyear,thecritdur,opath):
arcpy.env.workspace = opath
arcpy.env.scratchWorkspace = opath + "/aux_folder"
precdist = []
for i in range(0,11): # change range to number of prcp check boxes
thedur = durlist[i]
thefilename = pr_path+ "p" + theyear + "-" + thedur + "m"
basingrid = arcpy.sa.Times(opath + "/basingrid",thefilename)
precavg = arcpy.GetRasterProperties_management(basingrid, "MEAN")
precavg = float(precavg.getOutput(0))
theavg = precavg/1000 # precip raster values are in cm*1000 (convert to inches)
precdist.append(theavg)
if thedur == thecritdur:
thecritavg = theavg
return precdist, thecritavg
#*******************************************************************************************************
# Get year index and precipitation values spcfied by user in "Choose Storm Depths" in list form
#*******************************************************************************************************
def PrcpandYear(cb,prcp):
year_uSpecified = []
prec_uSpecified = []
for i, (a,b) in enumerate(zip(cb,prcp)):
if a.GetValue():
year_uSpecified.append(i)
prec_uSpecified.append(b)
return year_uSpecified, prec_uSpecified
#*******************************************************************************************************
# Adjust storm depths using Areal Correction formulae obtained from:
# http://www.gishydro.eng.umd.edu/HydroPanel/hydrology_panel_report_3rd_edition_final.pdf
#
# Variables and Formulae:
#
# a = Alpha [0.008245], b = Beta [0.558], p = Phi [0.01044]
# r = Rho [0.4], g = Gamma [0.005], k = Kappa [0.5169], A = area in sq. miles
#
# For 6 Hour:
# RF = 1 - a*((A)**b)
#
# For 12 Hour:
# RF = 1 - (a/2)*((A)**b) - ((p/2)*((A)**r))
#
# For 24 Hour:
# RF = 1 - p*((A)**r)
#
# For 48 Hour:
# RF = 1 - g*((A)**k)
#
#*******************************************************************************************************
def ArealRF(year_uSpecified,prec_uSpecified,Area):
RF_precip = []
for i, (y,p) in enumerate(zip(year_uSpecified,prec_uSpecified)):
durlist = ["05n","10n","15n","30n","01","02","03","06","12","24","48"]
thecritdur = durlist[(i)%4+7]
A = float(Area)
if thecritdur == "06":
RF_6hour = float(p)*(1 - (0.008245*(pow(A,0.558))))
RF_precip.append(RF_6hour)
if thecritdur == "12":
RF_12hour = float(p)*(1 - ((0.008245/2)*(pow(A,0.558))) - ((0.01044/2)*(pow(A,0.4))))
RF_precip.append(RF_12hour)
if thecritdur == "24":
RF_24hour = float(p)*(1 - (0.01044*(pow(A,0.4))))
RF_precip.append(RF_24hour)
if thecritdur == "48":
RF_48hour = float(p)*(1 - (0.005*(pow(A,0.5169))))
RF_precip.append(RF_48hour)
RF_precip = ["%0.2f" %x for x in RF_precip]
return RF_precip
#*******************************************************************************************************
# Count number of sub-lists in a list -- special case for multiple provinces
#*******************************************************************************************************
def FloodIntervalLists(lists):
count = 0
for item in lists:
if isinstance(item,list):
count += 1 + FloodIntervalLists(item)
else:
count += 1
return count
#*******************************************************************************************************
# Compute number of segments and area downstream
#*******************************************************************************************************
def AverageArea(upval,downval,arcid_global,opath):
value = []
area = []
theme_to_find = opath + "/vel_meth/LongPathSub" + str(arcid_global)
sc = arcpy.SearchCursor(theme_to_find,"","","VALUE;Drain_Area","")
for i in sc:
v = i.getValue("VALUE")
value.append(v)
a = i.getValue("Drain_Area")
area.append(a)
upindex = 0
downindex = 0
for idx,v in enumerate(value):
if value[v] == upval:
upindex = value[v]
if value[v] == downval:
downindex = value[v]
dsarea = 0
thecount = 0
areaconv = 0.000347492 # converts pixels to square miles: 900 square meters to square miles
for ix in xrange(upindex,downindex+1):
thearea = area[ix]
dsarea = dsarea + thearea
thecount = thecount + 1
avgarea = float((float(dsarea) / float(thecount)) * areaconv)
return avgarea, dsarea
#*******************************************************************************************************
# Using os.startfile and webbrowser.open in ArcGIS for Desktop
# Source: https://arcpy.wordpress.com/2013/10/25/using-os-startfile-and-webbrowser-open-in-arcgis-for-desktop/
#*******************************************************************************************************
# A decorator that will run its wrapped function in a new thread
def run_in_other_thread(function):
# functool.wraps will copy over the docstring and some other metadata
# from the original function
@functools.wraps(function)
def fn_(*args, **kwargs):
thread = threading.Thread(target=function, args=args, kwargs=kwargs)
thread.start()
thread.join()
return fn_
# Our new wrapped versions of os.startfile and webbrowser.open
startfile = run_in_other_thread(os.startfile)
openbrowser = run_in_other_thread(webbrowser.open)
#*******************************************************************************************************
# Following function does:
# 1) extracts lu code and lu description from lookup table text files
# 2) prepare a dictionary for matching each lu code to lu description, and
# 3) finally matched dictionary to lu list from area of interest land use to get
#*******************************************************************************************************
#+++++++++++++++++++++
# Two functions:
# i) with key and dict
# ii) with a new dictionary -- nlcd table has lu code values not
# in ascending order. Use 2nd function.
#
#+++++++++++++++++++++
# lu description list
##def lu_description(lut_path,lu_match):
## lu_code = []
## lu_name = []
## with open(lut_path, "r") as f:
## next(f)
## for line in f:
## lu_strip = line.split("\t")[0] # lu code
## lu_code.append(lu_strip)
## dc_strip = line.split("\t")[1] # lu description
## lu_name.append(dc_strip)
## #convert lu_code to integer
## lu_code = [int(i) for i in lu_code]
## d = dict(zip(lu_name, lu_code))
## pythonaddins.MessageBox(d,"dict","0")
## rev1 = {v:k for k,v in d.iteritems()}
## rev2 = [rev1.get(item,item) for item in lu_match]
##
## return rev2
def lu_description(lut_file,lu_match):
lu_code = []
lu_desc = []
with open(lut_file, "r") as f:
next(f)
for line in f:
lu_strip = line.split("\t")[0] # lu code
lu_code.append(lu_strip)
dc_strip = line.split("\t")[1] # lu description
lu_desc.append(dc_strip)
#convert lu_code to integer
lu_code = [int(i) for i in lu_code]
d_new = dict(zip(lu_code,lu_desc))
rev1 = {v:k for v,k in d_new.iteritems()}
rev2 = [rev1.get(item,item) for item in lu_match]
return rev2
#*******************************************************************************************************
# Compare Discharges: Carpenter (1980), Dillow (1996), Moglen et al (2006), and Moglen et al (2010)
#*******************************************************************************************************
# (a) Carpentar calculations
##############################
def CarpenterQ(theprov, regionfrac, DA, FC, HA, HD, p2yr, SLC, ST):
Qlist = []
QSlist = []
# These SE/mean values correspond to the gage database results (3/13/2006)
# These values are based on USGS method to estimate standard error - log-normal distribution
AS = [-999, -999, -999, 0.403, 0.432, 0.469, 0.535, 0.600, 0.671, -999, -999]
BS = [-999, -999, -999, 0.957, 0.923, 0.926, 0.956, 0.992, 1.039, -999, -999]
PS = [-999, -999, -999, 0.465, 0.496, 0.536, 0.601, 0.658, 0.722, -999, -999]
WS = [-999, -999, -999, 1.125, 1.030, 1.000, 0.995, 1.002, 1.079, -999, -999]
ES = [-999, -999, -999, 0.585, 0.604, 0.616, 0.645, 0.678, 0.722, -999, -999]
if theprov == "A" or theprov == "P" or theprov == "B":
Q2 = 142.0 * (DA**0.745) * ((FC + 10)**(-0.273)) * (p2yr**0.669)
Q5 = 120.0 * (DA**0.731) * ((FC + 10)**(-0.275)) * (p2yr**1.358)
Q10 = 106.0 * (DA**0.724) * ((FC + 10)**(-0.286)) * (p2yr**1.810)
Q25 = 90.1 * (DA**0.717) * ((FC + 10)**(-0.307)) * (p2yr**2.376)
Q50 = 78.5 * (DA**0.712) * ((FC + 10)**(-0.323)) * (p2yr**2.793)
Q100 = 66.6 * (DA**0.708) * ((FC + 10)**(-0.336)) * (p2yr**3.212)
elif theprov == "E":
Q2 = 28.6 * (DA**0.910) * (SLC**0.681) * ((ST + 10)**(-0.148)) * ((FC + 10)**(-0.647)) * ((HA + 10)**(-0.309)) * ((HD + 10)**(0.560))
Q5 = 119.0 * (DA**0.989) * (SLC**0.843) * ((ST + 10)**(-0.533)) * ((FC + 10)**(-0.731)) * ((HA + 10)**(-0.369)) * ((HD + 10)**(0.577))
Q10 = 306.0 * (DA**1.016) * (SLC**0.911) * ((ST + 10)**(-0.820)) * ((FC + 10)**(-0.804)) * ((HA + 10)**(-0.367)) * ((HD + 10)**(0.624))
Q25 = 936.0 * (DA**1.039) * (SLC**0.974) * ((ST + 10)**(-1.114)) * ((FC + 10)**(-0.868)) * ((HA + 10)**(-0.384)) * ((HD + 10)**(0.655))
Q50 = 2120.0 * (DA**1.051) * (SLC**1.009) * ((ST + 10)**(-1.321)) * ((FC + 10)**(-0.916)) * ((HA + 10)**(-0.396)) * ((HD + 10)**(0.676))
Q100 = 4800.0 * (DA**1.060) * (SLC**1.035) * ((ST + 10)**(-1.519)) * ((FC + 10)**(-0.963)) * ((HA + 10)**(-0.410)) * ((HD + 10)**(0.695))
elif theprov == "W":
Q2 = 55.1 * (DA ** 0.672)
Q5 = 112.0 * (DA ** 0.670)
Q10 = 172.0 * (DA ** 0.667)
Q25 = 280.0 * (DA ** 0.666)
Q50 = 394.0 * (DA ** 0.665)
Q100 = 548.0 * (DA ** 0.662)
else:
pythonaddins.MessageBox("Error in acquiring geomorphological area code","Province Code Error",0)
Qlist.append(-999) ; Qlist.append(-999) ; Qlist.append(-999) ; Qlist.append(Q2) ; Qlist.append(Q5) ; Qlist.append(Q10)
Qlist.append(Q25) ; Qlist.append(Q50) ; Qlist.append(Q100) ; Qlist.append(-999) ; Qlist.append(-999)
Qlist = ["-" if i == -999 else round(i*regionfrac,1) for i in Qlist]
# Hydro.CalcQS is implemented below
if theprov == "A":
QSin = AS
elif theprov == "B":
QSin = BS
elif theprov == "P":
QSin = PS
elif theprov == "W":
QSin = WS
elif theprov == "E":
QSin = ES
for i,j in zip(Qlist,QSin):
if i == "-":
QSlist.append("-")
else:
val = round(i*regionfrac * (1+j),1)
QSlist.append(val)
return Qlist, QSlist
# (b) Dillow calculations
# based on Dillow (1996) equations from Hydrology Panel 4th Edition, Appendix 8
################################################################################
def DillowQ(theprov, regionfrac, DA, FC, BR, LIME, RCN, ST):
Qlist = []
QSlist = []
# These SE/mean values correspond to the gage database results (3/13/2006)
# These values are based on USGS method to estimate standard error - log-normal distribution
AS = [-999, -999, -999, 0.201, 0.243, 0.296, 0.376, 0.442, 0.510, -999, 0.687]
BS = [-999, -999, -999, 0.613, 0.617, 0.631, 0.661, 0.694, 0.737, -999, 0.864]
PS = [-999, -999, -999, 0.376, 0.359, 0.377, 0.431, 0.492, 0.567, -999, 0.769]
WS = [-999, -999, -999, 0.640, 0.702, 0.765, 0.873, 0.921, 1.039, -999, 1.290]
ES = [-999, -999, -999, 0.773, 0.788, 0.796, 0.792, 0.840, 0.855, -999, 0.897]
if theprov == "A":
Q2 = 106.0 * (DA**0.851) * ((FC + 10)**(-0.223)) * (BR**0.056)
Q5 = 109.0 * (DA**0.858) * ((FC + 10)**(-0.143)) * (BR**0.064)
Q10 = 113.0 * (DA**0.859) * ((FC + 10)**(-0.106)) * (BR**0.072)
Q25 = 118.0 * (DA**0.858) * ((FC + 10)**(-0.072)) * (BR**0.087)
Q50 = 121.0 * (DA**0.858) * ((FC + 10)**(-0.051)) * (BR**0.099)
Q100 = 124.0 * (DA**0.858) * ((FC + 10)**(-0.033)) * (BR**0.111)
Q500 = 127.0 * (DA**0.859) * ((FC + 10)**(-0.004)) * (BR**0.140)
elif theprov == "B":
Q2 = 4260.0 * (DA**0.774) * ((LIME + 10)**(-0.549)) * (BR**0.405)
Q5 = 6670.0 * (DA**0.752) * ((LIME + 10)**(-0.564)) * (BR**0.354)
Q10 = 8740.0 * (DA**0.741) * ((LIME + 10)**(-0.579)) * (BR**0.326)
Q25 = 12000.0 * (DA**0.730) * ((LIME + 10)**(-0.602)) * (BR**0.295)
Q50 = 15100.0 * (DA**0.723) * ((LIME + 10)**(-0.620)) * (BR**0.276)
Q100 = 18900.0 * (DA**0.719) * ((LIME + 10)**(-0.639)) * (BR**0.261)
Q500 = 31800.0 * (DA**0.712) * ((LIME + 10)**(-0.686)) * (BR**0.241)
elif theprov == "P":
Q2 = 451.0 * (DA**0.635) * ((FC + 10)**(-0.266))
Q5 = 839.0 * (DA**0.606) * ((FC + 10)**(-0.248))
Q10 = 1210.0 * (DA**0.589) * ((FC + 10)**(-0.242))
Q25 = 1820.0 * (DA**0.574) * ((FC + 10)**(-0.239))
Q50 = 2390.0 * (DA**0.565) * ((FC + 10)**(-0.240))
Q100 = 3060.0 * (DA**0.557) * ((FC + 10)**(-0.241))
Q500 = 5190.0 * (DA**0.543) * ((FC + 10)**(-0.245))
elif theprov == "W":
Q2 = 1410.0 * (DA**0.761) * ((FC + 10)**(-0.782))
Q5 = 1780.0 * (DA**0.769) * ((FC + 10)**(-0.687))
Q10 = 1910.0 * (DA**0.771) * ((FC + 10)**(-0.613))
Q25 = 2000.0 * (DA**0.772) * ((FC + 10)**(-0.519))
Q50 = 2060.0 * (DA**0.771) * ((FC + 10)**(-0.452))
Q100 = 2140.0 * (DA**0.770) * ((FC + 10)**(-0.391))
Q500 = 2380.0 * (DA**0.765) * ((FC + 10)**(-0.263))
elif theprov == "E":
Q2 = 0.25 * (DA**0.591) * ((RCN - 33)**(1.70)) * (BR**0.310) * ((FC + 10)**(-0.464)) * ((ST + 10)**(-0.148))
Q5 = 1.05 * (DA**0.595) * ((RCN - 33)**(1.74)) * (BR**0.404) * ((FC + 10)**(-0.586)) * ((ST + 10)**(-0.498))
Q10 = 3.24 * (DA**0.597) * ((RCN - 33)**(1.71)) * (BR**0.436) * ((FC + 10)**(-0.667)) * ((ST + 10)**(-0.694))
Q25 = 13.1 * (DA**0.597) * ((RCN - 33)**(1.66)) * (BR**0.457) * ((FC + 10)**(-0.770)) * ((ST + 10)**(-0.892))
Q50 = 35.0 * (DA**0.594) * ((RCN - 33)**(1.62)) * (BR**0.465) * ((FC + 10)**(-0.847)) * ((ST + 10)**(-1.01))
Q100 = 87.6 * (DA**0.589) * ((RCN - 33)**(1.58)) * (BR**0.470) * ((FC + 10)**(-0.923)) * ((ST + 10)**(-1.11))
Q500 = 627.0 * (DA**0.573) * ((RCN - 33)**(1.49)) * (BR**0.478) * ((FC + 10)**(-1.10)) * ((ST + 10)**(-1.29))
else:
pythonaddins.MessageBox("Error in acquiring geomorphological area code","Province Code Error",0)
Qlist.append(-999) ; Qlist.append(-999) ; Qlist.append(-999) ; Qlist.append(Q2) ; Qlist.append(Q5) ; Qlist.append(Q10)
Qlist.append(Q25) ; Qlist.append(Q50) ; Qlist.append(Q100) ; Qlist.append(-999) ; Qlist.append(Q500)
Qlist = ["-" if i == -999 else round(i*regionfrac,1) for i in Qlist]
# Hydro.CalcQS is implemented below
if theprov == "A":
QSin = AS
elif theprov == "B":
QSin = BS
elif theprov == "P":
QSin = PS
elif theprov == "W":
QSin = WS
elif theprov == "E":
QSin = ES
for i,j in zip(Qlist,QSin):
if i == "-":
QSlist.append("-")
else:
val = round(i*regionfrac * (1+j),1)
QSlist.append(val)
return Qlist, QSlist
# (c) Fixed Region calculations
# based on Moglen et al. (2006) equations from Hydrology Panel 4th Edition,
# Appendix 8 using Old Lime shapefile
################################################################################
def FixedRegionQ2006(theprov, regionfrac, DA, SLL, LIME, FC, IA, HD, BR, HA):
Qlist = []
QSlist = []
# These SE/mean values correspond to the gage database results (3/13/2006)
# These values are based on USGS method to estimate standard error - log-normal distribution
AS = [0.236, 0.219, 0.212, 0.207, 0.216, 0.242, 0.291, 0.331, 0.374, 0.418, 0.480]
BS = [0.746, 0.671, 0.652, 0.640, 0.554, 0.525, 0.516, 0.525, 0.544, 0.574, 0.623]
PSrural = [0.390, 0.338, 0.321, 0.313, 0.256, 0.243, 0.253, 0.275, 0.306, 0.342, 0.397]
PSurban = [0.417, 0.369, 0.356, 0.351, 0.285, 0.262, 0.260, 0.277, 0.307, 0.348, 0.412]
WS = [0.389, 0.363, 0.356, 0.354, 0.363, 0.406, 0.489, 0.547, 0.657, 0.755, 0.898]
ES = [0.342, 0.337, 0.342, 0.349, 0.369, 0.382, 0.400, 0.417, 0.440, 0.465, 0.508]
if theprov == "A":
Q1p25 = 70.25 * (DA ** 0.837) * (SLL ** 0.327)
Q1p50 = 87.42 * (DA ** 0.837) * (SLL ** 0.321)
Q1p75 = 96.37 * (DA ** 0.836) * (SLL ** 0.307)
Q2 = 101.41 * (DA ** 0.834) * (SLL ** 0.300)
Q5 = 179.13 * (DA ** 0.826) * (SLL ** 0.314)
Q10 = 255.75 * (DA ** 0.821) * (SLL ** 0.340)
Q25 = 404.22 * (DA ** 0.812) * (SLL ** 0.393)
Q50 = 559.80 * (DA ** 0.806) * (SLL ** 0.435)
Q100 = 766.28 * (DA ** 0.799) * (SLL ** 0.478)
Q200 = 1046.9 * (DA ** 0.793) * (SLL ** 0.525)
Q500 = 1565.0 * (DA ** 0.784) * (SLL ** 0.589)
elif theprov == "B":
Q1p25 = 57.39 * (DA ** 0.784) * ((LIME + 1) ** (-0.190))
Q1p50 = 81.45 * (DA ** 0.764) * ((LIME + 1) ** (-0.193))
Q1p75 = 96.33 * (DA ** 0.755) * ((LIME + 1) ** (-0.194))
Q2 = 107.20 * (DA ** 0.750) * ((LIME + 1) ** (-0.194))
Q5 = 221.28 * (DA ** 0.710) * ((LIME + 1) ** (-0.202))
Q10 = 336.84 * (DA ** 0.687) * ((LIME + 1) ** (-0.207))
Q25 = 545.62 * (DA ** 0.660) * ((LIME + 1) ** (-0.214))
Q50 = 759.45 * (DA ** 0.641) * ((LIME + 1) ** (-0.219))
Q100 = 1034.7 * (DA ** 0.624) * ((LIME + 1) ** (-0.224))
Q200 = 1387.6 * (DA ** 0.608) * ((LIME + 1) ** (-0.229))
Q500 = 2008.6 * (DA ** 0.587) * ((LIME + 1) ** (-0.235))
elif theprov == "P":
if IA > 10.0:
Q1p25 = 17.85 * (DA ** 0.652) * ((IA + 1.0) ** 0.635)
Q1p50 = 24.66 * (DA ** 0.648) * ((IA + 1.0) ** 0.631)
Q1p75 = 30.82 * (DA ** 0.643) * ((IA + 1.0) ** 0.611)
Q2 = 37.01 * (DA ** 0.635) * ((IA + 1.0) ** 0.588)
Q5 = 94.76 * (DA ** 0.624) * ((IA + 1.0) ** 0.499)
Q10 = 169.2 * (DA ** 0.622) * ((IA + 1.0) ** 0.435)
Q25 = 341.0 * (DA ** 0.619) * ((IA + 1.0) ** 0.349)
Q50 = 562.4 * (DA ** 0.619) * ((IA + 1.0) ** 0.284)
Q100 = 898.3 * (DA ** 0.619) * ((IA + 1.0) ** 0.222)
Q200 = 1413 * (DA ** 0.621) * ((IA + 1.0) ** 0.160)
Q500 = 2529 * (DA ** 0.623) * ((IA + 1.0) ** 0.079)
else:
Q1p25 = 202.9 * (DA ** 0.682) * ((FC + 1.0) ** (-0.222))
Q1p50 = 262.0 * (DA ** 0.683) * ((FC + 1.0) ** (-0.217))
Q1p75 = 308.9 * (DA ** 0.679) * ((FC + 1.0) ** (-0.219))
Q2 = 349.0 * (DA ** 0.674) * ((FC + 1.0) ** (-0.224))
Q5 = 673.8 * (DA ** 0.659) * ((FC + 1.0) ** (-0.228))
Q10 = 992.6 * (DA ** 0.649) * ((FC + 1.0) ** (-0.230))
Q25 = 1556 * (DA ** 0.635) * ((FC + 1.0) ** (-0.231))
Q50 = 2146 * (DA ** 0.624) * ((FC + 1.0) ** (-0.235))
Q100 = 2897 * (DA ** 0.613) * ((FC + 1.0) ** (-0.238))
Q200 = 3847 * (DA ** 0.603) * ((FC + 1.0) ** (-0.239))
Q500 = 5519 * (DA ** 0.589) * ((FC + 1.0) ** (-0.242))
elif theprov == "W":
Q1p25 = 18.62 * (DA ** 0.611) * ((IA + 1) ** 0.419) * ((HD + 1) ** 0.165)
Q1p50 = 21.97 * (DA ** 0.612) * ((IA + 1) ** 0.399) * ((HD + 1) ** 0.226)
Q1p75 = 24.42 * (DA ** 0.612) * ((IA + 1) ** 0.391) * ((HD + 1) ** 0.246)
Q2 = 26.32 * (DA ** 0.612) * ((IA + 1) ** 0.386) * ((HD + 1) ** 0.256)
Q5 = 42.64 * (DA ** 0.607) * ((IA + 1) ** 0.347) * ((HD + 1) ** 0.340)
Q10 = 58.04 * (DA ** 0.603) * ((IA + 1) ** 0.323) * ((HD + 1) ** 0.382)
Q25 = 86.25 * (DA ** 0.582) * ((IA + 1) ** 0.295) * ((HD + 1) ** 0.421)
Q50 = 111.50 * (DA ** 0.584) * ((IA + 1) ** 0.270) * ((HD + 1) ** 0.457)
Q100 = 143.56 * (DA ** 0.586) * ((IA + 1) ** 0.260) * ((HD + 1) ** 0.469)
Q200 = 185.15 * (DA ** 0.580) * ((IA + 1) ** 0.243) * ((HD + 1) ** 0.488)
Q500 = 256.02 * (DA ** 0.573) * ((IA + 1) ** 0.222) * ((HD + 1) ** 0.510)
elif theprov == "E":
Q1p25 = 19.85 * (DA ** 0.796) * (BR ** 0.066) * ((HA + 1) ** (-0.106))
Q1p50 = 20.48 * (DA ** 0.795) * (BR ** 0.156) * ((HA + 1) ** (-0.140))
Q1p75 = 20.81 * (DA ** 0.799) * (BR ** 0.197) * ((HA + 1) ** (-0.146))
Q2 = 20.95 * (DA ** 0.803) * (BR ** 0.222) * ((HA + 1) ** (-0.144))
Q5 = 25.82 * (DA ** 0.793) * (BR ** 0.368) * ((HA + 1) ** (-0.190))
Q10 = 31.17 * (DA ** 0.777) * (BR ** 0.439) * ((HA + 1) ** (-0.215))
Q25 = 40.26 * (DA ** 0.751) * (BR ** 0.511) * ((HA + 1) ** (-0.242))
Q50 = 50.00 * (DA ** 0.732) * (BR ** 0.549) * ((HA + 1) ** (-0.261))
Q100 = 63.44 * (DA ** 0.711) * (BR ** 0.576) * ((HA + 1) ** (-0.279))
Q200 = 79.81 * (DA ** 0.689) * (BR ** 0.601) * ((HA + 1) ** (-0.296))
Q500 = 108.7 * (DA ** 0.660) * (BR ** 0.628) * ((HA + 1) ** (-0.316))
else:
pythonaddins.MessageBox("Error in acquiring geomorphological area code","Province Code Error",0)
Qlist.append(Q1p25) ; Qlist.append(Q1p50) ; Qlist.append(Q1p75) ; Qlist.append(Q2) ; Qlist.append(Q5) ; Qlist.append(Q10)
Qlist.append(Q25) ; Qlist.append(Q50) ; Qlist.append(Q100) ; Qlist.append(Q200) ; Qlist.append(Q500)
Qlist = ["-" if i == -999 else round(i*regionfrac,1) for i in Qlist]
# Hydro.CalcQS is implemented below
if theprov == "A":
QSin = AS
elif theprov == "B":
QSin = BS
elif theprov == "P":
if IA > 10.0:
QSin = PSurban
else:
QSin = PSrural
elif theprov == "W":
QSin = WS
elif theprov == "E":
QSin = ES
for i,j in zip(Qlist,QSin):
if i == "-":
QSlist.append("-")
else:
val = round(i*regionfrac * (1+j),1)
QSlist.append(val)
return Qlist, QSlist
# (d) Fixed Region calculations
# based on Moglen et al. (2010) equations from Hydrology Panel 4th Edition,
# Appendix 8 using New Lime shapefile
################################################################################
def FixedRegionQ2010(theprov, regionfrac, DA, LIME, FC, IA, HCD, HA, SLL):
Qlist = []
QSlist = []
# These SE/mean values correspond to the gage database results (3/13/2006)
# These values are based on USGS method to estimate standard error - log-normal distribution
AS = [0.236, 0.219, 0.212, 0.207, 0.216, 0.242, 0.291, 0.331, 0.374, 0.418, 0.480]
BS = [0.421, 0.376, 0.652, 0.356, 0.314, 0.309, 0.322, 0.345, 0.375, 0.410, 0.463]
PSrural = [0.421, 0.376, 0.321, 0.356, 0.314, 0.309, 0.322, 0.345, 0.375, 0.410, 0.463]
PSurban = [0.417, 0.369, 0.356, 0.351, 0.285, 0.262, 0.260, 0.277, 0.307, 0.348, 0.412]
WS = [0.390, 0.364, 0.356, 0.332, 0.382, 0.427, 0.481, 0.540, 0.612, 0.696, 0.825]
ES = [0.324, 0.323, 0.342, 0.328, 0.351, 0.367, 0.393, 0.416, 0.442, 0.472, 0.516]
if theprov == "A":
Q1p25 = 70.25 * (DA ** 0.837) * (SLL ** 0.327)
Q1p50 = 87.42 * (DA ** 0.837) * (SLL ** 0.321)
Q2 = 101.41 * (DA ** 0.834) * (SLL ** 0.300)
Q5 = 179.13 * (DA ** 0.826) * (SLL ** 0.314)
Q10 = 255.75 * (DA ** 0.821) * (SLL ** 0.340)
Q25 = 404.22 * (DA ** 0.812) * (SLL ** 0.393)
Q50 = 559.80 * (DA ** 0.806) * (SLL ** 0.435)
Q100 = 766.28 * (DA ** 0.799) * (SLL ** 0.478)
Q200 = 1046.9 * (DA ** 0.793) * (SLL ** 0.525)
Q500 = 1565.0 * (DA ** 0.784) * (SLL ** 0.589)
elif theprov == "B":
Q1p25 = 287.1 * (DA ** 0.774) * ((LIME + 1) ** (-0.118)) * ((FC + 1) ** (-0.418))
Q1p50 = 327.3 * (DA ** 0.758) * ((LIME + 1) ** (-0.121)) * ((FC + 1) ** (-0.358))
Q2 = 396.9 * (DA ** 0.743) * ((LIME + 1) ** (-0.124)) * ((FC + 1) ** (-0.332))
Q5 = 592.5 * (DA ** 0.705) * ((LIME + 1) ** (-0.133)) * ((FC + 1) ** (-0.237))
Q10 = 751.1 * (DA ** 0.682) * ((LIME + 1) ** (-0.138)) * ((FC + 1) ** (-0.183))
Q25 = 996.0 * (DA ** 0.655) * ((LIME + 1) ** (-0.145)) * ((FC + 1) ** (-0.122))
Q50 = 1218.8 * (DA ** 0.635) * ((LIME + 1) ** (-0.150)) * ((FC + 1) ** (-0.082))
Q100 = 1471.1 * (DA ** 0.617) * ((LIME + 1) ** (-0.154)) * ((FC + 1) ** (-0.045))
Q200 = 1760.7 * (DA ** 0.600) * ((LIME + 1) ** (-0.159)) * ((FC + 1) ** (-0.009))
Q500 = 2215.4 * (DA ** 0.577) * ((LIME + 1) ** (-0.165)) * ((FC + 1) ** (-0.035))
elif theprov == "P":
if IA > 10.0:
Q1p25 = 17.85 * (DA ** 0.652) * ((IA + 1.0) ** 0.635)
Q1p50 = 24.66 * (DA ** 0.648) * ((IA + 1.0) ** 0.631)
Q2 = 37.01 * (DA ** 0.635) * ((IA + 1.0) ** 0.588)
Q5 = 94.76 * (DA ** 0.624) * ((IA + 1.0) ** 0.499)
Q10 = 169.2 * (DA ** 0.622) * ((IA + 1.0) ** 0.435)
Q25 = 341.0 * (DA ** 0.619) * ((IA + 1.0) ** 0.349)
Q50 = 562.4 * (DA ** 0.619) * ((IA + 1.0) ** 0.284)
Q100 = 898.3 * (DA ** 0.619) * ((IA + 1.0) ** 0.222)
Q200 = 1413.0 * (DA ** 0.621) * ((IA + 1.0) ** 0.160)
Q500 = 2529.0 * (DA ** 0.623) * ((IA + 1.0) ** 0.079)
else:
Q1p25 = 287.1 * (DA ** 0.774) * ((LIME + 1) ** (-0.118)) * ((FC + 1) ** (-0.418))
Q1p50 = 327.3 * (DA ** 0.758) * ((LIME + 1) ** (-0.121)) * ((FC + 1) ** (-0.358))
Q2 = 396.9 * (DA ** 0.743) * ((LIME + 1) ** (-0.124)) * ((FC + 1) ** (-0.332))
Q5 = 592.5 * (DA ** 0.705) * ((LIME + 1) ** (-0.133)) * ((FC + 1) ** (-0.237))
Q10 = 751.1 * (DA ** 0.682) * ((LIME + 1) ** (-0.138)) * ((FC + 1) ** (-0.183))
Q25 = 996.0 * (DA ** 0.655) * ((LIME + 1) ** (-0.145)) * ((FC + 1) ** (-0.122))
Q50 = 1218.8 * (DA ** 0.635) * ((LIME + 1) ** (-0.150)) * ((FC + 1) ** (-0.082))
Q100 = 1471.1 * (DA ** 0.617) * ((LIME + 1) ** (-0.154)) * ((FC + 1) ** (-0.045))
Q200 = 1760.7 * (DA ** 0.600) * ((LIME + 1) ** (-0.159)) * ((FC + 1) ** (-0.009))
Q500 = 2215.4 * (DA ** 0.577) * ((LIME + 1) ** (-0.165)) * ((FC + 1) ** (-0.035))
elif theprov == "W":
Q1p25 = 5.18 * (DA ** 0.694) * ((IA + 1) ** 0.382) * ((HCD + 1) ** 0.414)
Q1p50 = 6.73 * (DA ** 0.682) * ((IA + 1) ** 0.374) * ((HCD + 1) ** 0.429)
Q2 = 7.61 * (DA ** 0.678) * ((IA + 1) ** 0.362) * ((HCD + 1) ** 0.475)
Q5 = 10.5 * (DA ** 0.665) * ((IA + 1) ** 0.290) * ((HCD + 1) ** 0.612)
Q10 = 13.1 * (DA ** 0.653) * ((IA + 1) ** 0.270) * ((HCD + 1) ** 0.669)
Q25 = 17.5 * (DA ** 0.634) * ((IA + 1) ** 0.264) * ((HCD + 1) ** 0.719)
Q50 = 21.2 * (DA ** 0.621) * ((IA + 1) ** 0.263) * ((HCD + 1) ** 0.751)
Q100 = 25.6 * (DA ** 0.608) * ((IA + 1) ** 0.262) * ((HCD + 1) ** 0.781)
Q200 = 30.5 * (DA ** 0.596) * ((IA + 1) ** 0.261) * ((HCD + 1) ** 0.812)
Q500 = 37.9 * (DA ** 0.579) * ((IA + 1) ** 0.261) * ((HCD + 1) ** 0.849)
elif theprov == "E":
Q1p25 = 24.44 * (DA ** 0.815) * ((HA + 1) **(-0.139)) * (SLL ** 0.115)
Q1p50 = 32.14 * (DA ** 0.824) * ((HA + 1) **(-0.144)) * (SLL ** 0.194)
Q2 = 42.48 * (DA ** 0.836) * ((HA + 1) **(-0.158)) * (SLL ** 0.249)
Q5 = 81.20 * (DA ** 0.847) * ((HA + 1) **(-0.184)) * (SLL ** 0.385)
Q10 = 119.3 * (DA ** 0.844) * ((HA + 1) **(-0.196)) * (SLL ** 0.445)
Q25 = 186.7 * (DA ** 0.834) * ((HA + 1) **(-0.212)) * (SLL ** 0.499)
Q50 = 254.7 * (DA ** 0.824) * ((HA + 1) **(-0.222)) * (SLL ** 0.531)
Q100 = 340.4 * (DA ** 0.812) * ((HA + 1) **(-0.230)) * (SLL ** 0.557)
Q200 = 450.5 * (DA ** 0.800) * ((HA + 1) **(-0.237)) * (SLL ** 0.582)
Q500 = 638.7 * (DA ** 0.783) * ((HA + 1) **(-0.247)) * (SLL ** 0.610)
else:
pythonaddins.MessageBox("Error in acquiring geomorphological area code","Province Code Error",0)
Qlist.append(Q1p25) ; Qlist.append(Q1p50) ; Qlist.append(-999) ; Qlist.append(Q2) ; Qlist.append(Q5) ; Qlist.append(Q10)
Qlist.append(Q25) ; Qlist.append(Q50) ; Qlist.append(Q100) ; Qlist.append(Q200) ; Qlist.append(Q500)
Qlist = ["-" if i == -999 else round(i*regionfrac,1) for i in Qlist]
# Hydro.CalcQS is implemented below
if theprov == "A":
QSin = AS
elif theprov == "B":
QSin = BS
elif theprov == "P":
if IA > 10.0:
QSin = PSurban
else:
QSin = PSrural
elif theprov == "W":
QSin = WS
elif theprov == "E":
QSin = ES
for i,j in zip(Qlist,QSin):
if i == "-":
QSlist.append("-")
else:
val = round(i*regionfrac * (1+j),1)
QSlist.append(val)
return Qlist, QSlist
#*******************************************************************************************************
# Smooth transect plot
#*******************************************************************************************************
def smoothListGaussian(lst,strippedXs=False,degree=3):
window=degree*2-1
weight=np.array([1.0]*window)
weightGauss=[]
for i in range(window):
i=i-degree+1
frac=i/float(window)
gauss=1/(np.exp((4*(frac))**2))
weightGauss.append(gauss)
weight=np.array(weightGauss)*weight
smoothed=[0.0]*(len(lst)-window)
for i in range(len(smoothed)):
smoothed[i]=sum(np.array(lst[i:i+window])*weight)/sum(weight)
return smoothed
#*******************************************************************************************************
# Expand width and height of polygon for Will Thomas Tc calculation
# Original author:
# Tim Loesch
# GIS Program Manager
# MN Department of Natural Resources
# 651-259-5475
# tim.loesch@dnr.state.mn.us
#
# Modified for use in GISHydroNXT.
#*******************************************************************************************************
def expandExtent(theExtent,expPercent):
expandWidth = theExtent.width * expPercent
expandHeight = theExtent.height * expPercent
xMin = theExtent.XMin - expandWidth
xMax = theExtent.XMax + expandWidth
yMin = theExtent.YMin - expandHeight
yMax = theExtent.YMax + expandHeight
return [xMin,yMin,xMax,yMax]