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
