Select tree species basal area raster data were retrieved and summed into 1 raster using Map Algebra. The following python code was used to estimate an average basal area for select tree species per NHDPlus reach catchment:
import os, sys, arcpy
from arcpy import env
from arcpy import sa
from arcpy.sa import *
from datetime import datetime
wks = r"E:/basal/NLCD4"
env.workspace = wks
startTime = datetime.now()
# Check out Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")
from arcpy.sa import *
from arcpy import env
arcpy.env.overwriteOutput = True
env.qualifiedFieldNames = False
# Set the Geoprocessing environment...
arcpy.scratchWorkspace = "E:/basal/NLCD4"
arcpy.env.geographicTransformations = ""
arcpy.env.outputCoordinateSystem = "PROJCS['NAD_1983_Albers',GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Albers'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-96.0],PARAMETER['Standard_Parallel_1',29.5],PARAMETER['Standard_Parallel_2',45.5],PARAMETER['Latitude_Of_Origin',23.0],UNIT['Meter',1.0]]"
arcpy.env.pyramid = "NONE"
NHD_dir = "Q:/WORK/NHDPlusV2"
tv = "tv"
tvStat = "tvStat"
tmpTable2 = "temp2.dbf"
#Local variables...
inputs = {'CA':['18'],'CO':['14','15'],'GB':['16'],'GL':['04'],'MA':['02'],'MS':['05','06','07','08','10L','10U','11'],'NE':['01'],'PN':['17'],\
'RG':['13'],'SA':['03N','03S','03W'],'SR':['09'],'TX':['12']}
for regions in inputs.keys():
for hydro in inputs[regions]:
print ('Region: ' + regions + ' Production unit: ' + hydro)
inZoneData = "%s/NHDPlus%s/NHDPlus%s/NHDPlusCatchment/cat"%(NHD_dir,regions,hydro)
print (inZoneData)
zoneField = "VALUE"
print (zoneField)
valGrid = r"E:/basal/NLCD4/NLCD4.gdb/Nbasal_area"
print (valGrid)
outTable = "temp0.dbf"
print (outTable)
arcpy.env.cellSize = inZoneData
arcpy.env.Extent = inZoneData
arcpy.env.SnapRaster = inZoneData
outTable1 = "E:/basal/dbfs/Basal_PU%s.dbf"%(hydro)
outTable1 = os.path.join(wks,outTable1)
catShape= "%s/NHDPlus%s/NHDPlus%s/NHDPlusCatchment/catchment.shp"%(NHD_dir,regions,hydro)
print (catShape)
print (outTable1)
outZSaT = ZonalStatisticsAsTable(inZoneData, zoneField, valGrid, outTable, "DATA", "MEAN")
# join catchment shape attributes to stats and write out
arcpy.MakeTableView_management(catShape,tv)
arcpy.MakeTableView_management(outTable,tvStat)
#arcpy.AddIndex_management(tvStat, "VALUE","IDX")
arcpy.AddJoin_management (tv, "GRIDCODE", tvStat, "VALUE")
arcpy.CopyRows_management(tv,tmpTable2)
print (catShape)
# Set up desired fields and order for output
arcpy.CreateTable_management(os.path.dirname(outTable1), os.path.basename(outTable1))
arcpy.AddField_management(outTable1,"FEATUREID","LONG")
arcpy.AddField_management(outTable1,"AREASQKM","DOUBLE")
arcpy.AddField_management(outTable1,"basin_area","DOUBLE")
arcpy.AddField_management(outTable1,"nodata_are","DOUBLE")
arcpy.AddField_management(outTable1,"NODATA","DOUBLE")
arcpy.AddField_management(outTable1,"AREA","DOUBLE")
arcpy.AddField_management(outTable1,"MEAN","DOUBLE") # statistic value
arcpy.AddField_management(outTable1,"BASAL","DOUBLE") # raster name #statistic value
arcpy.AddField_management(outTable1,"COMID","LONG")
# delete dummy field added by CreateTable
arcpy.DeleteField_management(outTable1,"ID")
# Copy matching fields to the output table
arcpy.Append_management(tmpTable2,outTable1,"NO_TEST")
# tag non-matches with the zonal stat = -9999
arcpy.MakeTableView_management(outTable1,tv)
arcpy.CalculateField_management(tv,'COMID', '!FEATUREID!', "PYTHON_9.3")
arcpy.CalculateField_management(tv,'basin_area', '!AREASQKM! * 1000000', "PYTHON_9.3")
arcpy.CalculateField_management(tv,'BASAL', '!MEAN!', "PYTHON_9.3")
arcpy.CalculateField_management(tv,'nodata_are', '!basin_area! - !AREA!', "PYTHON_9.3")
arcpy.CalculateField_management(tv,'NODATA', '(!nodata_are! /!basin_area!) * 100', "PYTHON_9.3")
where = "AREA = 0" # if output is GDB use "AREA IS NULL"
arcpy.SelectLayerByAttribute_management(tv,"",where)
arcpy.CalculateField_management(tv,"NODATA","100", "PYTHON_9.3")
arcpy.CalculateField_management(tv,"BASAL","0", "PYTHON_9.3")
