'Name: NHDBASIN.ProjectCalculateBasinCharacteristics
'Purpose: Calculate the project specific basin charactistics
'Author: Jen Hill, Horizon Systems Corporation
'Date: March 2002
'Input Parameters: None
'Input Files: None
'Output Parameters (Return Values): None
'Output Files: None
'Calls: None
'Called By:  NHDBASIN.btnBasinCharacteristicsClick
'Globals: None
'Locals: 
'Notes: None

'------------------------------------------------------------------------
'Revision Author: Peter Steeves, US Geological Survey
'Revision Date: 9/30/02
'Revision Notes: The 5 sample characteristic summaries below are a diverse
' range of calculation types to give users a few examples for calculating their
' own characteristics.  Each characteristic summary is headed by a text section
' explaining what is being calculated, and how the user could modify the characteristic
' for their own project.  
' These characteristics were the initial set used in Massachusetts as variables for a regression
' equation that determines whether a stream is perennial or intermittent.
'------------------------------------------------------------------------
'
'------------------------------------------------------------------------
'Revision Author:
'Revision Date:
'Revision Notes:
'------------------------------------------------------------------------
'


av.ShowMsg("Calculating basin characteristics")

viewCurrent   = av.GetActiveDoc
strWatershed  = "Watershed Tool Results"
listReach     =  av.run("NHDZCOMMON.RetrieveValue","NHDWSSELECTEDREACH")
'listReach contains:
'  0   record number
'  1   com_id
'  2   rch_code
'  3   measure
'  4   theme name
                      
'Gather some information that will be necessary for several of the basin charactistics                      

'Get information about the delineated watershed 
   themeWS = viewCurrent.FindTheme(strWatershed)
   vtabWS  = themeWS.GetFtab
   'There will be only one record
   for each objRec in vtabWS
      'Get the sq_km 
      numWSSq_km = vtabWS.ReturnValue(vtabWS.FindField("sq_km"),objRec)
      'Get the shape
      shpWS = vtabWS.ReturnValue(vtabWS.FindField("shape"),objRec)
   end
   'msgbox.info("1","")
   
'Get information about the location of the basin characteristics supporting data location
   fnamePath = av.Run("NHDZCOMMON.RetrieveValue","NHDBASINDATALOCATION")
   if (fnamePath = Nil) then
      return nil
   end
   'msgbox.info("2","")

'Set the analysis environment  
   ae=viewCurrent.GetExtension(AnalysisEnvironment)
   ae.setExtent(#ANALYSISENV_VALUE,shpWS.ReturnExtent)
   ae.Activate
   Grid.SetAnalysisCellSize(#GRID_ENVTYPE_VALUE,10)
   'msgbox.info("3","")

'1.  Area in square miles

' By default NHD Watershed should be calculating basin area in square 
' kilometers.  This example shows a simple conversion of that value to square miles.
' The 'drainageArea' is the only line of code in this section, and would need to be 
' modified for other area measurements (i.e. square feet).


   'Convert from sq_km to sq miles
   'To convert sq_km into sq miles multiply sq_km by .3861

   drainageArea = numWSSq_km * 0.3861 'conversion to square miles
   'msgbox.info("4","")

'2. Area of Stratified Drift
'   This data is in GRID format.  drift = 2 / other = 1

   'Create theme and grid objects for the  surficial geology grid
   srcnamesg = grid.MakeSrcName(fnamePath.AsString +"\sg_gr")
   themesg = theme.make( srcnamesg )
   gridsg = themesg.GetGrid
   'viewCurrent.AddTheme(themesg)
   'msgbox.info("4a","")
   
   'Get the intersection of the watershed polygon and the grid
   gridPartial = gridsg.ExtractByPolygon(shpWS, prj.MakeNull, false)  
   themePartial = gtheme.make(gridPartial)
   'viewCurrent.AddTheme(themePartial)
   'msgbox.info("4b","")
   
   'Get the counts for values for 1 and 2
   ftabPartial = themePartial.GetVtab
   num1 = 0
   num2 = 0
   for each objRec in ftabPartial
      if (ftabPartial.ReturnValue(ftabPartial.FindField("value"),objRec) = 1) then
         num1 = ftabPartial.ReturnValue(ftabPartial.FindField("count"),objRec)
      end
      if (ftabPartial.ReturnValue(ftabPartial.FindField("value"),objRec) = 2) then
         num2 = ftabPartial.ReturnValue(ftabPartial.FindField("count"),objRec)
      end
   end
      
   driftPercent = (num2 / (num1 + num2) ) * 100
   'msgbox.info("5","")

   

'3. Drainage Density   '   

'   Get the length from the NHD Navigation tool summary of stream length
'   Convert length to miles (the perennial/intermittent project data is in units meters)
'   Get basin area from the NHD Watershed tool summary table

'JRH - Using NHD Navigate to get the upstream mainstem distance 
'listNHDNAVNavOptions = list.Make
'listNHDNAVNavOptions.Add(nil)     'No maximum travel distance
'listNHDNAVNavOptions.Add("MI")    'Distances in miles
'listNHDNAVNavOptions.Add("L")     'Level Path
'listNHDNAVNavOptions.Add("P")     'Point on reach start mode
'listNHDNAVNavOptions.Add(False)   'Do not stop at WS boundary
'listReach = av.run("NHDZCOMMON.RetrieveValue","NHDWSSELECTEDREACH")
'listNHDNAVSelectedPoints = listReach
''To display or not display the navigation results, use the flag in the last parameter  
'av.run("NHDNAV.CoordinateNavigation",{listNHDNAVSelectedPoints,"U", 1, listNHDNAVNavOptions, false} )

strOutFile = "output.txt"
strWorkDir = system.GetEnvVar("TEMP")+"\"

'Create a filename object for the navigation results (using the base name and view name)
fnameResults = ( strWorkDir+"\"+ viewCurrent.GetName.Substitute(" ","_") + strOutFile).asfilename
'Create the vtab object for the navigation results
vtabResults=vtab.make (fnameResults, false, false)   
numMeters = 0
for each objRec in vtabResults
   numRec = objRec
   'Meters if using the results from watershed delineation
   numMeters = numMeters + vtabResults.ReturnValue(vtabResults.FindField("rchlength"),numRec)
end
'From upstream MS  - when navigation occurs using MI as the units
'numMiles = vtabResults.ReturnValue(vtabResults.FindField("distance"),numRec)
'Convert meters to miles
numMiles = numMeters / 1609
drainageDensity = numMiles / drainageArea
    'msgbox.info("6","")

 
    
'4. Mean Basin Elevation
'   This data is in GRID format.

   'Create theme and grid objects for the  elevation grid
   srcnameelev = grid.MakeSrcName(fnamePath.AsString +"\ned_gr")
   themeelev = theme.make( srcnameelev )
   gridelev = themeelev.GetGrid
   'viewCurrent.AddTheme(themeelev)
   
   'Get the intersection of the watershed polygon and the grid
   gridPartial = gridelev.ExtractByPolygon(shpWS, prj.MakeNull, false)  
   themePartial = gtheme.make(gridPartial)
   'viewCurrent.AddTheme(themePartial)
   
   listStatistics = gridPartial.GetStatistics
   
   meanBasinElevation_meters = listStatistics.Get(2)
   meanBasinElevation = meanBasinElevation_meters * 3.2808
   'msgbox.info("7","")
   
'5. Area of Forest
'   
'   Any Forest dataset (or other landcover) can be used once it is in GRID format
'   This data is in GRID format. forest = 2 / other = 1

   'Create theme and grid objects for the  forest grid
   srcnameforest = grid.MakeSrcName(fnamePath.AsString +"\forest_gr")
   themeforest = theme.make( srcnameforest )
   gridforest = themeforest.GetGrid
   'viewCurrent.AddTheme(themeforest)
   
   'Get the intersection of the watershed polygon and the grid
   gridPartial = gridforest.ExtractByPolygon(shpWS, prj.MakeNull, false)  
   themePartial = gtheme.make(gridPartial)
   'viewCurrent.AddTheme(themePartial)
   
      'Get the counts for values for 1 and 2
   ftabPartial = themePartial.GetVtab
   num1 = 0
   num2 = 0
   for each objRec in ftabPartial
      if (ftabPartial.ReturnValue(ftabPartial.FindField("value"),objRec) = 1) then
         num1 = ftabPartial.ReturnValue(ftabPartial.FindField("count"),objRec)
      end
      if (ftabPartial.ReturnValue(ftabPartial.FindField("value"),objRec) = 2) then
         num2 = ftabPartial.ReturnValue(ftabPartial.FindField("count"),objRec)
      end
   end
      
   forestPercent = (num2 / (num1 + num2) ) * 100
   'msgbox.info("8","")


listBasinCharacteristics = list.Make
listBasinCharacteristics.Add(drainageArea)
listBasinCharacteristics.Add(driftPercent)
listBasinCharacteristics.Add(drainageDensity)
listBasinCharacteristics.Add(meanBasinElevation)
listBasinCharacteristics.Add(forestPercent)

   
' for all lines below, switch commented lines to uncommented for looping, and visa versa
' for manual method.  Also, make the switch for one line at the bottom of PI.CalcProbability ...
' (msgbox.info(z.asstring,"")  
msgbox.info("Square Miles: "       + drainageArea.AsString      + NL +
            "Percent Stratified Drift: " + driftPercent.AsString     + NL +
            "Drainage Density:  "       + drainageDensity.AsString     + NL +
            "Mean Basin Elevation: "           + meanBasinElevation.AsString      + NL +
            "Percent Forest: " + forestPercent.AsString,"")      


av.ShowMsg("")

'return listBasinCharacteristics


passlist = {drainageArea, driftPercent, drainageDensity, meanBasinElevation, forestPercent}
av.Run ("PI.CalcProbability", passlist)