jwsury12

Creating a trend plane from a raster layer

Discussion created by jwsury12 on Dec 14, 2011
Is there any way to use a RasterLayer as the input to the IInterpolationOp.Trend function.  I can access my RasterLayer through the IGeoDataset Interface however the Trend function doesn't seem to like what it gets and throws an "ERROR 010085: SPOT is an invalid field name".  All of the help material suggests that the function expects a FeatureClass or a FeatureClassDescriptor.

I have added my attempts below.  The first works and generates a trend plane from LiDAR XY points.  The points are very cumbersome to work with and I was hoping in the second attempt (which doesn't work) to use the interpolated raster as my input elevations.

Your wise insights would be appreciated!

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

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

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

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 pSlopeRaster = pSurfaceOp.Slope(pOutputRaster, esriGeoAnalysisSlopePercentrise)
    
    Set pArea = pDAFeat.Shape
    Set pPoint = pArea.Centroid
        
    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

Public Sub CreateTrendPlaneUsingLiDAR()
        Dim pMxDoc As IMxDocument
        Dim pMap As IMap
        Dim pDALayer As IFeatureLayer
        Dim pLiDARLayer As IRasterLayer
        Dim pRasterGDS As IGeoDataset
        Dim lAspectFieldIndex As Long
        Dim lSlopeFieldIndex As Long
        Dim iCount As Integer
        Dim iNumber As Integer
        Dim pDAFeatCur As IFeatureCursor
        Dim pDAFeat As IFeature
        Dim pInterpolationOp As IInterpolationOp
        Dim pRasterAnalysisEnv As IRasterAnalysisEnvironment
        Dim pOutputRaster As IRaster
        Dim pSurfaceOp As ISurfaceOp
        Dim pAspectRaster As IRaster
        Dim pSlopeRaster As IRaster
        Dim pArea As IArea
        Dim pPoint As IPoint
        
        Set pMxDoc = ThisDocument
        Set pMap = pMxDoc.FocusMap

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

        Set pRasterGDS = pLiDARLayer

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

        iCount = 1
        iNumber = pDALayer.FeatureClass.FeatureCount(Nothing)

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

        While Not pDAFeat Is Nothing

            Set pInterpolationOp = New RasterInterpolationOp
            
            Set pRasterAnalysisEnv = pInterpolationOp

            pRasterAnalysisEnv.SetCellSize esriRasterEnvValue, 15
            pRasterAnalysisEnv.SetExtent esriRasterEnvValue, pDAFeat.Shape.Envelope

            Set pOutputRaster = pInterpolationOp.Trend(pRasterGDS, esriGeoAnalysisLinearTrend, 1)

            Set pSurfaceOp = New RasterSurfaceOp

            Set pAspectRaster = pSurfaceOp.Aspect(pOutputRaster)
            Set pSlopeRaster = pSurfaceOp.Slope(pOutputRaster, esriGeoAnalysisSlopePercentrise)

            Set pArea = pDAFeat.Shape
            Set pPoint = pArea.Centroid

            pDAFeat.Value(lAspectFieldIndex) = IdentifyPixelValue(pAspectRaster, pPoint.X, pPoint.Y)
            pDAFeat.Store

            pDAFeat.Value(lSlopeFieldIndex) = IdentifyPixelValue(pSlopeRaster, pPoint.X, pPoint.Y)
            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

Outcomes