Mean slope calculation code for VBScript or Python

[TysonHart]

I've been working tirelessly on trying to develop a code that will return an average slope for a selected area.  I can build a model to give me the answer, but not as concise as I want it.

The ideal code would be VBScript in which the user select an existing shapfile, runs the code, then a message box appears stating "Average Slope = XX.XX%".

I have an elevation database on C:\geodata\elevation\ containing a 10m DEM and a slope file derived using spatial analyst of that DEM.  The problem I see is using the extent/mask of a selected shapefile to calculate the mean pixel value of the slope raster.

We use a code for finding shape acres (code listed below) and I've been trying to piggyback off of it with no avail.  I've also tried using python but am having trouble making a scipt that will list a message box stating "Average Slope =".  When I use this method I use zonal statistics and it provides a nice large table, but I only want the Mean value to be printed out as a message box.

Thanks in advance for any advice anyone may have,

- Tyson

Public Sub ShapeAcres()

Dim pMxDoc As IMxDocument
Dim pArea As IArea
Dim pAcres As Double, pFeet As Double, pMeters As Double
Dim i As Integer, Looper As Integer
Dim pPoly As IFeature
Dim pMap As IMap
Dim pFSelection As IEnumFeature
Dim pActiveView As IActiveView
Dim pContentsView As IContentsView

Set pMxDoc = ThisDocument
Set pMap = pMxDoc.FocusMap
Set pActiveView = pMap
Set pContentsView = pMxDoc.CurrentContentsView
Set pFSelection = pMap.FeatureSelection
i = pMap.SelectionCount

'-- Make sure element is selected
If i = 0 Then
    MsgBox "Please select one or more shapes"
    Exit Sub
End If

Do Until Looper = i
    Set pPoly = pFSelection.Next
    Looper = Looper + 1
    Set pArea = pPoly.Shape
    pMeters = (pMeters + pArea.Area)
Loop

pFeet = pMeters * 10.76391042
pAcres = pFeet / 43560
MsgBox "Percent Slope = " & Format(pAcres, "##,###.0")

End Sub

[JeremySury]

This may not be what you're looking for and its certainly not what I would call concise however the attached code generates a trend plane through all the points within each drainage area.  It then calculates the slope and aspect of the trend plane and enters them into the drainage area feature.

I used point geometry as the input to the IInterpolationOp.Trend method but I believe you may be able to use a grid instead.

Cheers!

Public Sub CreateTrendPlane()
Dim pMxDoc As IMxDocument
Dim pMap As IMap
Dim pSF As ISpatialFilter
Dim pLiDARPointsLayer  As IFeatureLayer
Dim pDALayer As IFeatureLayer
Dim pDAFeatCur As IFeatureCursor
Dim pDAFeat As IFeature
Dim pSelectionSet As ISelectionSet
Dim pArea As IArea
Dim pPoint As IPoint
Dim pArray As IArray
Dim pRasterIdentify As IRasterIdentifyObj
Dim pSurfaceOp As ISurfaceOp
Dim pRasterLayer As IRasterLayer
Dim pOutputRaster As IRaster
Dim pAspectRaster As IRaster
Dim pSlopeRaster As IRaster
Dim pIdentify As IIdentify
Dim lAspectFieldIndex As Long
Dim lSlopeFieldIndex As Long

Const ElevationFieldName = "GRID_CODE"

Set pMxDoc = ThisDocument
Set pMap = pMxDoc.FocusMap

Set pDALayer = pMap.Layer(0) ' put DA layer at top of TOC
Set pLiDARPointsLayer = pMap.Layer(1) ' put LidarPoints layer next

lAspectFieldIndex = pDALayer.FeatureClass.FindField("MyAspect")
lSlopeFieldIndex = pDALayer.FeatureClass.FindField("MySlope")

Dim iCount As Integer
Dim iNumber As Integer
iCount = 1
iNumber = pDALayer.FeatureClass.FeatureCount(Nothing)

Set pDAFeatCur = pDALayer.FeatureClass.Search(Nothing, False)
Set pDAFeat = pDAFeatCur.NextFeature

While Not pDAFeat Is Nothing

    Set pSF = New SpatialFilter
    Set pSF.Geometry = pDAFeat.Shape
    pSF.GeometryField = pLiDARPointsLayer.FeatureClass.ShapeFieldName
    pSF.SpatialRel = esriSpatialRelContains
    Set pSelectionSet = pLiDARPointsLayer.FeatureClass.Select(pSF, esriSelectionTypeIDSet, esriSelectionOptionNormal, Nothing)
        
    Dim pFDescr As IFeatureClassDescriptor
    Set pFDescr = New FeatureClassDescriptor
    pFDescr.CreateFromSelectionSet pSelectionSet, Nothing, ElevationFieldName
    
    Dim pInterpolationOp As IInterpolationOp
    Set pInterpolationOp = New RasterInterpolationOp
    
    Dim pRasterAnalysisEnv As IRasterAnalysisEnvironment
    Set pRasterAnalysisEnv = pInterpolationOp
      
    pRasterAnalysisEnv.SetCellSize esriRasterEnvValue, CVar(15)
    pRasterAnalysisEnv.SetExtent esriRasterEnvMaxOf
    
    Dim pInputPoints As IGeoDataset
    Set pInputPoints = pFDescr
    
    Set pOutputRaster = pInterpolationOp.Trend(pInputPoints, esriGeoAnalysisLinearTrend, 1)
    
    Set pSurfaceOp = New RasterSurfaceOp
    
    Set pAspectRaster = pSurfaceOp.Aspect(pOutputRaster)
    Set pSlopeRaster = pSurfaceOp.Slope(pOutputRaster, esriGeoAnalysisSlopePercentrise)
    
    Set pRasterLayer = New rasterLayer
    pRasterLayer.CreateFromRaster pAspectRaster
        
    Set pIdentify = pRasterLayer
    Set pArea = pDAFeat.Shape
    Set pPoint = pArea.Centroid
    Set pArray = pIdentify.Identify(pPoint)
    Set pRasterIdentify = pArray.Element(0)
  
    pDAFeat.Value(lAspectFieldIndex) = pRasterIdentify.Name
    pDAFeat.Store
        
    Set pRasterLayer = New rasterLayer
    pRasterLayer.CreateFromRaster pSlopeRaster
        
    Set pIdentify = pRasterLayer
    Set pArray = pIdentify.Identify(pPoint)
    Set pRasterIdentify = pArray.Element(0)
  
    pDAFeat.Value(lSlopeFieldIndex) = pRasterIdentify.Name
    pDAFeat.Store
        
    Set pDAFeat = pDAFeatCur.NextFeature
    
    ThisDocument.Parent.StatusBar.Message(0) = "Calculated " & CStr(iCount) & " of " & CStr(iNumber)
    iCount = iCount + 1
    
Wend

MsgBox "Done"

End Sub

[TysonHart]

Thanks for the reply.

I guess with my limited scripting knowledge I couldn't get the code to run.  There were several arguments that i wasn't able to debug.  But, in looking around I think I've found the code I want but I can't figure out how to get it to run.  Its probably very basic but I'm really struggling with making it run.  Here's the code, I just need to figure out how to implement it and adjust the resulting min/max to mean from the dbf table:

Public Sub GetLayerZonalMinMax(pRasLayer As IRasterLayer, ByRef min As Double, ByRef max As Double)
    If pRasLayer.Raster Is Nothing Then Exit Sub
    
    Dim pMxDoc As IMxDocument
    Dim pMap As IMap
    Dim pCompositeLayer As ICompositeLayer
    Dim pLayer As ILayer
    Set pMxDoc = ThisDocument
    Set pMap = pMxDoc.FocusMap
    Set pCompositeLayer = pMap.Layer(0)
    Set pLayer = pCompositeLayer.Layer(0)
    Dim pFeaturelayer As IFeatureLayer
    Set pFeaturelayer = pLayer
    Dim pDataset As IDataset
    Set pDataset = pFeaturelayer.FeatureClass
    Dim sPath As String
    sPath = pDataset.Workspace.pathname & "\" & pDataset.Name & ".shp"
    
    Dim outputTable As String
    Dim tablePath As String
    tablePath = "C:\gis\temp"
    Dim tableName As String
    tableName = "output.dbf"
    'tableName = pRasLayer.Name + ".dbf"
    outputTable = tablePath + "\" + tableName
    Dim rasPath As String
    rasPath = pRasLayer.FilePath
    
    Dim gp As Object
    Set gp = CreateObject("esriGeoprocessing.GPDispatch.1")
    gp.CheckOutExtension "spatial"
    
    ' Load required toolboxes...
    gp.AddToolbox "C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Spatial Analyst Tools.tbx"
    gp.ZonalStatisticsAsTable_sa sPath, "ZONE", rasPath, outputTable, "DATA"
    
    Dim pWSF As IWorkspaceFactory
    Dim pFeatureWS As IFeatureWorkspace
    Dim pTable As ITable
    
    Set pWSF = New ShapefileWorkspaceFactory
    Set pFeatureWS = pWSF.OpenFromFile(tablePath, 0)
    Set pTable = pFeatureWS.OpenTable(tableName)
    Dim myMin As Double
    Dim myMax As Double
    myMin = pTable.GetRow(0).Value(4)
    myMax = pTable.GetRow(0).Value(5)
    Set pDataset = pTable
    pDataset.Delete
    
End Sub