# gSSURGO_ValuTable.py
#
# Steve Peaslee, National Soil Survey Center
# 2014-09-27
# Adapted from HorizonData_Mapping
# Calculates root zone depth, root zone available water supply and a set of
# available water supply values at various defined depth ranges.
#
# The root zone available water supply is calculated down to the root zone for
# each component and then a weighted average is calculated for the map unit.
#
# The rest of the AWS columns are calculated for all components at a variety of depths.
# Output is to a geodatabase table (Valu2)
#
# To do list for FY2016....
# In order for this script to meet the 2015 standard for the Valu1 table it needs to:
# skip components where compkind is null
# stop horizon calculations for all data where there is a discrepancy in horizon depths
# null map unit AWS, SOC and PWSL values where sum of all components > 100
# null Rootzone Depth, Rootzone AWS and NCCPI* where sum of major components > 100
# perform independent validation of SOC results
# compare metadata dates to make sure they are properly coded and inserted
# look at array processing for calculations based upon depth ranges
## ===================================================================================
class MyError(Exception):
pass
## ===================================================================================
def errorMsg():
try:
tb = sys.exc_info()[2]
tbinfo = traceback.format_tb(tb)[0]
theMsg = tbinfo + " \n" + str(sys.exc_type)+ ": " + str(sys.exc_value) + " \n"
PrintMsg(theMsg, 2)
except:
PrintMsg("Unhandled error in errorMsg method", 2)
pass
## ===================================================================================
def PrintMsg(msg, severity=0):
# Adds tool message to the geoprocessor
#
#Split the message on \n first, so that if it's multiple lines, a GPMessage will be added for each line
try:
for string in msg.split('\n'):
#Add a geoprocessing message (in case this is run as a tool)
if severity == 0:
arcpy.AddMessage(string)
elif severity == 1:
arcpy.AddWarning(string)
elif severity == 2:
arcpy.AddMessage(" ")
arcpy.AddError(string)
except:
pass
## ===================================================================================
def Number_Format(num, places=0, bCommas=True):
try:
# Format a number according to locality and given places
#locale.setlocale(locale.LC_ALL, "")
locale.getlocale()
if bCommas:
theNumber = locale.format("%.*f", (places, num), True)
else:
theNumber = locale.format("%.*f", (places, num), False)
return theNumber
except:
PrintMsg("Unhandled exception in Number_Format function (" + str(num) + ")", 2)
return False
## ===================================================================================
def GetLastDate(inputDB):
# Get the most recent date 'YYYYMMDD' from SACATALOG.SAVEREST and use it to populate metadata
#
try:
tbl = os.path.join(inputDB, "SACATALOG")
today = ""
with arcpy.da.SearchCursor(tbl, ['SAVEREST'], "", "", [None, "ORDER_BY SAVEREST DESC"]) as cur:
for rec in cur:
#lastDate = rec[0].split(" ")[0].replace("-", "")
lastDate = rec[0].strftime('%Y%m%d')
break
return lastDate
except MyError, e:
# Example: raise MyError("this is an error message")
PrintMsg(str(e) + " \n", 2)
return ""
except:
errorMsg()
return ""
## ===================================================================================
def CreateQueryTables(inputDB, outputDB, maxD):
#
try:
# Problems with query table not getting all records where join table record is null
# FoQueryTable_CRr instance, I appear to be missing a component record when there is no chorizon data,
# such as for 'Miscellaneous area'
env.workspace = inputDB
queryMU = "MU"
queryCO = "CO"
queryHZ = "HZ"
queryCR = "CR"
queryCT = "CT"
queryTemp = "Tmp" # empty table uses as template
# Create an empty table using the original queryHZ as a template
# Mapunit query
PrintMsg(" \nReading MAPUNIT table...", 0)
whereClause = ""
fldMu = [["mukey", "mukey"], ["musym", "musym"], ["muname", "muname"]]
#arcpy.MakeQueryTable_management(["mapunit"], queryMU, "USE_KEY_FIELDS", "#", fldMu, whereClause)
fldMu2 = list()
dMu = dict()
#sqlClause = (None, "ORDER BY cokey, resdept_r ASC")
sqlClause = (None, "ORDER BY mukey")
for fld in fldMu:
fldMu2.append(fld[0])
muList = list()
#with arcpy.da.SearchCursor(queryMU, fldMu2, "", "", "", sqlClause) as mcur:
muTbl = os.path.join(inputDB, "mapunit")
tmukey = 515455
with arcpy.da.SearchCursor(muTbl, fldMu2, "", "", "", sqlClause) as mcur:
for mrec in mcur:
rec = list(mrec)
mukey = int(rec[0])
rec.pop(0)
dMu[mukey] = rec
muList.append(mukey)
#if mukey == tmukey:
# PrintMsg(" \nFound mukey " + str(tmukey) + " in mapunit table", 1)
muList.sort()
# Component query
PrintMsg(" \nReading COMPONENT table...", 0)
fldCo = [["mukey", "mukey"], ["cokey", "cokey"], ["comppct_r", "comppct_r"], ["majcompflag", "majcompflag"], \
["compname", "compname"], ["compkind", "compkind"], ["taxorder", "taxorder"], ["taxsubgrp", "taxsubgrp"], \
["localphase", "localphase"], ["otherph", "otherph"], ["hydricrating", "hydricrating"], ["drainagecl", "drainagecl"]]
fldCo2 = list()
dCo = dict()
whereClause = "comppct_r is not NULL"
sqlClause = (None, "ORDER BY cokey, comppct_r DESC")
coTbl = os.path.join(inputDB, "component")
for fld in fldCo:
fldCo2.append(fld[0])
with arcpy.da.SearchCursor(coTbl, fldCo2, whereClause, "", "", sqlClause) as ccur:
for crec in ccur:
rec = list(crec)
mukey = int(rec.pop(0)) # get rid of mukey from component record
try:
# Add next component record to list
dCo[mukey].append(rec)
except:
# initialize list of records
dCo[mukey] = [rec]
# HORIZON TABLE
PrintMsg(" \nReading HORIZON table...", 0)
fldHz = [["cokey", "cokey"], ["chkey", "chkey"], ["hzname", "hzname"], ["desgnmaster", "desgnmaster"], \
["hzdept_r", "hzdept_r"], ["hzdepb_r", "hzdepb_r"], ["sandtotal_r", "sandtotal_r"], \
["silttotal_r", "silttotal_r"], ["claytotal_r", "claytotal_r"], ["om_r", "om_r"], \
["dbthirdbar_r", "dbthirdbar_r"], ["ec_r", "ec_r"], ["ph1to1h2o_r", "ph1to1h2o_r"], \
["awc_r", "awc_r"]]
fldHz2 = list()
dHz = dict()
whereClause = "hzdept_r is not NULL and hzdepb_r is not NULL"
sqlClause = (None, "ORDER BY chkey, hzdept_r ASC")
for fld in fldHz:
fldHz2.append(fld[0])
hzTbl = os.path.join(inputDB, "chorizon")
with arcpy.da.SearchCursor(hzTbl, fldHz2, whereClause, "", "", sqlClause) as hcur:
for hrec in hcur:
rec = list(hrec)
cokey = int(rec.pop(0))
try:
# Add next horizon record to list
dHz[cokey].append(rec)
except:
# initialize list of horizon records
dHz[cokey] = [rec]
# HORIZON TEXTURE QUERY
#
PrintMsg(" \nMaking query table (CT) for texture information", 0)
inputTbls = list()
tbls = ["chtexturegrp", "chtexture"]
for tbl in tbls:
inputTbls.append(os.path.join(inputDB, tbl))
txList1 = [["chtexturegrp.chkey", "chkey"], ["chtexturegrp.texture", "texture"], ["chtexture.lieutex", "lieutex"]]
whereClause = "chtexturegrp.chtgkey = chtexture.chtgkey and chtexturegrp.rvindicator = 'Yes'"
arcpy.MakeQueryTable_management(inputTbls, queryCT, "USE_KEY_FIELDS", "#", txList1, whereClause)
# Read texture query into dictionary
txList2 = ["chtexturegrp.chkey", "chtexturegrp.texture", "chtexture.lieutex"]
dTexture = dict()
ctCnt = int(arcpy.GetCount_management(queryCT).getOutput(0))
arcpy.SetProgressor ("step", "Getting horizon texture information for QueryTable_HZ...", 0, ctCnt, 1)
# Have been running out of memory here if other applications are running. 5.66GB
with arcpy.da.SearchCursor(queryCT, txList2) as cur:
for rec in cur:
dTexture[int(rec[0])] = [rec[1], rec[2]]
arcpy.SetProgressorPosition()
arcpy.Delete_management(queryCT)
# COMPONENT RESTRICTIONS which will be saved to a gdb table
#
crTbl = os.path.join(inputDB, "corestrictions")
PrintMsg(" \nWriting component restrictions to " + os.path.join(outputDB, "QueryTable_CR"), 0)
fldCr = [["cokey", "cokey"], \
["reskind", "reskind"],\
["reshard", "reshard"],\
["resdept_r", "resdept_r"]]
whereClause = "OBJECTID = 1"
arcpy.MakeQueryTable_management(crTbl, queryCR, "USE_KEY_FIELDS", "#", fldCr, whereClause)
arcpy.CreateTable_management(outputDB, "QueryTable_CR", queryCR)
arcpy.Delete_management(queryCR)
fldCr2 = list()
dCr = dict()
sqlClause = (None, "ORDER BY cokey, resdept_r ASC")
for fld in fldCr:
fldCr2.append(fld[0])
whereClause = "resdept_r is not NULL and resdept_r < " + str(maxD)
outputCR = os.path.join(outputDB, "QueryTable_CR")
crList = list() # list of components with a restriction
crCnt = int(arcpy.GetCount_management(crTbl).getOutput(0))
arcpy.SetProgressor ("step", "Saving component restriction information...", 0, crCnt, 1)
with arcpy.da.SearchCursor(crTbl, fldCr2, whereClause, "", "", sqlClause) as crcur:
for crrec in crcur:
rec = list(crrec)
cokey = int(rec[0])
arcpy.SetProgressorPosition()
if not cokey in crList:
# Only save the highest level restriction above 150cm
crList.append(cokey)
dCr[cokey] = rec
PrintMsg(" \nCreating QueryTable_HZ in " + outputDB, 0)
fldCo.pop(0)
fldCo2.pop(0)
fldHz.pop(0)
fldHz2.pop(0)
# Create list of fields for query table
fldAll = list()
# Create list of fields for output cursor
fldAll2 = list()
for fld in fldMu:
fldAll.append(["mapunit." + fld[0], fld[1]])
fldAll2.append(fld[1])
for fld in fldCo:
fldAll.append(["component." + fld[0], fld[1]])
fldAll2.append(fld[1])
for fld in fldHz:
fldAll.append(["chorizon." + fld[0], fld[1]])
fldAll2.append(fld[1])
# Texture fields:
fldAll2.append("texture")
fldAll2.append("lieutex")
# Select component-horizon data for ALL components that have horizon data. Lack of horizon data
# will cause some components to be missing from the PWSL.
#
# Later on in the actual calculations for RZAWS, only the major-earthy components will be used. But
# all components are in this table!
whereClause = "mapunit.mukey = component.mukey and \
component.cokey = chorizon.cokey and mapunit.objectid = 1"
outputTable = os.path.join(outputDB, "QueryTable_HZ")
PrintMsg(" \nCreating table " + outputTable, 0)
arcpy.MakeQueryTable_management(['mapunit', 'component', 'chorizon'], queryTemp, "USE_KEY_FIELDS", "#", fldAll, whereClause)
arcpy.CreateTable_management(outputDB, "QueryTable_HZ", queryTemp)
arcpy.AddField_management(outputTable, "texture", "TEXT", "", "", "30", "texture")
arcpy.AddField_management(outputTable, "lieutex", "TEXT", "", "", "254", "lieutex")
arcpy.Delete_management(queryTemp)
# Process dictionaries and use them to write out the new QueryTable_HZ table
#
# Open output table
outFld2 = arcpy.Describe(outputTable).fields
outFlds = list()
for fld in outFld2:
outFlds.append(fld.name)
outFlds.pop(0)
# Create empty lists to replace missing data
missingCo = ["", None, None, None, None, None, None, None, None, None, None]
missingHz = ["", None, None, None, None, None, None, None, None, None, None, None, None]
missingTx = [None, None]
# Save information on mapunits or components with bad or missing data
#badMu = list() # list of mapunits with no components
muNoCo = list()
dNoCo = dict()
coNoHz = list() # list of components with no horizons
dNoHz = dict() # component data for those components in coNoHz
arcpy.SetProgressor ("step", "Writing data to " + outputTable + "...", 0, len(muList), 1)
with arcpy.da.InsertCursor(outputTable, fldAll2) as ocur:
for mukey in muList:
mrec = dMu[mukey]
arcpy.SetProgressorPosition()
try:
coVals = dCo[mukey] # got component records for this mapunit
# Sort lists by comppct_r
coList = sorted(coVals, key = lambda x: int(x[1]))
for corec in coList:
cokey = int(corec[0])
#if mukey == tmukey:
# PrintMsg("\t\tCO1: " + str(cokey) + ": " + str(corec), 1)
try:
hzVals = dHz[cokey] # horizon records for this component
# Sort record by hzdept_r
hzList = sorted(hzVals, key = lambda x: int(x[3]))
for hzrec in hzList:
chkey = int(hzrec[0])
#if mukey == tmukey:
# PrintMsg("\t\t\tHZ1: " + str(chkey) + ": " + str(hzrec), 1)
try:
# Get horizon texture
txrec = dTexture[chkey]
except:
txrec = missingTx
# Combine all records and write to table
newrec = [mukey]
newrec.extend(mrec)
newrec.extend(corec)
newrec.extend(hzrec)
newrec.extend(txrec)
#if mukey == tmukey:
# PrintMsg("\t\t\tHZ2: " + str(cokey) + ": " + str(corec), 1)
ocur.insertRow(newrec)
except KeyError:
# No horizon records for this component
comppct = corec[1]
compname = corec[3]
compkind = corec[4]
mjrcomp = corec[2]
#PrintMsg("Major compflag = " + str(corec), 1)
hzrec = missingHz
txrec = missingTx
newrec = [mukey]
newrec.extend(mrec)
newrec.extend(corec)
newrec.extend(hzrec)
newrec.extend(txrec)
#if mukey == tmukey:
#PrintMsg("\t\tNo horizon data for mapunit:component: " + str(mukey) + ":" + str(cokey), 1)
#PrintMsg("\t\t\tHZ3: " + str(cokey) + ": " + str(newrec), 1)
ocur.insertRow(newrec)
if not (compname in ["NOTCOM", "NOTPUB"] or compkind == 'Miscellaneous area'):
badComp = [mukey, str(cokey), compname, compkind, mjrcomp, str(comppct)]
coNoHz.append(str(cokey)) # add cokey to list of components with no horizon data
dNoHz[cokey] = badComp # add component information to dictionary
#PrintMsg(" \nMissing horizon data: " + str(corec), 1)
except:
PrintMsg(" \nhzVals error for " + str(mukey) + ":" + str(cokey) + ": " + str(txrec), 2)
PrintMsg(" \n" + str(fldAll2), 1)
errorMsg()
except:
# No component records for this map unit
corec = missingCo
hzrec = missingHz
txrec = missingTx
newrec = [mukey]
newrec.extend(mrec)
newrec.extend(corec)
newrec.extend(hzrec)
newrec.extend(txrec)
#if mukey == tmukey:
#PrintMsg("\t\tNo component data for mapunit:component: " + str(mukey) + ":" + str(cokey), 1)
ocur.insertRow(newrec)
if not mrec[0] in ['NOTCOM', 'NOTPUB']:
# skip map units that should never have component data
#
muNoCo.append(str(mukey))
dNoCo[str(mukey)] = [mrec[0], mrec[1]] # Save map unit name for the report
#PrintMsg(" \n\n** No component data for " + str(mrec[1]), 2)
# Run through QueryTbl_HZ table, checking for inconsistencies in horizon depths
# Create a dictionary containing a list of top and bottom of each horizon in each component
# dictionary key = cokey
# list contains tuples of hzdept_r, hzdepb_r, hzname, mukey, compname, localphase
dHZ = dict()
# Exclude horizon data with null hzdep with whereclause
wc = "hzdept_r is not null and hzdepb_r is not null"
arcpy.ResetProgressor()
arcpy.SetProgressorLabel("Looking for inconsistencies in horizon depths...")
with arcpy.da.SearchCursor(outputTable, ['mukey', 'cokey','hzdept_r','hzdepb_r', 'hzname', 'compname', 'localphase', 'majcompflag'], wc) as cur:
for rec in cur:
mukey, cokey, top, bot, hzname, compname, localphase, majcomp = rec
cokey = int(cokey)
if cokey in dHZ:
vals = dHZ[cokey]
vals.append([top, bot, hzname, mukey, compname, localphase, majcomp])
dHZ[cokey] = vals
else:
dHZ[cokey] = [[top, bot, hzname, mukey, compname, localphase, majcomp]]
# Number of items in each component value = len(rec)
# top of first horizon = rec[0][0]
# bottom of last horizon = rec[len(rec) - 1][1]
# component thickness #1 = rec[len(rec) - 1][1] - rec[0][0]
# Read each entry in the dictionary and check for gaps and overlaps in the horizons
#
badCoHz = list()
badHorizons = list()
for cokey, vals in dHZ.items():
# Process each component
#
hzSum = 0 # itialize sum of horizon thicknesses
lb = vals[0][0] # initialize last bottom to the top of the first horizon
localphase = vals[0][5]
if localphase is None:
localphase = ""
else:
localphase = " " + localphase
for v in vals:
# Process each horizon in the component record
#
# sum of bottom - top for each horizon
hzSum += (v[1] - v[0])
# Check for consistency between tops and bottoms for each consecutive horizon
if v[0] != lb:
diff = v[0] - lb
badCoHz.append(str(cokey))
badHorizons.append(v[3] + ", " + str(cokey) + ", " + v[4] + localphase + ", " + majcomp + ", " + str(v[2]) + ", " + str(v[0]) + ", " + str(diff) )
lb = v[1] # update last bottom depth
PrintMsg(" \nWriting component restrictions to " + outputCR, 0)
arcpy.SetProgressor ("step", "Writing component restriction data to " + outputCR + "...", 0, len(dCr), 1)
with arcpy.da.InsertCursor(outputCR, fldCr2) as ocur:
for cokey, crrec in dCr.items():
ocur.insertRow(crrec)
# Save data issues to permanent files for later review
#
# Don't report these errors from this script tool. This code has been moved to another tool
if 1 == 2: # Skip this code. May use some of it to null out some valu table results.
#
#if len(muNoCo) > 0 or len(coNoHz) or len(badCoHz) > 0:
PrintMsg(" \nCreating log file: " + logFile, 1)
now = datetime.now()
fh = open(logFile, "w")
fh.write("\n" + inputDB + "\n")
fh.write("\nProcessed on " + now.strftime('%A %x %X') + "\n\n")
fh.write("This log file containins record of any data inconsistencies found while processing the Valu2 table \n\n")
fh.close()
# Save data issues (mapunits with no components) to log file for later review
# Some of these will be NOTCOMs
if len(muNoCo) > 0:
fh = open(logFile, "a")
#fh.write(inputDB + "\n")
fh.write("\nQuery for map units missing component data\n")
fh.write("====================================================================================\n")
fh.write("MUKEY IN ('" + "', '".join(muNoCo) + "') \n\n")
fh.write("\n\nTable of map units missing component data\n")
fh.write("\nMUKEY, MUSYM, MUNAME\n")
for mukey in muNoCo:
fh.write(str(mukey) + ", " + dNoCo[mukey][0] + ", " + dNoCo[mukey][1] + "\n")
fh.close()
PrintMsg(" \nMap units missing component data (" + Number_Format(len(muNoCo), 0, True) + ") saved to:\t" + logFile, 0)
# Save data issues (components with no horizons) to log file for later review
# Note; these COKEYs will work with gSSURGO but not Soil Data Access
# Some of these may be NOTCOMs
if len(coNoHz) > 0:
PrintMsg(" \nQuery for components missing horizon data (" + Number_Format(len(coNoHz), 0, True) + ") saved to:\t" + logFile, 0)
fh = open(logFile, "a")
fh.write("\n\nQuery for components with no horizon data\n")
fh.write("====================================================================================\n")
fh.write("COKEY IN ('" + "', '".join(coNoHz) + "') \n\n")
fh.write("Table of map-unit\components missing horizon data\n")
fh.write("\nMUKEY, COKEY, COMPNAME, COMPKIND, MAJCOMPFLAG, COMPPCT, COMPKIND\n")
for mukey, compInfo in dNoHz.items():
# mukey, str(cokey), compname, compkind, mjrcomp, str(comppct)
mukey, cokey, compname, compkind, majcomp, comppct = compInfo
#compInfo = mukey, str(cokey), compname, compkind, str(comppct)
fh.write(str(mukey) + ", " + str(cokey) + ", " + compname + ", " + str(compkind) + ", " + majcomp + ", " + str(comppct) + ", " + str(compkind) + "\n")
fh.close()
# These components have inconsistent horizon depths that overlap or gap.
# Note; these COKEYs will work with gSSURGO but not Soil Data Access
if len(badCoHz) > 0:
PrintMsg(" \nComponents with horizon gaps or overlaps (" + Number_Format(len(badCoHz), 0, True) + ") saved to:\t" + logFile, 0)
fh = open(logFile, "a")
fh.write("\n\nQuery for components with horizon gaps or overlaps\n")
fh.write("====================================================================================\n")
fh.write("COKEY IN ('" + "', '".join(badCoHz) + "') \n\n")
fh.write("\nTable of components with horizon gaps or overlaps\n")
fh.write("MUKEY, COKEY, COMPNAME, MJRCOMP, HZNAME, HZDEPT, DIFF\n")
for h in badHorizons:
fh.write(h + "\n")
fh.close()
arcpy.ResetProgressor()
env.workspace = outputDB
return True
except MyError, e:
# Example: raise MyError("this is an error message")
PrintMsg(str(e) + " \n", 2)
return False
except:
errorMsg()
return False
## ===================================================================================
def CreateOutputTableMu(theMuTable, depthList, dPct):
# Create the mapunit level table
#
try:
# Create the output tables and add required fields
try:
# Try to handle existing output table if user has added it to ArcMap from a previous run
if arcpy.Exists(theMuTable):
arcpy.Delete_management(theMuTable)
except:
raise MyError, "Previous output table (" + theMuTable + ") is in use and cannot be removed"
return False
PrintMsg(" \nCreating new output table (" + os.path.basename(theMuTable) + ") for mapunit level data", 0)
outputDB = os.path.dirname(theMuTable)
tmpTable = os.path.join("IN_MEMORY", os.path.basename(theMuTable))
#arcpy.CreateTable_management(outputDB, os.path.basename(theMuTable))
arcpy.CreateTable_management("IN_MEMORY", os.path.basename(theMuTable))
# Add fields for AWS
for rng in depthList:
# Create the AWS fields in a loop
#
td = rng[0]
bd = rng[1]
awsField = "AWS" + str(td) + "_" + str(bd)
arcpy.AddField_management(tmpTable, awsField, "SHORT", "", "", "", awsField)
for rng in depthList:
# Create the AWS fields in a loop
#
td = rng[0]
bd = rng[1]
awsField = "TK" + str(td) + "_" + str(bd) + "A"
arcpy.AddField_management(tmpTable, awsField, "FLOAT", "", "", "", awsField)
arcpy.AddField_management(tmpTable, "MUSUMCPCTA", "SHORT", "", "", "")
# Add Fields for SOC
for rng in depthList:
# Create the SOC fields in a loop
#
td = rng[0]
bd = rng[1]
socField = "SOC" + str(td) + "_" + str(bd)
arcpy.AddField_management(tmpTable, socField, "LONG", "", "", "", socField)
for rng in depthList:
# Create the SOC fields in a loop
#
td = rng[0]
bd = rng[1]
socField = "TK" + str(td) + "_" + str(bd) + "S"
arcpy.AddField_management(tmpTable, socField, "FLOAT", "", "", "", socField)
arcpy.AddField_management(tmpTable, "MUSUMCPCTS", "SHORT", "", "", "")
# Add fields for NCCPI
arcpy.AddField_management(tmpTable, "NCCPI2CS", "FLOAT", "", "", "")
arcpy.AddField_management(tmpTable, "NCCPI2SG", "FLOAT", "", "", "")
arcpy.AddField_management(tmpTable, "NCCPI2CO", "FLOAT", "", "", "")
arcpy.AddField_management(tmpTable, "NCCPI2ALL", "FLOAT", "", "", "")
# Add fields for root zone depth and root zone available water supply
arcpy.AddField_management(tmpTable, "PCTEARTHMC", "SHORT", "", "", "")
arcpy.AddField_management(tmpTable, "ROOTZNEMC", "SHORT", "", "", "")
arcpy.AddField_management(tmpTable, "ROOTZNAWS", "SHORT", "", "", "")
# Add field for droughty soils
arcpy.AddField_management(tmpTable, "DROUGHTY", "SHORT", "", "", "")
# Add field for potential wetland soils
arcpy.AddField_management(tmpTable, "PWSL1POMU", "SHORT", "", "", "")
# Add field for mapunit-sum of ALL component-comppct_r values
arcpy.AddField_management(tmpTable, "MUSUMCPCT", "SHORT", "", "", "")
# Add Mukey field (primary key)
arcpy.AddField_management(tmpTable, "MUKEY", "TEXT", "", "", "30", "MUKEY")
# Convert IN_MEMORY table to a permanent table
arcpy.CreateTable_management(outputDB, os.path.basename(theMuTable), tmpTable)
# Add attribute indexes for key fields
arcpy.AddIndex_management(theMuTable, "MUKEY", "Indx_ResMukey", "NON_UNIQUE", "NON_ASCENDING")
arcpy.Delete_management(os.path.join("IN_MEMORY", os.path.basename(theMuTable)))
# Reading from the original mapunit table, populate the output Valu2 table with mukey and musumpct
#PrintMsg(" \n\tPopulating " + theMuTable + " with mukey values", 1)
with arcpy.da.SearchCursor(os.path.join(outputDB, "mapunit"), ["mukey"]) as incur:
outcur = arcpy.da.InsertCursor(theMuTable, ["mukey", "musumcpct"])
for inrec in incur:
mukey = inrec[0]
try:
sumPct = dPct[mukey][0]
except:
sumPct = 0
inrec = [mukey, sumPct]
outcur.insertRow(inrec)
return True
except MyError, e:
# Example: raise MyError("this is an error message")
PrintMsg(str(e) + " \n", 2)
return False
except:
errorMsg()
return False
## ===================================================================================
def CreateOutputTableCo(theCompTable, depthList):
# Create the component level table
# The new input field is created using adaptive code from another script.
#
try:
# Create two output tables and add required fields
try:
# Try to handle existing output table if user has added it to ArcMap from a previous run
if arcpy.Exists(theCompTable):
arcpy.Delete_management(theCompTable)
except:
raise MyError, "Previous output table (" + theCompTable + ") is in use and cannot be removed"
return False
PrintMsg(" \nCreating new output table (" + os.path.basename(theCompTable) + ") for component level data", 0)
outputDB = os.path.dirname(theCompTable)
tmpTable = os.path.join("IN_MEMORY", os.path.basename(theCompTable))
arcpy.CreateTable_management("IN_MEMORY", os.path.basename(theCompTable))
# Add fields appropriate for the component level restrictions
# mukey,cokey, compName, localphase, compPct, comppct, resdept, restriction
arcpy.AddField_management(tmpTable, "COKEY", "TEXT", "", "", "30", "COKEY")
arcpy.AddField_management(tmpTable, "COMPNAME", "TEXT", "", "", "60", "COMPNAME")
arcpy.AddField_management(tmpTable, "LOCALPHASE", "TEXT", "", "", "40", "LOCALPHASE")
arcpy.AddField_management(tmpTable, "COMPPCT_R", "SHORT", "", "", "", "COMPPCT_R")
for rng in depthList:
# Create the AWS fields in a loop
#
td = rng[0]
bd = rng[1]
awsField = "AWS" + str(td) + "_" + str(bd)
arcpy.AddField_management(tmpTable, awsField, "FLOAT", "", "", "", awsField)
for rng in depthList:
# Create the AWS fields in a loop
#
td = rng[0]
bd = rng[1]
awsField = "TK" + str(td) + "_" + str(bd) + "A"
arcpy.AddField_management(tmpTable, awsField, "FLOAT", "", "", "", awsField)
arcpy.AddField_management(tmpTable, "MUSUMCPCTA", "SHORT", "", "", "")
for rng in depthList:
# Create the SOC fields in a loop
#
td = rng[0]
bd = rng[1]
awsField = "SOC" + str(td) + "_" + str(bd)
arcpy.AddField_management(tmpTable, awsField, "FLOAT", "", "", "")
for rng in depthList:
# Create the rest of the SOC thickness fields in a loop
#
td = rng[0]
bd = rng[1]
awsField = "TK" + str(td) + "_" + str(bd) + "S"
arcpy.AddField_management(tmpTable, awsField, "FLOAT", "", "", "")
arcpy.AddField_management(tmpTable, "MUSUMCPCTS", "SHORT", "", "", "")
# Root Zone and root zone available water supply
arcpy.AddField_management(tmpTable, "PCTEARTHMC", "SHORT", "", "", "")
arcpy.AddField_management(tmpTable, "ROOTZNEMC", "SHORT", "", "", "")
arcpy.AddField_management(tmpTable, "ROOTZNAWS", "SHORT", "", "", "")
arcpy.AddField_management(tmpTable, "RESTRICTION", "TEXT", "", "", "254", "RESTRICTION")
# Droughty soils
arcpy.AddField_management(tmpTable, "DROUGHTY", "SHORT", "", "", "")
# Add field for potential wetland soils
arcpy.AddField_management(tmpTable, "PWSL1POMU", "SHORT", "", "", "")
# Add primary key field
arcpy.AddField_management(tmpTable, "MUKEY", "TEXT", "", "", "30", "MUKEY")
# Convert IN_MEMORY table to a permanent table
arcpy.CreateTable_management(outputDB, os.path.basename(theCompTable), tmpTable)
# add attribute indexes for key fields
arcpy.AddIndex_management(theCompTable, "MUKEY", "Indx_Res2Mukey", "NON_UNIQUE", "NON_ASCENDING")
arcpy.AddIndex_management(theCompTable, "COKEY", "Indx_ResCokey", "UNIQUE", "NON_ASCENDING")
# populate table with mukey values
#PrintMsg(" \n\tPopulating " + theCompTable + " with basic component values", 1)
with arcpy.da.SearchCursor(os.path.join(outputDB, "component"), ["mukey", "cokey", "compname", "localphase", "comppct_r"]) as incur:
outcur = arcpy.da.InsertCursor(theCompTable, ["mukey", "cokey", "compname", "localphase", "comppct_r"])
for inrec in incur:
outcur.insertRow(inrec)
except MyError, e:
# Example: raise MyError("this is an error message")
PrintMsg(str(e) + " \n", 2)
return False
except:
errorMsg()
return False
## ===================================================================================
def CheckTexture(mukey, cokey, desgnmaster, om, texture, lieutex, taxorder, taxsubgrp):
# Is this an organic horizon? Look at desgnmaster and OM first. If those
# don't help, look at chtexturegrp.texture next.
#
# if True: Organic, exclude from root zone calculations unless it is 'buried'
# if False: Mineral, include in root zone calculations
#
# 01-26-2015
#
# According to Bob, if TAXORDER = 'Histosol' and DESGNMASTER = 'O' or 'L' then it should NOT be included in the RZAWS calculations
#
# If desgnmast = 'O' or 'L' and not (TAXORDER = 'Histosol' OR TAXSUBGRP like 'Histic%') then exclude this horizon from all RZAWS calcualtions.
#
# lieutext values: Slightly decomposed plant material, Moderately decomposed plant material,
# Bedrock, Variable, Peat, Material, Unweathered bedrock, Sand and gravel, Mucky peat, Muck,
# Highly decomposed plant material, Weathered bedrock, Cemented, Gravel, Water, Cobbles,
# Stones, Channers, Parachanners, Indurated, Cinders, Duripan, Fragmental material, Paragravel,
# Artifacts, Boulders, Marl, Flagstones, Coprogenous earth, Ashy, Gypsiferous material,
# Petrocalcic, Paracobbles, Diatomaceous earth, Fine gypsum material, Undecomposed organic matter
# According to Bob, any of the 'decomposed plant material', 'Muck, 'Mucky peat, 'Peat', 'Coprogenous earth' LIEUTEX
# values qualify.
#
# This function does not determine whether the horizon might be a buried organic. That is done in CalcRZAWS1.
#
lieuList = ['Slightly decomposed plant material', 'Moderately decomposed plant material', \
'Highly decomposed plant material', 'Undecomposed plant material', 'Muck', 'Mucky peat', \
'Peat', 'Coprogenous earth']
txList = ["CE", "COP-MAT", "HPM", "MPM", "MPT", "MUCK", "PDOM", "PEAT", "SPM", "UDOM"]
try:
if str(taxorder) == 'Histosols' or str(taxsubgrp).lower().find('histic') >= 0:
# Always treat histisols and histic components as having all mineral horizons
#if mukey == tmukey:
# PrintMsg("\tHistisol or histic: " + cokey + ", " + str(taxorder) + ", " + str(taxsubgrp), 1)
return False
elif desgnmaster in ["O", "L"]:
# This is an organic horizon according to CHORIZON.DESGNMASTER OR OM_R
#if mukey == tmukey:
# PrintMsg("\tO: " + cokey + ", " + str(taxorder) + ", " + str(taxsubgrp), 1)
return True
#elif om > 19:
# This is an organic horizon according to CHORIZON.DESGNMASTER OR OM_R
# if mukey == tmukey:
# PrintMsg("\tHigh om_r: " + cokey + ", " + str(taxorder) + ", " + str(taxsubgrp), 1)
# return True
elif str(texture) in txList:
# This is an organic horizon according to CHTEXTUREGRP.TEXTURE
#if mukey == tmukey:
# PrintMsg("\tTexture: " + cokey + ", " + str(taxorder) + ", " + str(taxsubgrp), 1)
return True
elif str(lieutex) in lieuList:
# This is an organic horizon according to CHTEXTURE.LIEUTEX
#if mukey == tmukey:
# PrintMsg("\tLieutex: " + cokey + ", " + str(taxorder) + ", " + str(taxsubgrp), 1)
return True
else:
# Default to mineral horizon if it doesn't match any of the criteria
#if mukey == tmukey:
# PrintMsg("\tDefault mineral: " + cokey + ", " + str(taxorder) + ", " + str(taxsubgrp), 1)
return False
except MyError, e:
# Example: raise MyError("this is an error message")
PrintMsg(str(e) + " \n", 2)
return False
except:
errorMsg()
return False
## ===================================================================================
def CheckBulkDensity(sand, silt, clay, bd, mukey, cokey):
# Bob's check for a dense layer
# If sand, silt or clay are missing then we default to Dense layer = False
# If the sum of sand, silt, clay are less than 100 then we default to Dense layer = False
# If a single sand, silt or clay value is NULL, calculate it
try:
#if mukey == tmukey:
# PrintMsg("\tCheck for Dense: " + str(mukey) + ", " + str(cokey) + ", " + \
# str(sand) + ", " + str(silt) + ", " + str(clay) + ", " + str(bd), 1)
txlist = [sand, silt, clay]
if bd is None:
# This is not a Dense Layer
#if mukey == tmukey:
# PrintMsg("\tMissing bulk density", 1)
return False
if txlist.count(None) == 1:
# Missing a single total_r value, calculate it
if txlist[0] is None:
sand = 100.0 - silt - clay
elif silt is None:
silt = 100.0 - sand - clay
else:
clay = 100.0 - sand - silt
txlist = [sand, silt, clay]
if txlist.count(None) > 0:
# Null values for more than one, return False
#if mukey == tmukey:
# PrintMsg("\tDense layer with too many null texture values", 1)
return False
if round(sum(txlist), 1) <> 100.0:
# Cannot run calculation, default value is False
#if mukey == tmukey:
# PrintMsg("\tTexture values do not sum to 100", 1)
return False
# All values required to run the Dense Layer calculation are available
a = bd - ((( sand * 1.65 ) / 100.0 ) + (( silt * 1.30 ) / 100.0 ) + (( clay * 1.25 ) / 100.0))
b = ( 0.002081 * sand ) + ( 0.003912 * silt ) + ( 0.0024351 * clay )
if a > b:
# This is a Dense Layer
#if mukey == tmukey:
# PrintMsg("\tDense layer: a = " + str(a) + " and b = " + str(b) + " and BD = " + str(bd), 1)
return True
else:
# This is not a Dense Layer
#if mukey == tmukey:
# PrintMsg("\tNot a Dense layer: a = " + str(a) + " and b = " + str(b) + " and BD = " + str(bd), 1)
return False
except MyError, e:
# Example: raise MyError("this is an error message")
PrintMsg(str(e) + " \n", 2)
return False
except:
errorMsg()
return False
## ===================================================================================
def CalcRZDepth(inputDB, outputDB, theCompTable, theMuTable, maxD, dPct, dCR):
#
# Look at soil horizon properties to adjust the root zone depth.
# This is in addition to the standard component restrictions
#
# Read the component restrictions into a dictionary, then read through the
# QueryTable_Hz table, calculating the final component rootzone depth
#
# Components with COMPKIND = 'Miscellaneous area' or NULL are filtered out.
# Components with no horizon data are assigned a root zone depth of zero.
#
# Horizons with NULL hzdept_r or hzdepb_r are filtered out
# Horizons with hzdept_r => hzdepb_r are filtered out
# O horizons or organic horizons from the surface down to the first mineral horizon
# are filtered out.
#
# Horizon data below 150cm or select component restrictions are filtered out.
# A Dense layer calculation is also included as an additional horizon-specific restriction.
try:
dComp = dict() # component level data for all component restrictions
dComp2 = dict() # store all component level data plus default values
coList = list()
#bInput = CreateQueryTables1(inputDB, outputDB)
queryTbl = os.path.join(outputDB, "QueryTable_Hz")
#if bInput == False:
# raise MyError, "Problem querying input database"
PrintMsg(" \nProcessing input table (" + queryTbl + ")")
# Create dictionaries and lists
dMapunit = dict() # store mapunit weighted restriction depths
# FIELDS LIST FOR INPUT TABLE
# areasymbol, mukey, musym, muname, mukname,
# cokey, compct, compname, compkind, localphase,
# taxorder, taxsubgrp, ec, pH, dbthirdbar, hzname,
# hzdesgn, hzdept, hzdepb, hzthk, sand,
# silt, clay, om, reskind, reshard,
# resdept, resthk, texture, lieutex
# All reskind values: Strongly contrasting textural stratification, Lithic bedrock, Densic material,
# Ortstein, Permafrost, Paralithic bedrock, Cemented horizon, Undefined, Fragipan, Plinthite,
# Abrupt textural change, Natric, Petrocalcic, Duripan, Densic bedrock, Salic,
# Human-manufactured materials, Sulfuric, Placic, Petroferric, Petrogypsic
#
# Using these restrictions:
# Lithic bedrock, Paralithic bedrock, Densic bedrock, Fragipan, Duripan, Sulfuric
# Other restrictions include pH < 3.5 and EC > 16
resTbl = os.path.join(outputDB, "QueryTable_CR")
crFlds = ["cokey","reskind", "reshard", "resdept_r"]
sqlClause = (None, "ORDER BY cokey, resdept_r ASC")
# ********************************************************
#
# Read the QueryTable_HZ and adjust the component restrictions for additional
# issues such as pH, EC, etc.
#
# Save these new restriction values to dComp dictionary
#
# Only process major-earthy components...
whereClause = "component.compkind <> 'Miscellaneous area' and component.compkind is not Null and component.majcompflag = 'Yes'"
sqlClause = (None, "ORDER BY mukey, comppct_r DESC, cokey, hzdept_r ASC")
curFlds = ["mukey", "cokey", "compname", "compkind", "localphase", "comppct_r", "taxorder", "taxsubgrp", "hzname", "desgnmaster", "hzdept_r", "hzdepb_r", "sandtotal_r", "silttotal_r", "claytotal_r", "om_r", "dbthirdbar_r", "ph1to1h2o_r", "ec_r", "awc_r", "texture", "lieutex"]
resList = ['Lithic bedrock','Paralithic bedrock','Densic bedrock', 'Fragipan', 'Duripan', 'Sulfuric']
lastCokey = "xxxx"
lastMukey = 'xxxx'
# Display status of processing input table containing horizon data and component restrictions
inCnt = int(arcpy.GetCount_management(queryTbl).getOutput(0))
if inCnt > 0:
arcpy.SetProgressor ("step", "Processing input table...", 0, inCnt, 1)
else:
raise MyError, "Input table contains no data"
with arcpy.da.SearchCursor(queryTbl, curFlds, whereClause, "", "", sqlClause) as cur:
# Reading horizon-level data
for rec in cur:
# ********************************************************
#
# Read QueryTable_HZ record
mukey, cokey, compName, compKind, localPhase, compPct, taxorder, taxsubgrp, hzname, desgnmaster, hzDept, hzDepb, sand, silt, clay, om, bd, pH, ec, awc, texture, lieutex = rec
# Initialize component restriction depth to maxD
dComp2[cokey] = [mukey, compName, localPhase, compPct, maxD, ""]
if lastCokey != cokey:
# Accumulate a list of components for future use
lastCokey = cokey
coList.append(cokey)
if hzDept < maxD:
# ********************************************************
# For horizons above the floor level (maxD), look for other restrictive
# layers based on horizon properties such as pH, EC and bulk density.
# Start with the top horizons and work down.
# initialize list of restrictions
resKind = ""
restriction = list()
bOrganic = CheckTexture(mukey, cokey, desgnmaster, om, texture, lieutex, taxorder, taxsubgrp)
if not bOrganic:
# calculate alternate dense layer per Dobos
bDense = CheckBulkDensity(sand, silt, clay, bd, mukey, cokey)
if bDense:
# use horizon top depth for the dense layer
restriction.append("Dense")
resDept = hzDept
# Not sure whether these horizon property checks should be skipped for Organic
# Bob said to only skip Dense Layer check, but VALU table RZAWS looks like all
# horizon properties were skipped.
#
# If we decide to skip EC and pH horizon checks for histosols/histic, use this query
# Example Pongo muck in North Carolina that have low pH but no other restriction
#
if str(taxorder) != 'Histosols' and str(taxsubgrp).lower().find('histic') == -1:
# Only non histosols/histic soils will be checked for pH or EC restrictive horizons
if pH <= 3.5 and pH is not None:
restriction.append("pH")
resDept = hzDept
#if mukey == tmukey:
# PrintMsg("\tpH restriction at " + str(resDept) + "cm", 1)
if ec >= 16.0 and ec is not None:
# Originally I understood that EC > 12 is a restriction, but Bob says he is
# now using 16.
restriction.append("EC")
resDept = hzDept
#if mukey == tmukey:
# PrintMsg("\tEC restriction at " + str(resDept) + "cm", 1)
#if bd >= 1.8:
# restriction.append("BD")
# resDept = hzDept
#if awc is None:
# restriction.append("AWC")
# resDept = hzDept
# ********************************************************
#
# Finally, check for one of the standard component restrictions
#
if cokey in dCR:
resDepth2, resKind = dCR[cokey]
if hzDept <= resDepth2 < hzDepb:
# This restriction may not be at the top of the horizon, thus we
# need to override this if one of the other restrictions exists for this
# horizon
if len(restriction) == 0:
# If this is the only restriction, set the restriction depth
# to the value from the corestriction table.
resDept = resDepth2
# Adding this restriction name to the list even if there are others
# May want to take this out later
restriction.append(resKind)
# ********************************************************
#
if len(restriction) > 0:
# Found at least one restriction for this horizon
if not cokey in dComp:
# if there are no higher restrictions for this component, save this one
# to the dComp dictionary as the top-most restriction
#
dComp[cokey] = [mukey, compName, localPhase, compPct, resDept, restriction]
arcpy.SetProgressorPosition()
arcpy.ResetProgressor()
# Load restrictions from dComp into dComp2 so that there is complete information for all components
for cokey in dComp2:
try:
dComp2[cokey] = dComp[cokey]
except:
pass
# Return the dictionary containing restriction depths and the dictionary containing defaults
return dComp2
except MyError, e:
# Example: raise MyError("this is an error message")
PrintMsg(str(e) + " \n", 2)
return dComp2
except:
errorMsg()
return dComp2
## ===================================================================================
def GetCoRestrictions(outputDB, maxD, resList):
#
# Returns a dictionary of top component restrictions for root growth
#
# resList is a comma-delimited string of reskind values, surrounded by parenthesis
#
# Get component root zone depth from QueryTable_CR and load into dictionary (dCR)
# This is NOT the final root zone depth. This information will be compared with the
# horizon soil properties to determine the final root zone depth.
try:
rSQL = "resdept_r < " + str(maxD) + " and reskind in " + resList
sqlClause = (None, "ORDER BY cokey, resdept_r ASC")
resTbl = os.path.join(outputDB, "QueryTable_CR")
#PrintMsg("\tGetting corestrictions matching: " + resList, 1)
if not arcpy.Exists(resTbl):
raise MyError, "Missing required input table (" + resTbl + ")"
dRestrictions = dict()
# Get the top component restriction from the sorted table
with arcpy.da.SearchCursor(resTbl, ["cokey", "resdept_r", "reskind"], rSQL, "", "", sqlClause) as cur:
for rec in cur:
cokey, resDept, reskind = rec
#PrintMsg("Restriction: " + str(rec), 1)
if not cokey in dRestrictions:
dRestrictions[cokey] = resDept, reskind
return dRestrictions
except MyError, e:
# Example: raise MyError("this is an error message")
PrintMsg(str(e) + " \n", 2)
return dict()
except:
errorMsg()
return dict()
## ===================================================================================
def CalcRZAWS(inputDB, outputDB, td, bd, theCompTable, theMuTable, dRestrictions, maxD, dPct):
# Create a component-level summary table
# Calculate mapunit-weighted average for each mapunit and write to a mapunit-level table
# Need to filter out compkind = 'Miscellaneous area' for RZAWS
# dRestrictions[cokey] = [mukey, compName, localPhase, compPct, resDept, restriction]
try:
import decimal
env.workspace = outputDB
# Using the same component horizon table that has been
queryTbl = os.path.join(outputDB, "QueryTable_Hz")
numRows = int(arcpy.GetCount_management(queryTbl).getOutput(0))
PrintMsg(" \nCalculating Root Zone AWS for " + str(td) + " to " + str(bd) + "cm...", 0)
# QueryTable_HZ fields
qFieldNames = ["mukey", "cokey", "comppct_r", "compname", "localphase", "majcompflag", "compkind", "taxorder", "taxsubgrp", "desgnmaster", "om_r", "awc_r", "hzdept_r", "hzdepb_r", "texture", "lieutex"]
#arcpy.SetProgressorLabel("Creating output tables using dominant component...")
#arcpy.SetProgressor("step", "Calculating root zone available water supply..." , 0, numRows, 1)
# Open edit session on geodatabase to allow multiple update cursors
with arcpy.da.Editor(inputDB) as edit:
# initialize list of components with horizon overlaps
#badCo = list()
# Output fields for root zone and droughty
muFieldNames = ["mukey", "pctearthmc", "rootznemc", "rootznaws", "droughty"]
muCursor = arcpy.da.UpdateCursor(theMuTable, muFieldNames)
# Open component-level output table for updates
#coCursor = arcpy.da.InsertCursor(theCompTable, coFieldNames)
coFieldNames = ["mukey", "cokey", "compname", "localphase", "comppct_r", "pctearthmc", "rootznemc", "rootznaws", "restriction"]
coCursor = arcpy.da.UpdateCursor(theCompTable, coFieldNames)
# Process query table using cursor, write out horizon data for each major component
postFix = "order by mukey, comppct_r DESC, cokey, hzdept_r ASC"
iCnt = int(arcpy.GetCount_management(queryTbl).getOutput(0))
# For root zone calculations, we only want earthy, major components
#PrintMsg(" \nFiltering components in Query_HZ for CalcRZAWS1 function", 1)
#
# Major-Earthy Components
#hzSQL = "component.compkind <> 'Miscellaneous area' and component.compkind is not NULL and component.majcompflag = 'Yes'"
# All Components
hzSQL = ""
inCur = arcpy.da.SearchCursor(queryTbl, qFieldNames, hzSQL, "", False, (None, postFix))
arcpy.SetProgressor("step", "Reading query table...", 0, iCnt, 1)
# Create dictionaries to handle the mapunit and component summaries
dMu = dict()
dComp = dict()
# I may have to pull the sum of component percentages out of this function?
# It seems to work OK for the earthy-major components, but will not work for
# the standard AWS calculations. Those 'Miscellaneous area' components with no horizon data
# are excluded from the Query table because it does not support Outer Joins.
#
mCnt = 0
#PrintMsg("\tmukey, cokey, comppct, top, bottom, resdepth, thickness, aws", 0)
# TEST: keep list of cokeys as a way to track the top organic horizons
skipList = list()
for rec in inCur:
# read each horizon-level input record from QueryTable_HZ ...
#
mukey, cokey, compPct, compName, localPhase, mjrFlag, cKind, taxorder, taxsubgrp, desgnmaster, om, awc, top, bot, texture, lieutex = rec
if mjrFlag == "Yes" and cKind != "Miscellaneous area" and cKind is not None:
# For major-earthy components
# Get restriction information from dictionary
# For non-Miscellaneous areas with no horizon data, set hzdepth values to zero so that
# PWSL and Droughty will get populated with zeros instead of NULL.
if top is None and bot is None:
if not cokey in dComp:
dComp[cokey] = mukey, compName, localPhase, compPct, 0, 0, ""
try:
# mukey, compName, localPhase, compPct, resDept, restriction
# rDepth is the component restriction depth or calculated horizon restriction from CalcRZDepth1 function
# mukey, compName, localPhase, compPct, resDept, restriction] = dRestrictions
d1, d2, d3, d4, rDepth, restriction = dRestrictions[cokey]
cBot = min(rDepth, bot, maxD) # 01-05-2015 Added maxD because I found 46 CONUS mapunits with a ROOTZNEMC > 150
#if mukey == tmukey and rDepth != 150:
# PrintMsg("\tRestriction, " + str(mukey) + ", " + str(cokey) + ", " + str(rDepth) + ", " + str(restriction), 1)
except:
#errorMsg()
cBot = min(maxD, bot)
restriction = []
rDepth = maxD
#if mukey == tmukey:
# PrintMsg("RestrictionError, " + str(mukey) + ", " + str(cokey) + ", " + str(rDepth) + ", " + str(restriction), 1)
bOrganic = CheckTexture(mukey, cokey, desgnmaster, om, texture, lieutex, taxorder, taxsubgrp)
#if mukey == tmukey and bOrganic:
# PrintMsg("Organic: " + str(mukey) + ", " + str(cokey) )
# fix awc_r to 2 decimal places
if awc is None:
awc = 0.0
else:
awc = round(awc, 2)
# Reasons for skipping RZ calculations on a horizon:
# 1. Desgnmaster = O, L and Taxorder != Histosol and is at the surface
# 2. Do I need to convert null awc values to zero?
# 3. Below component restriction or horizon restriction level
if bOrganic and not cokey in skipList:
# Organic surface horizon - Not using this horizon in the calculations
useHz = False
#if mukey == tmukey:
# PrintMsg("Organic, " + str(mukey) + ", " + str(cokey) + ", " + str(compPct) + ", " + str(desgnmaster) + ", " + taxorder + ", " + str(top) + ", " + str(bot) + ", " + str(cBot) + ", " + str(awc) + ", " + str(useHz), 1)
else:
# Mineral, Histosol, buried Organic, Bedrock or there is a horizon restriction (EC, pH - Using this horizon in the calculations
useHz = True
skipList.append(cokey)
# Looking for problems
#if mukey == tmukey:
# PrintMsg("Mineral, " + str(mukey) + ", " + str(cokey) + ", " + str(compPct) + ", " + str(desgnmaster) + ", " + str(taxorder) + ", " + str(top) + ", " + str(bot) + ", " + str(cBot) + ", " + str(awc) + ", " + str(useHz), 1)
# Attempt to fix component with a surface-level restriction that might be in an urban soil
if not cokey in dComp and cBot == 0:
dComp[cokey] = mukey, compName, localPhase, compPct, 0, 0, restriction
# Looking for problems
#if mukey == tmukey:
# PrintMsg("MUKEY2: " + str(mukey) + ", " + str(top) + ", " + str(bot) + ", " + str(cBot) + ", " + str(useHz), 1)
if top < cBot and useHz == True:
# If the top depth is less than the bottom depth, proceed with the calculation
# Calculate sum of horizon thickness and sum of component ratings for all horizons above bottom
hzT = cBot - top
aws = float(hzT) * float(awc) * 10.0
# Looking for problems
#if mukey == tmukey:
# PrintMsg("MUKEY3: " + str(mukey) + ", " + str(top) + ", " + str(bot) + ", " + str(cBot) + ", " + str(useHz), 1)
if cokey in dComp:
# accumulate total thickness and total rating value by adding to existing component values
mukey, compName, localPhase, compPct, dHzT, dAWS, restriction = dComp[cokey]
dAWS = dAWS + aws
dHzT += hzT
dComp[cokey] = mukey, compName, localPhase, compPct, dHzT, dAWS, restriction
# Looking for problems
#if mukey == tmukey:
# PrintMsg("MUKEY4: " + str(mukey) + ", " + str(top) + ", " + str(bot) + ", " + str(cBot) + ", " + str(useHz), 1)
#if dHzT > 150.0:
# overlap = dHzT - 150
# Found some components where there are overlapping horizons
#badCo.append(str(cokey))
else:
# Create initial entry for this component using the first horizon
dComp[cokey] = mukey, compName, localPhase, compPct, hzT, aws, restriction
# Looking for problems
#if mukey == tmukey:
# PrintMsg("MUKEY5: " + str(mukey) + ", " + str(top) + ", " + str(bot) + ", " + str(cBot) + ", " + str(useHz), 1)
else:
# Do not include this horizon in the rootzone calculations
pass
else:
# Not a major-earthy component, so write out everything BUT rzaws-related data (last values)
dComp[cokey] = mukey, compName, localPhase, compPct, None, None, None, None
arcpy.SetProgressorPosition()
# end of processing major-earthy components
arcpy.ResetProgressor()
# get the total number of major-earthy components from the dictionary count
iComp = len(dComp)
# Read through the component-level data and summarize to the mapunit level
if iComp > 0:
PrintMsg(" \nSaving component average RZAWS to table... (" + str(iComp) + ")", 0 )
arcpy.SetProgressor("step", "Saving component data...", 0, iComp, 1)
iCo = 0 # count component records written to theCompTbl
for corec in coCursor:
mukey, cokey, compName, localPhase, compPct, pctearthmc, rDepth, aws, restrictions = corec
try:
# get sum of component percent for the mapunit
pctearthmc = float(dPct[mukey][1]) # sum of comppct_r for all major components Test 2014-10-07
# get rootzone data from dComp
mukey1, compName1, localPhase1, compPct1, hzT, awc, restriction = dComp[cokey]
except:
pctearthmc = 0
hzT = None
rDepth = None
awc = None
restriction = []
# calculate component percentage adjustment
if pctearthmc > 0 and not awc is None:
# If there is no data for any of the component horizons, could end up with 0 for
# sum of comppct_r
adjCompPct = float(compPct) / float(pctearthmc)
# adjust the rating value down by the component percentage and by the sum of the usable horizon thickness for this component
aws = adjCompPct * float(awc) # component rating
if restriction is None:
restrictions = ''
elif len(restriction) > 0:
restrictions = ",".join(restriction)
else:
restrictions = ''
corec = mukey, cokey, compName, localPhase, compPct, pctearthmc, hzT, aws, restrictions
coCursor.updateRow(corec)
iCo += 1
# Weight hzT for ROOTZNEMC by component percent
hzT = (float(hzT) * float(compPct) / pctearthmc)
if mukey in dMu:
val1, val2, val3 = dMu[mukey]
dMu[mukey] = pctearthmc, (hzT + val2), (aws + val3)
else:
# first entry for map unit ratings
dMu[mukey] = pctearthmc, hzT, aws
else:
# Populate component level record for a component with no AWC
corec = mukey, cokey, compName, localPhase, compPct, None, None, None, ""
coCursor.updateRow(corec)
iCo += 1
arcpy.SetProgressorPosition()
arcpy.ResetProgressor()
else:
raise MyError, "No component data in dictionary dComp"
if len(dMu) > 0:
PrintMsg(" \nSaving map unit average RZAWS to table...(" + str(len(dMu)) + ")", 0 )
else:
raise MyError, "No map unit information in dictionary dMu"
# Save root zone available water supply and droughty soils to output map unit table
#
for murec in muCursor:
mukey, pctearthmc, rootznemc, rootznaws, droughty = murec
try:
rec = dMu[mukey]
pct, rootznemc, rootznaws = rec
pctearthmc = dPct[mukey][1]
if rootznemc > 150.0:
# This is a bandaid for components that have horizon problems such
# overlapping that causes the calculated total to exceed 150cm.
rootznemc = 150.0
rootznaws = round(rootznaws, 0)
rootznemc = round(rootznemc, 0)
if rootznaws > 152:
droughty = 0
else:
droughty = 1
except:
pctearthmc = 0
rootznemc = None
rootznaws = None
murec = mukey, pctearthmc, rootznemc, rootznaws, droughty
muCursor.updateRow(murec)
#if mukey == tmukey:
# values at this point seem to be correct
# fldnames = muCursor.fields
# PrintMsg(str(fldnames), 1)
# PrintMsg(str(murec), 1)
# Save data issues to permanent files for later review
#if len(badCo) > 0:
# fileCo = os.path.basename(inputDB)[:-4] + "_OverlappingHz.txt"
# fileCo = os.path.join(os.path.dirname(inputDB), fileCo)
# fh = open(fileCo, "w")
# fh.write(inputDB + "\n")
# fh.write("Components with overlapping horizons\n\n")
# fh.write("COKEY IN ('" + "', '".join(badCo) + "') \n")
# fh.close()
# PrintMsg(" \nComponents with overlapping horizons (" + Number_Format(len(badCo), 0, True) + ") saved to:\t" + fileCo, 0)
PrintMsg("", 0)
return True
except MyError, e:
# Example: raise MyError("this is an error message")
PrintMsg(str(e) + " \n", 2)
return False
except:
errorMsg()
return False
## ===================================================================================
def CalcAWS(inputDB, outputDB, theCompTable, theMuTable, dPct):
# Create a component-level summary table
# Calculate the standard mapunit-weighted available waters supply for each mapunit and
# add it to the map unit-level table.
#
# 12-08 I see that for mukey='2479901' my rating is
try:
# Using the same component horizon table that has been
queryTbl = os.path.join(outputDB, "QueryTable_HZ")
numRows = int(arcpy.GetCount_management(queryTbl).getOutput(0))
# mukey, cokey, compPct,val, top, bot
qFieldNames = ["mukey", "cokey", "comppct_r", "awc_r", "hzdept_r", "hzdepb_r"]
# Track map units that are missing data
missingList = list()
minusList = list()
PrintMsg(" \nCalculating standard available water supply for:", 0)
for rng in depthList:
# Calculating and updating just one AWS column at a time
#
td = rng[0]
bd = rng[1]
#outputFields = "AWS" + str(td) + "_" + str(bd), "TK" + str(td) + "_" + str(bd) + "A"
# Open output table Mu...All in write mode
muFieldNames = ["MUKEY", "MUSUMCPCTA", "AWS" + str(td) + "_" + str(bd), "TK" + str(td) + "_" + str(bd) + "A"]
coFieldNames = ["COKEY", "AWS" + str(td) + "_" + str(bd), "TK" + str(td) + "_" + str(bd) + "A"]
# Create dictionaries to handle the mapunit and component summaries
dMu = dict()
dComp = dict()
dSum = dict() # store sum of comppct_r and total thickness for the component
dHz = dict() # Trying a new dictionary that will s
arcpy.SetProgressorLabel("Creating output tables using dominant component...")
arcpy.SetProgressor("step", "Aggregating data for the dominant component..." , 0, numRows, 1)
# Open edit session on geodatabase to allow multiple insert cursors
with arcpy.da.Editor(inputDB) as edit:
# Open output mapunit-level table in update mode
# MUKEY, AWS
muCursor = arcpy.da.UpdateCursor(theMuTable, muFieldNames)
# Open output component-level table in write mode
# MUKEY, AWS
coCursor = arcpy.da.UpdateCursor(theCompTable, coFieldNames)
# Process query table using a searchcursor, write out horizon data for each component
# At this time, almost all components are being used! There is no filter.
postFix = "order by mukey, comppct_r DESC, cokey, hzdept_r ASC"
#hzSQL = "compkind is not null and hzdept_r is not null" # prevent divide-by-zero errors
hzSQL = "hzdept_r is not null" # prevent divide-by-zero errors by skipping components with no horizons
iCnt = int(arcpy.GetCount_management(queryTbl).getOutput(0))
inCur = arcpy.da.SearchCursor(queryTbl, qFieldNames, hzSQL, "", False, (None, postFix))
arcpy.SetProgressor("step", "Reading QueryTable_HZ ...", 0, iCnt, 1)
for rec in inCur:
# read each horizon-level input record from the query table ...
mukey, cokey, compPct, awc, top, bot = rec
if awc is not None:
# Calculate sum of horizon thickness and sum of component ratings for all horizons above bottom
hzT = min(bot, bd) - max(top, td) # usable thickness from this horizon
if hzT > 0:
aws = float(hzT) * float(awc) * 10
if not cokey in dComp:
# Create initial entry for this component using the first horiozon CHK
dComp[cokey] = (mukey, compPct, hzT, aws)
else:
# accumulate total thickness and total rating value by adding to existing component values CHK
mukey, compName, dHzT, dAWS = dComp[cokey]
dAWS = dAWS + aws
dHzT = dHzT + hzT
dComp[cokey] = (mukey, compPct, dHzT, dAWS)
arcpy.SetProgressorPosition()
# get the total number of major components from the dictionary count
iComp = len(dComp)
# Read through the component-level data and summarize to the mapunit level
if iComp > 0:
PrintMsg("\t" + str(td) + " - " + str(bd) + "cm (" + Number_Format(iComp, 0, True) + " components)" , 0)
arcpy.SetProgressor("step", "Saving map unit and component AWS data...", 0, iComp, 1)
for corec in coCursor:
# get component level data CHK
cokey = corec[0]
if cokey in dComp:
dRec = dComp[corec[0]]
mukey, compPct, hzT, awc = dRec
# get sum of component percent for the mapunit CHK
try:
# Value[0] is for all components,
# Value[1] is just for major-earthy components,
# Value[2] is all major components
# Value[3] is earthy components
sumCompPct = float(dPct[mukey][0])
#sumCompPct = float(dPct[mukey][1])
except:
# set the component percent to zero if it is not found in the
# dictionary. This is probably a 'Miscellaneous area' not included in the CHK
# data or it has no horizon information.
sumCompPct = 0
#missingList.append("'" + mukey + "'")
# calculate component percentage adjustment
if sumCompPct > 0:
# If there is no data for any of the component horizons, could end up with 0 for
# sum of comppct_r
#PrintMsg(" \nMUKEY " + mukey + " - " + compName + " has zero percent Sum Comppct", 1)
#adjCompPct = float(compPct) / sumCompPct # WSS method
adjCompPct = compPct / 100.0 # VALU table method
# adjust the rating value down by the component percentage and by the sum of the usable horizon thickness for this component
aws = round((adjCompPct * awc), 2) # component rating
corec[1] = aws
hzT = hzT * adjCompPct # Adjust component share of horizon thickness by comppct
corec[2] = hzT # This is new for the TK0_5A column
coCursor.updateRow(corec)
# Update component values in component dictionary CHK
# Not sure what dComp is being used for ???
dComp[cokey] = mukey, compPct, hzT, aws
# Try to fix high mapunit aggregate HZ by weighting with comppct
# Testing new mapunit aggregation 09-08-2014
# Trying to replace dMu dictionary
if mukey in dMu:
val1, val2, val3 = dMu[mukey]
#dMu[mukey] = (compPct + val1, hzT + val2, aws + val2)
compPct = compPct + val1
hzT = hzT + val2
aws = aws + val3
#else:
dMu[mukey] = (compPct, hzT, aws)
#if mukey == '2479892':
arcpy.SetProgressorPosition()
arcpy.ResetProgressor()
else:
PrintMsg("\t" + Number_Format(iComp, 0, True) + " components for " + str(td) + " - " + str(bd) + "cm", 1)
# Write out map unit aggregated AWS
#
for murec in muCursor:
mukey = murec[0]
if mukey in dMu:
compPct, hzT, aws = dMu[mukey]
murec[1] = compPct
murec[2] = aws
murec[3] = round(hzT, 2) # sometimes this ends up being 2 or 3X what it should
muCursor.updateRow(murec)
if len(missingList) > 0:
missingList = list(set(missingList))
PrintMsg(" \nFollowing mapunits have no comppct_r: " + ", ".join(missingList), 1)
PrintMsg("", 0)
return True
except MyError, e:
# Example: raise MyError("this is an error message")
PrintMsg(str(e) + " \n", 2)
return False
except:
errorMsg()
return False
## ===================================================================================
def CalcSOC(inputDB, outputDB, theCompTable, theMuTable, dPct, dFrags, depthList, dRestrictions, maxD):
# Create a component-level summary table
# Calculate the standard mapunit-weighted available waters supply for each mapunit and
# add it to the map unit-level table
#
# Current problems: have MuSOC with multiple records for mukey=2479888
try:
# Using the same component horizon table that has been
queryTbl = os.path.join(outputDB, "QueryTable_HZ")
numRows = int(arcpy.GetCount_management(queryTbl).getOutput(0))
# mukey, cokey, compPct,val, top, bot
#qFieldNames = ["mukey", "cokey", "comppct_r", "hzdept_r", "hzdepb_r", "om_r", "dbthirdbar_r"]
qFieldNames = ["mukey","cokey","comppct_r","compname","localphase","chkey","om_r","dbthirdbar_r", "hzdept_r","hzdepb_r"]
# Track map units that are missing data
missingList = list()
minusList = list()
PrintMsg(" \nCalculating soil organic carbon for:", 0)
for rng in depthList:
# Calculating and updating just one SOC column at a time
#
td = rng[0]
bd = rng[1]
# Open output table Mu...All in write mode
# I lately added the "MUSUMCPCTS" to the output. Need to check output because
# it will be writing this out for every range. Lots more overhead.
#
# Create dictionaries to handle the mapunit and component summaries
dMu = dict()
dComp = dict()
#dSumPct = dict() # store the sum of comppct_r for each mapunit to use in the calculations
dSum = dict() # store sum of comppct_r and total thickness for the component
dHz = dict() # Trying a new dictionary that will s
mCnt = 0
arcpy.SetProgressorLabel("Creating output tables using dominant component...")
arcpy.SetProgressor("step", "Aggregating data for the dominant component..." , 0, numRows, 1)
# Open edit session on geodatabase to allow multiple insert cursors
with arcpy.da.Editor(inputDB) as edit:
# Open output mapunit-level table in update mode
muFieldNames = ["MUKEY", "MUSUMCPCTS", "SOC" + str(td) + "_" + str(bd), "TK" + str(td) + "_" + str(bd) + "S"]
muCursor = arcpy.da.UpdateCursor(theMuTable, muFieldNames)
# Open output component-level table in write mode
coFieldNames = ["COKEY", "SOC" + str(td) + "_" + str(bd), "TK" + str(td) + "_" + str(bd) + "S"]
coCursor = arcpy.da.UpdateCursor(theCompTable, coFieldNames)
# Process query table using a searchcursor, write out horizon data for each component
# At this time, almost all components are being used! There is no filter.
postFix = "order by mukey, comppct_r DESC, cokey, hzdept_r ASC"
#hzSQL = "compkind is not null and hzdept_r is not null" # prevent divide-by-zero errors
hzSQL = "hzdept_r is not null" # prevent divide-by-zero errors by skipping components with no horizons
iCnt = int(arcpy.GetCount_management(queryTbl).getOutput(0))
inCur = arcpy.da.SearchCursor(queryTbl, qFieldNames, hzSQL, "", False, (None, postFix))
arcpy.SetProgressor("step", "Reading QueryTable_HZ ...", 0, iCnt, 1)
for rec in inCur:
# read each horizon-level input record from the query table ...
mukey, cokey, compPct, compName, localPhase, chkey, om, db3, top, bot = rec
if om is not None and db3 is not None:
# Calculate sum of horizon thickness and sum of component ratings for
# that portion of the horizon that is with in the td-bd range
top = max(top, td)
bot = min(bot, bd)
om = round(om, 3)
try:
rz, resKind = dRestrictions[cokey]
except:
rz = maxD
resKind = ""
# Now check for horizon restrictions within this range
if top < rz < bot:
# restriction found in this horizon, use it to set a new depth
#PrintMsg("\t\t" + resKind + " restriction for " + mukey + ":" + cokey + " at " + str(rz) + "cm", 1)
cBot = rz
else:
cBot = min(rz, bot)
# Calculate initial usable horizon thickness
hzT = cBot - top
if hzT > 0 and top < cBot:
# get horizon fragment volume
try:
fragvol = dFrags[chkey]
except:
fragvol = 0.0
pass
# Calculate SOC using horizon thickness, OM, BD, FragVol, CompPct
soc = ( (hzT * ( ( om * 0.58 ) * db3 )) / 100.0 ) * ((100.0 - fragvol) / 100.0) * ( compPct * 100 )
if not cokey in dComp:
# Create initial entry for this component using the first horizon CHK
dComp[cokey] = (mukey, compPct, hzT, soc)
else:
# accumulate total thickness and total rating value by adding to existing component values CHK
mukey, compName, dHzT, dSOC = dComp[cokey]
dSOC = dSOC + soc
dHzT = dHzT + hzT
dComp[cokey] = (mukey, compPct, dHzT, dSOC)
arcpy.SetProgressorPosition()
# get the total number of major components from the dictionary count
iComp = len(dComp)
# Read through the component-level data and summarize to the mapunit level
#
if iComp > 0:
PrintMsg("\t" + str(td) + " - " + str(bd) + "cm (" + Number_Format(iComp, 0, True) + " components)", 0)
arcpy.SetProgressor("step", "Saving map unit and component SOC data...", 0, iComp, 1)
for corec in coCursor:
# Could this be where I am losing minor components????
#
# get component level data CHK
cokey = corec[0]
if cokey in dComp:
# get SOC-related data from dComp by cokey
# dComp soc = ( (hzT * ( ( om * 0.58 ) * db3 )) / 100.0 ) * ~
# ((100.0 - fragvol) / 100.0) * ( compPct * 100 )
mukey, compPct, hzT, soc = dComp[corec[0]]
# get sum of component percent for the mapunit (all components???)
# Value[0] is for all components,
# Value[1] is just for major-earthy components,
# Value[2] is all major components
# Value[3] is earthy components
try:
sumCompPct = float(dPct[mukey][0]) # Sum comppct for ALl components
# Sum of comppct for Major-earthy components
#sumCompPct = float(dPct[mukey][1])
except:
# set the component percent to zero if it is not found in the
# dictionary. This is probably a 'Miscellaneous area' not included in the CHK
# data or it has no horizon information.
sumCompPct = 0.0
#missingList.append("'" + mukey + "'")
# calculate component percentage adjustment
if sumCompPct > 0:
# adjust the rating value down by the component percentage and by the sum of the usable horizon thickness for this component
#soc = round((adjCompPct * om), 2) # component rating
#soc = round((100.0 * om / sumCompPct), 0) # Try adjusting soc down by the sum of the component percents
adjCompPct = float(compPct) / sumCompPct #
# write the new component-level SOC data to the Co_VALU table
corec[1] = soc # Test
hzT = hzT * compPct / 100.0 # Adjust component share of horizon thickness by comppct
corec[2] = hzT # This is new for the TK0_5A column
coCursor.updateRow(corec)
# Update component values in component dictionary CHK
dComp[cokey] = mukey, compPct, hzT, soc
# Try to fix high mapunit aggregate HZ by weighting with comppct
# Testing new mapunit aggregation 09-08-2014
# Trying to replace dMu dictionary
if mukey in dMu:
val1, val2, val3 = dMu[mukey]
#dMu[mukey] = (compPct + val1, hzT + val2, aws + val2)
compPct = compPct + val1
hzT = hzT + val2
soc = soc + val3
dMu[mukey] = (compPct, hzT, soc) # SOC is a little low
#dMu[mukey] = (compPct, hzT, ( soc / adjCompPct))
arcpy.SetProgressorPosition()
arcpy.ResetProgressor()
else:
PrintMsg("\t" + Number_Format(iComp, 0, True) + " components for " + str(td) + " - " + str(bd) + "cm", 1)
# Write out map unit aggregated AWS
#
for murec in muCursor:
mukey = murec[0]
if mukey in dMu:
compPct, hzT, soc = dMu[mukey]
murec[1] = compPct
murec[2] = round(soc, 0)
murec[3] = round(hzT, 0) # this value appears to be low sometimes
muCursor.updateRow(murec)
#if len(missingList) > 0:
# missingList = list(set(missingList))
# PrintMsg(" \nFollowing mapunits have no comppct_r: " + ", ".join(missingList), 1)
#PrintMsg("", 0)
return True
except MyError, e:
# Example: raise MyError("this is an error message")
PrintMsg(str(e) + " \n", 2)
return False
except:
errorMsg()
return False
## ===================================================================================
def GetFragVol(inputDB):
# Get the horizon summary of rock fragment volume (percent)
# load sum of comppct_r into a dictionary by chkey. This
# value will be used to reduce amount of SOC for each horizon
# If not all horizons are not present in the dictionary, failover to
# zero for the fragvol value.
try:
fragSQL = ""
fragFlds = ["chkey", "fragvol_r"]
dFrags = dict()
with arcpy.da.SearchCursor(os.path.join(inputDB, "chfrags"), fragFlds, fragSQL) as fragCur:
for rec in fragCur:
chkey, fragvol = rec
if chkey in dFrags:
# This horizon has a volume already
# Get the existing values
# limit total fragvol to 100 or we will get negative SOC values where there
# are problems with fragvol data
val = dFrags[chkey]
dFrags[chkey] = min(val + max(fragvol, 0), 100)
else:
# this is the first component for this map unit
dFrags[chkey] = min(max(fragvol, 0), 100)
# in the rare case where fragvol sum is greater than 100%, return 100
return dFrags
except MyError, e:
# Example: raise MyError("this is an error message")
PrintMsg(str(e) + " \n", 2)
return dict()
except:
errorMsg()
return dict()
## ===================================================================================
def GetSumPct(inputDB):
# Get map unit - sum of component percent for all components and also for major-earthy components
# load sum of comppct_r into a dictionary.
# Value[0] is for all components,
# Value[1] is just for major-earthy components,
# Value[2] is all major components
# Value[3] is earthy components
#
# Do I need to add another option for earthy components?
# WSS and SDV use all components with data for AWS.
try:
pctSQL = "comppct_r is not null"
pctFlds = ["mukey", "compkind", "majcompflag", "comppct_r"]
dPct = dict()
with arcpy.da.SearchCursor(os.path.join(inputDB, "component"), pctFlds, pctSQL) as pctCur:
for rec in pctCur:
mukey, compkind, flag, comppct = rec
m = 0 # major component percent
me = 0 # major-earthy component percent
e = 0 # earthy component percent
if flag == 'Yes':
# major component percent
m = comppct
if not compkind in ["Miscellaneous area", ""]:
# major-earthy component percent
me = comppct
e = comppct
else:
me = 0
elif not compkind in ["Miscellaneous area", ""]:
e = comppct
if mukey in dPct:
# This mapunit has a pair of values already
# Get the existing values from the dictionary
#pctAll, pctMjr = dPct[mukey] # all components, major-earthy
pctAll, pctME, pctMjr, pctE = dPct[mukey]
dPct[mukey] = (pctAll + comppct, pctME + me, pctMjr + m, pctE + e)
else:
# this is the first component for this map unit
dPct[mukey] = (comppct, me, m, e)
return dPct
except MyError, e:
# Example: raise MyError("this is an error message")
PrintMsg(str(e) + " \n", 2)
return dict()
except:
errorMsg()
return dict()
## ===================================================================================
def MakeNCCPIQueryTable(inputDB, qTable):
# create query table containing information from component and chorizon tables
# return name of querytable. Failure returns an empty string for the table name.
#
# COINTERP.RULENAME CHOICES:
#
# 'NCCPI - National Commodity Crop Productivity Index (Ver 2.0)'
# 'NCCPI - NCCPI Corn and Soybeans Submodel (II)'
# 'NCCPI - NCCPI Cotton Submodel (II)'
# 'NCCPI - NCCPI Small Grains Submodel (II)'
#
try:
# Join chorizon table with component table
inTables = [os.path.join(inputDB, "component"), os.path.join(inputDB, "cointerp")]
# interphr is the fuzzy value
theFields = [["COMPONENT.MUKEY", "MUKEY"], \
["COMPONENT.COKEY", "COKEY"], \
["COMPONENT.COMPPCT_R", "COMPPCT_R"], \
["COINTERP.RULENAME", "RULENAME"], \
["COINTERP.RULEDEPTH", "RULEDEPTH"], \
["COINTERP.INTERPHR", "INTERPHR"]]
rule = 'NCCPI - National Commodity Crop Productivity Index (Ver 2.0)'
theSQL = "COMPONENT.COMPPCT_R > 0 AND COMPONENT.MAJCOMPFLAG = 'Yes' AND COMPONENT.COKEY = COINTERP.COKEY AND COINTERP.MRULENAME = '" + rule + "'"
PrintMsg(" \nCalculating NCCPI weighted averages for all major components...", 0)
# Things to be aware of with MakeQueryTable:
# USE_KEY_FIELDS does not create OBJECTID field. Lack of OBJECTID precludes sorting on Mukey.
# ADD_VIRTUAL_KEY_FIELD creates OBJECTID, but qualifies field names using underscore (eg. COMPONENT_COKEY)
#
arcpy.MakeQueryTable_management(inTables, qTable, "ADD_VIRTUAL_KEY_FIELD","",theFields, theSQL)
if arcpy.Exists(qTable):
iCnt = int(arcpy.GetCount_management(qTable).getOutput(0))
if iCnt == 0:
raise MyError, "Failed to retrieve NCCPI data"
return True
except MyError, e:
# Example: raise MyError("this is an error message")
PrintMsg(str(e) + " \n", 2)
return False
except:
errorMsg()
return False
## ===================================================================================
def CalcNCCPI(inputDB, qTable, dPct):
#
#
try:
# and write to Mu_NCCPI2 table
#
PrintMsg(" \nAggregating data to mapunit level...", 0)
# Alternate component fields for all NCCPI values
cFlds = ["MUKEY","COKEY","COMPPCT_R","COMPNAME","LOCALPHASE", "DUMMY"]
mFlds = ["MUKEY","COMPPCT_R","COMPNAME","LOCALPHASE", "DUMMY"]
# Create dictionary key as MUKEY:INTERPHRC
# Need to look through the component rating class for ruledepth = 0
# and sum up COMPPCT_R for each key value
#
dVals = dict() # dictionary containing sum of comppct for each MUKEY:RATING combination
# Get sum of component percent for each map unit. There are different options:
# 1. Use all major components
# 2. Use all components that have an NCCPI rating
# 3. Use all major-earthy components. This one is not currently available.
# 4. Use all components (that have a component percent)
#
# Query table fields
qFields = ["COMPONENT_MUKEY", "COMPONENT_COKEY", "COMPONENT_COMPPCT_R", "COINTERP_RULEDEPTH", "COINTERP_RULENAME", "COINTERP_INTERPHR"]
sortFields = "ORDER BY COMPONENT_COKEY ASC, COMPONENT_COMPPCT_R DESC"
querytblSQL = "COMPONENT_COMPPCT_R IS NOT NULL"
iCnt = int(arcpy.GetCount_management(qTable).getOutput(0))
noVal = list() # Get a list of components with no overall index rating
PrintMsg(" \nReading query table with " + Number_Format(iCnt, 0, True) + " records...", 0)
arcpy.SetProgressor("step", "Reading query table...", 0,iCnt, 1)
with arcpy.da.SearchCursor(qTable, qFields, querytblSQL, "","", (None, sortFields)) as qCursor:
for qRec in qCursor:
# qFields = MUKEY, COKEY, COMPPCT_R, COMPNAME, LOCALPHASE, RULEDEPTH, RULENAME, INTERPHR
mukey, cokey, comppct, ruleDepth, ruleName, fuzzyValue = qRec
# Dictionary order: All, CS, CT, SG
if not mukey in dVals:
# Initialize mukey NCCPI values
dVals[mukey] = [None, None, None, None]
if not fuzzyValue is None:
if ruleDepth == 0:
# This is NCCPI Overall Index
oldVal = dVals[mukey][0]
if oldVal is None:
dVals[mukey][0] = fuzzyValue * comppct
else:
dVals[mukey][0] = (oldVal + (fuzzyValue * comppct))
elif ruleName == "NCCPI - NCCPI Corn and Soybeans Submodel (II)":
oldVal = dVals[mukey][1]
if oldVal is None:
dVals[mukey][1] = fuzzyValue * comppct
else:
dVals[mukey][1] = (oldVal + (fuzzyValue * comppct))
elif ruleName == "NCCPI - NCCPI Cotton Submodel (II)":
oldVal = dVals[mukey][2]
if oldVal is None:
dVals[mukey][2] = fuzzyValue * comppct
else:
dVals[mukey][2] = (oldVal + (fuzzyValue * comppct))
elif ruleName == "NCCPI - NCCPI Small Grains Submodel (II)":
oldVal = dVals[mukey][3]
if oldVal is None:
dVals[mukey][3] = fuzzyValue * comppct
else:
dVals[mukey][3] = (oldVal + (fuzzyValue * comppct))
elif ruleName == "NCCPI - National Commodity Crop Productivity Index (Ver 2.0)":
# This component does not have an NCCPI rating
#PrintMsg(" \n" + mukey + ":" + cokey + ", " + str(comppct) + "% has no NCCPI rating", 1)
noVal.append("'" + cokey + "'")
arcpy.SetProgressorPosition()
#
# End of query table iteration
#
#if len(noVal) > 0:
# PrintMsg(" \nThe following components had no NCCPI overall index: " + ", ".join(noVal), 1)
iCnt = len(dVals)
if iCnt > 0:
PrintMsg(" \nSaving map unit weighted NCCPI data (" + Number_Format(iCnt, 0, True) + " records) to " + os.path.basename(theMuTable) + "..." , 0)
# Write map unit aggregate data to Mu_NCCPI2 table
#
# theMuTable is a global variable. Need to check this out in the gSSURGO_ValuTable script
with arcpy.da.UpdateCursor(theMuTable, ["mukey", "NCCPI2CS", "NCCPI2CO","NCCPI2SG", "NCCPI2ALL"]) as muCur:
arcpy.SetProgressor("step", "Saving map unit weighted NCCPI data to VALU table...", 0, iCnt, 0)
for rec in muCur:
mukey = rec[0]
try:
# Get output values from dVals and dPct dictionaries
val = dVals[mukey]
ovrall, cs, co, sg = val
sumPct = dPct[mukey][2] # sum of major-earthy components
rec = mukey, round(cs / sumPct, 2), round(co / sumPct, 2), round(sg / sumPct, 2), round(ovrall / sumPct, 2)
muCur.updateRow(rec)
#PrintMsg(mukey + ", " + str(round(cs / sumPct, 2)) + ", " + str(round(co / sumPct, 2)) + ", " + str(round(sg / sumPct, 2)) + ", " + str(round(ovrall / sumPct, 2)), 1)
except:
# Miscellaneous map unit encountered with no comppct_r?
pass
arcpy.SetProgressorPosition()
arcpy.Delete_management(qTable)
return True
else:
raise MyError, "No NCCPI data processed"
except MyError, e:
# Example: raise MyError("this is an error message")
PrintMsg(str(e) + " \n", 2)
return False
except:
errorMsg()
return False
## ===================================================================================
def CalcPWSL(inputDB, outputDB, theMuTable, dPct):
# Get potential wet soil landscape rating for each map unit
# Assuming that all components (with comppct_r) will be processed
#
# Sharon: I treat all map unit components the same, so if I find 1% water I think
# it should show up as 1% PWSL. If the percentage of water is >= 80% then I class
# it into the water body category or 999.
# Steve: In looking at some of Sharon's PWSL=999, I find some mapunits where the
# COMPNAME='Water' at 100%, and HYDRIC='Unranked'. My Valu2 table PWSL is null. That isn't right.
try:
env.workspace = outputDB
# Using the same component horizon table as always
queryTbl = os.path.join(outputDB, "QueryTable_Hz")
numRows = int(arcpy.GetCount_management(queryTbl).getOutput(0))
PrintMsg(" \nCalculating Potential Wet Soil Landscapes using " + queryTbl + "...", 0)
qFieldNames = ["mukey", "muname", "cokey", "comppct_r", "compname", "localphase", "otherph", "majcompflag", "compkind", "hydricrating", "drainagecl"]
pwSQL = ""
compList = list()
dMu = dict()
drainList = ["Poorly drained", "Very poorly drained"]
phaseList = ["drained", "undrained", "channeled", "protected", "ponded", "flooded"]
iCnt = int(arcpy.GetCount_management(queryTbl).getOutput(0))
lastCokey = 'xxx'
arcpy.SetProgressor("step", "Reading query table table for wetland information...", 0, iCnt, 1)
with arcpy.da.SearchCursor(queryTbl, qFieldNames, pwSQL) as pwCur:
for rec in pwCur:
mukey, muname, cokey, comppct_r, compname, localphase, otherph, majcompflag, compkind, hydricrating, drainagecl = rec
#if mukey == tmukey:
# PrintMsg(" \nGetting PWSL data for test map unit: " + str(mukey), 1)
if cokey != lastCokey:
# only process first horizon record for each component
compList.append(cokey)
# Only check the first horizon record, really only need component level
# Not very efficient, should problably create a new query table
#
# Need to split up these tests so that None types can be handled
# Sharon says that if the hydricrating for a component is 'No', don't
# look at it any further. If it is unranked, go ahead and look at
# other properties.
#
pw = False
if muname == 'Water' or compname == 'Water':
# Check for water before looking at Hydric rating
# Probably won't catch everything. Waiting for Sharon's criteria.
if comppct_r >= 80:
# Flag this mapunit with a '999'
# Not necessarily catching map unit with more than one Water component that
# might sum to >= 80.
pw = False
dMu[mukey] = 999
else:
pw = True
try:
sumPct = dMu[mukey]
if sumPct != 999:
dMu[mukey] = sumPct + comppct_r
except:
dMu[mukey] = comppct_r
elif hydricrating == 'No':
# Added this bit so that other properties cannot override hydricrating = 'No'
pw = False
elif hydricrating == 'Yes':
# This is always a Hydric component
# Get component percent and add to map unit total PWSL
pw = True
#if mukey == tmukey:
# PrintMsg("\tHydric percent = " + str(comppct_r), 1)
try:
sumPct = dMu[mukey]
if sumPct != 999:
dMu[mukey] = sumPct + comppct_r
except:
dMu[mukey] = comppct_r
elif hydricrating == 'Unranked':
# Not sure how Sharon is handling NULL hydric
#
# Unranked hydric from here on down, looking at other properties such as:
# Local phase
# Other phase
# Drainage class
# Map unit name strings
#if localphase is not None:
if [d for d in phaseList if str(localphase).lower().find(d) >= 0]:
pw = True
#if mukey == tmukey:
# PrintMsg("\tLocalPhase:" + localphase, 1)
try:
sumPct = dMu[mukey]
dMu[mukey] = sumPct + comppct_r
except:
dMu[mukey] = comppct_r
# otherphase
#if pw == False and otherph is not None:
elif [d for d in phaseList if str(otherph).lower().find(d) >= 0]:
pw = True
#if mukey == tmukey:
# PrintMsg("\tOtherPhase:" + localphase, 1)
try:
sumPct = dMu[mukey]
dMu[mukey] = sumPct + comppct_r
except:
dMu[mukey] = comppct_r
# look for specific strings in the map unit name
#if pw == False and muname is not None:
elif [d for d in phaseList if muname.find(d) >= 0]:
pw = True
#if mukey == tmukey:
# PrintMsg("\tMuname = " + muname, 1)
try:
sumPct = dMu[mukey]
dMu[mukey] = sumPct + comppct_r
except:
dMu[mukey] = comppct_r
elif str(drainagecl) in drainList:
pw = True
#if mukey == tmukey:
# PrintMsg("\tDrainageClass = " + drainagecl, 1)
try:
sumPct = dMu[mukey]
dMu[mukey] = sumPct + comppct_r
except:
dMu[mukey] = comppct_r
lastCokey = cokey # use this to skip the rest of the horizons for this component
arcpy.SetProgressorPosition()
if len(dMu) > 0:
arcpy.SetProgressor("step", "Populating " + os.path.basename(theMuTable) + "...", 0, len(dMu), 1)
# Populate the PWSL1POMU column in the map unit level table
muFlds = ["mukey", "pwsl1pomu"]
with arcpy.da.UpdateCursor(theMuTable, muFlds) as muCur:
for rec in muCur:
mukey = rec[0]
try:
rec[1] = dMu[mukey]
muCur.updateRow(rec)
except:
pass
arcpy.SetProgressorPosition()
arcpy.ResetProgressor()
return True
except MyError, e:
# Example: raise MyError("this is an error message")
PrintMsg(str(e) + " \n", 2)
return False
except:
errorMsg()
return False
## ===================================================================================
def StateNames():
# Create dictionary object containing list of state abbreviations and their names that
# will be used to name the file geodatabase.
# For some areas such as Puerto Rico, U.S. Virgin Islands, Pacific Islands Area the
# abbrevation is
# NEED TO UPDATE THIS FUNCTION TO USE THE LAOVERLAP TABLE AREANAME. AREASYMBOL IS STATE ABBREV
try:
stDict = dict()
stDict["AL"] = "Alabama"
stDict["AK"] = "Alaska"
stDict["AS"] = "American Samoa"
stDict["AZ"] = "Arizona"
stDict["AR"] = "Arkansas"
stDict["CA"] = "California"
stDict["CO"] = "Colorado"
stDict["CT"] = "Connecticut"
stDict["DC"] = "District of Columbia"
stDict["DE"] = "Delaware"
stDict["FL"] = "Florida"
stDict["GA"] = "Georgia"
stDict["HI"] = "Hawaii"
stDict["ID"] = "Idaho"
stDict["IL"] = "Illinois"
stDict["IN"] = "Indiana"
stDict["IA"] = "Iowa"
stDict["KS"] = "Kansas"
stDict["KY"] = "Kentucky"
stDict["LA"] = "Louisiana"
stDict["ME"] = "Maine"
stDict["MD"] = "Maryland"
stDict["MA"] = "Massachusetts"
stDict["MI"] = "Michigan"
stDict["MN"] = "Minnesota"
stDict["MS"] = "Mississippi"
stDict["MO"] = "Missouri"
stDict["MT"] = "Montana"
stDict["NE"] = "Nebraska"
stDict["NV"] = "Nevada"
stDict["NH"] = "New Hampshire"
stDict["NJ"] = "New Jersey"
stDict["NM"] = "New Mexico"
stDict["NY"] = "New York"
stDict["NC"] = "North Carolina"
stDict["ND"] = "North Dakota"
stDict["OH"] = "Ohio"
stDict["OK"] = "Oklahoma"
stDict["OR"] = "Oregon"
stDict["PA"] = "Pennsylvania"
stDict["PRUSVI"] = "Puerto Rico and U.S. Virgin Islands"
stDict["RI"] = "Rhode Island"
stDict["Sc"] = "South Carolina"
stDict["SD"] ="South Dakota"
stDict["TN"] = "Tennessee"
stDict["TX"] = "Texas"
stDict["UT"] = "Utah"
stDict["VT"] = "Vermont"
stDict["VA"] = "Virginia"
stDict["WA"] = "Washington"
stDict["WV"] = "West Virginia"
stDict["WI"] = "Wisconsin"
stDict["WY"] = "Wyoming"
return stDict
except:
PrintMsg("\tFailed to create list of state abbreviations (CreateStateList)", 2)
return stDict
## ===================================================================================
def UpdateMetadata(outputWS, target, surveyInfo):
# Update metadata for target object (VALU1 table)
#
try:
# Clear process steps from the VALU1 table. Mostly AddField statements.
#
# Different path for ArcGIS 10.2.2??
#
#
if not arcpy.Exists(target):
target = os.path.join(outputWS, target)
# Remove geoprocessing history
remove_gp_history_xslt = os.path.join(os.path.dirname(sys.argv[0]), "remove geoprocessing history.xslt")
out_xml = os.path.join(env.scratchFolder, "xxClean.xml")
if arcpy.Exists(out_xml):
arcpy.Delete_management(out_xml)
# Using the stylesheet, write 'clean' metadata to out_xml file and then import back in
arcpy.XSLTransform_conversion(target, remove_gp_history_xslt, out_xml, "")
arcpy.MetadataImporter_conversion(out_xml, os.path.join(outputWS, target))
# Set metadata translator file
dInstall = arcpy.GetInstallInfo()
installPath = dInstall["InstallDir"]
prod = r"Metadata/Translator/ARCGIS2FGDC.xml"
mdTranslator = os.path.join(installPath, prod)
# Define input and output XML files
#mdExport = os.path.join(env.scratchFolder, "xxExport.xml") # initial metadata exported from current data data
xmlPath = os.path.dirname(sys.argv[0])
mdExport = os.path.join(xmlPath, "gSSURGO_ValuTable.xml") # template metadata stored in ArcTool folder
mdImport = os.path.join(env.scratchFolder, "xxImport.xml") # the metadata xml that will provide the updated info
# Cleanup XML files from previous runs
if os.path.isfile(mdImport):
os.remove(mdImport)
# Start editing metadata using search and replace
#
stDict = StateNames()
st = os.path.basename(outputWS)[8:-4]
if st in stDict:
# Get state name from the geodatabase name
mdState = stDict[st]
else:
mdState = ""
# Set date strings for metadata, based upon today's date
#
d = datetime.now()
#today = d.strftime('%Y%m%d')
# Alternative to using today's date. Use the last SAVEREST date
today = GetLastDate(outputWS)
# Set fiscal year according to the current month. If run during January thru September,
# set it to the current calendar year. Otherwise set it to the next calendar year.
#
if d.month > 9:
fy = "FY" + str(d.year + 1)
else:
fy = "FY" + str(d.year)
# Convert XML from template metadata to tree format
tree = ET.parse(mdExport)
root = tree.getroot()
# new citeInfo has title.text, edition.text, serinfo/issue.text
citeInfo = root.findall('idinfo/citation/citeinfo/')
if not citeInfo is None:
# Process citation elements
# title
#
# PrintMsg("citeInfo with " + str(len(citeInfo)) + " elements : " + str(citeInfo), 1)
for child in citeInfo:
PrintMsg("\t\t" + str(child.tag), 0)
if child.tag == "title":
child.text = os.path.basename(target).title()
if mdState != "":
child.text = child.text + " - " + mdState
elif child.tag == "edition":
if child.text.find('xxFYxx') >= 0:
child.text = child.text.replace('xxFYxx', fy)
else:
PrintMsg(" \n\tEdition: " + child.text, 1)
if child.text.find('xxTODAYxx') >= 0:
child.text = child.text.replace('xxTODAYxx', today)
elif child.tag == "serinfo":
for subchild in child.iter('issue'):
if subchild.text == "xxFYxx":
subchild.text = fy
if child.text.find('xxTODAYxx') >= 0:
child.text = child.text.replace('xxTODAYxx', today)
# Update place keywords
#PrintMsg("\tplace keywords", 0)
ePlace = root.find('idinfo/keywords/theme')
if ePlace is not None:
for child in ePlace.iter('themekey'):
if child.text == "xxSTATExx":
#PrintMsg("\tReplaced xxSTATExx with " + mdState)
child.text = mdState
elif child.text == "xxSURVEYSxx":
#child.text = "The Survey List"
child.text = surveyInfo
else:
PrintMsg("\tsearchKeys not found", 1)
idDescript = root.find('idinfo/descript')
if not idDescript is None:
for child in idDescript.iter('supplinf'):
#id = child.text
#PrintMsg("\tip: " + ip, 1)
if child.text.find("xxTODAYxx") >= 0:
#PrintMsg("\t\tip", 1)
child.text = child.text.replace("xxTODAYxx", today)
if child.text.find("xxFYxx") >= 0:
#PrintMsg("\t\tip", 1)
child.text = child.text.replace("xxFYxx", fy)
if not idDescript is None:
for child in idDescript.iter('purpose'):
#ip = child.text
#PrintMsg("\tip: " + ip, 1)
if child.text.find("xxFYxx") >= 0:
#PrintMsg("\t\tip", 1)
child.text = child.text.replace("xxFYxx", fy)
if child.text.find("xxTODAYxx") >= 0:
#PrintMsg("\t\tip", 1)
child.text = child.text.replace("xxTODAYxx", today)
idAbstract = root.find('idinfo/descript/abstract')
if not idAbstract is None:
iab = idAbstract.text
if iab.find("xxFYxx") >= 0:
#PrintMsg("\t\tip", 1)
idAbstract.text = iab.replace("xxFYxx", fy)
#PrintMsg("\tAbstract", 0)
# Use contraints
#idConstr = root.find('idinfo/useconst')
#if not idConstr is None:
# iac = idConstr.text
#PrintMsg("\tip: " + ip, 1)
# if iac.find("xxFYxx") >= 0:
# idConstr.text = iac.replace("xxFYxx", fy)
# PrintMsg("\t\tUse Constraint: " + idConstr.text, 0)
# Update credits
eIdInfo = root.find('idinfo')
if not eIdInfo is None:
for child in eIdInfo.iter('datacred'):
sCreds = child.text
if sCreds.find("xxTODAYxx") >= 0:
#PrintMsg("\tdata credits1", 1)
sCreds = sCreds.replace("xxTODAYxx", today)
if sCreds.find("xxFYxx") >= 0:
#PrintMsg("\tdata credits2", 1)
sCreds = sCreds.replace("xxFYxx", fy)
child.text = sCreds
#PrintMsg("\tCredits: " + sCreds, 1)
# create new xml file which will be imported, thereby updating the table's metadata
tree.write(mdImport, encoding="utf-8", xml_declaration=None, default_namespace=None, method="xml")
# import updated metadata to the geodatabase table
arcpy.MetadataImporter_conversion(mdExport, target)
arcpy.ImportMetadata_conversion(mdImport, "FROM_FGDC", target, "DISABLED")
# delete the temporary xml metadata files
if os.path.isfile(mdImport):
os.remove(mdImport)
#if os.path.isfile(mdExport):
# os.remove(mdExport)
return True
except:
errorMsg()
False
## ===================================================================================
## ====================================== Main Body ==================================
# Import modules
import os, sys, string, re, locale, arcpy, traceback, collections
from operator import itemgetter, attrgetter
import xml.etree.cElementTree as ET
from datetime import datetime
from arcpy import env
# Original input table fields:
# areasymbol, mukey, musym, muname, muname, cokey, compct, compname, compkind, localphase,
# taxorder, taxsubgrp, ec, pH, dbthirdbar, hzname, hzdesgn, hzdept, hzdepb, hzthk, sand, silt,
# clay, om, reskind, reshard, resdept, resthk, texture, lieutex
# Compkind values for major components: Miscellaneous area, Series, Taxon above family,
# Family, Taxadjunct, Variant
# chtexturegrp.texture: 3,468 unique values, too many to list here. Bob's queries only
# look at excluding 'COP-MAT' texture from the Dense Layer calculation. He identifies the rest
# of the organic horizons using lieutext.
# lieutext values: Slightly decomposed plant material, Moderately decomposed plant material,
# Bedrock, Variable, Peat, Material, Unweathered bedrock, Sand and gravel, Mucky peat, Muck,
# Highly decomposed plant material, Weathered bedrock, Cemented, Gravel, Water, Cobbles,
# Stones, Channers, Parachanners, Indurated, Cinders, Duripan, Fragmental material, Paragravel,
# Artifacts, Boulders, Marl, Flagstones, Coprogenous earth, Ashy, Gypsiferous material,
# Petrocalcic, Paracobbles, Diatomaceous earth, Fine gypsum material, Undecomposed organic matter
#
# reskind values: Strongly contrasting textural stratification, Lithic bedrock, Densic material,
# Ortstein, Permafrost, Paralithic bedrock, Cemented horizon, Undefined, Fragipan, Plinthite,
# Abrupt textural change, Natric, Petrocalcic, Duripan, Densic bedrock, Salic,
# Human-manufactured materials, Sulfuric, Placic, Petroferric, Petrogypsic
# Dobos - NASIS lieutext matches: mpm, mpt, muck, peat, spm, udom, pdom, hpm
#
# Paul's organic horizon filters: chtexturegrp.texture <> 'SPM', <>'UDOM'???, NOT LIKE: '%MPT%', '%MUCK', '%PEAT%' (from SVI query)
#
# Desgnmaster filter: hzdesn IN ["O", "O'", "L"]
#
# Taxonomic Order filter: upper(taxorder) LIKE 'HISTOSOLS'. Is this being used?
# Perhaps I should be using taxsubgrp.lower() like '%histic%' or use both!!!
#
# compkind filter for earthy components: <> 'Miscellaneous area'
try:
# Diagnostic map unit mukey:
tmukey = 'xx'
# Create geoprocessor object
#gp = arcgisscripting.create(9.3)
arcpy.OverwriteOutput = True
inputDB = arcpy.GetParameterAsText(0) # Input gSSURGO database
# Set location for temporary tables
#outputDB = "IN_MEMORY"
outputDB = env.scratchGDB
# Name of mapunit level output table (global variable)
theMuTable = os.path.join(inputDB, "Valu2")
# Name of component level output table (global variable)
theCompTable = os.path.join(inputDB, "Co_VALU")
# Set output workspace to same as the input table
#env.workspace = os.path.dirname(arcpy.Describe(queryTbl).catalogPath)
env.workspace = inputDB
# Save record of any issues to a text file
logFile = os.path.basename(inputDB)[:-4] + "_Problems.txt"
logFile = os.path.join(os.path.dirname(inputDB), logFile)
# Get the mapunit - sum of component percent for calculations
dPct = GetSumPct(inputDB)
if len(dPct) == 0:
raise MyError, ""
# Create initial set of query tables used for RZAWS, AWS and SOC
if CreateQueryTables(inputDB, outputDB, 150.0) == False:
raise MyError, ""
# Create permanent output tables for the map unit and component levels
depthList = [(0,5), (5, 20), (20, 50), (50, 100), (100, 150), (150, 999), (0, 20), (0, 30), (0, 100), (0, 150), (0, 999)]
if CreateOutputTableMu(theMuTable, depthList, dPct) == False:
raise MyError, ""
if CreateOutputTableCo(theCompTable, depthList) == False:
raise MyError, ""
# Store component restrictions for root growth in a dictionary
resListAWS = "('Lithic bedrock','Paralithic bedrock','Densic bedrock', 'Densic material', 'Fragipan', 'Duripan', 'Sulfuric')"
dRZRestrictions = GetCoRestrictions(outputDB, 150.0, resListAWS)
# Find the top restriction for each component, both from the corestrictions table and the horizon properties
dComp2 = CalcRZDepth(inputDB, outputDB, theCompTable, theMuTable, 150.0, dPct, dRZRestrictions)
# Calculate root zone available water capacity using a floor of 150cm or a root restriction depth
#
# dComp2[cokey] = [mukey, compName, localPhase, compPct, resDept, restriction]
if CalcRZAWS(inputDB, outputDB, 0.0, 150.0, theCompTable, theMuTable, dComp2, 150.0, dPct) == False:
raise MyError, ""
# Calculate standard available water supply
if CalcAWS(inputDB, outputDB, theCompTable, theMuTable, dPct) == False:
raise MyError, ""
# Run SOC calculations
# Seems to be a problem with SOC calculations, numbers are high
maxD = 999.0
# Get bedrock restrictions for SOC and write them to the output tables
resListSOC = "('Lithic bedrock', 'Paralithic bedrock', 'Densic bedrock')"
dSOCRestrictions = GetCoRestrictions(outputDB, maxD, resListSOC)
# Store horizon fragment volumes in a dictionary and use in the root zone SOC calculations
dFrags = GetFragVol(inputDB)
if len(dFrags) == 0:
raise MyError, "No fragment volume information"
# Calculate soil organic carbon for all the different depth ranges
depthList = [(0,5), (5, 20), (20, 50), (50, 100), (100, 150), (150, 999), (0, 20), (0, 30), (0, 100), (0, 150), (0, 999)]
if CalcSOC(inputDB, outputDB, theCompTable, theMuTable, dPct, dFrags, depthList, dSOCRestrictions, maxD) == False:
raise MyError, ""
# Calculate NCCPI
# Create query table using component and chorizon tables
nccpiTbl = "NCCPI_Table"
if MakeNCCPIQueryTable(inputDB, nccpiTbl) == False:
raise MyError, ""
if CalcNCCPI(inputDB, nccpiTbl, dPct) == False:
raise MyError, ""
if CalcPWSL(inputDB, outputDB, theMuTable, dPct) == False:
raise MyError, ""
PrintMsg(" \nAll calculations complete", 0)
# Create metadata for the VALU table
# Query the output SACATALOG table to get list of surveys that were exported to the gSSURGO
#
saTbl = os.path.join(inputDB, "sacatalog")
expList = list()
queryList = list()
with arcpy.da.SearchCursor(saTbl, ["AREASYMBOL", "SAVEREST"]) as srcCursor:
for rec in srcCursor:
expList.append(rec[0] + " (" + str(rec[1]).split()[0] + ")")
queryList.append("'" + rec[0] + "'")
surveyInfo = ", ".join(expList)
queryInfo = ", ".join(queryList)
# Update metadata for the geodatabase and all featureclasses
PrintMsg(" \nUpdating " + os.path.basename(theMuTable) + " metadata...", 0)
bMetadata = UpdateMetadata(inputDB, theMuTable, surveyInfo)
PrintMsg("\tMetadata complete", 0)
PrintMsg(" \nValu2 table complete for " + inputDB + " \n ", 0)
except MyError, e:
# Example: raise MyError("this is an error message")
PrintMsg(str(e) + " \n", 2)
except:
errorMsg()
