|
' Dim pGeometry As IGeometry
' Dim pacGeo As IGeoFeatureLayer
' Set pacGeo = m_pActiveLayer
' If pacGeo.FeatureClass.ShapeType = esriGeometryPolygon Then
Set pGeometry = Me.MapControl1.TrackPolygon
' End If
If pGeometry Is Nothing Then Exit Sub
''''''''''''''''''''''''''''''''''''''''''''''''
'/////////////1.求出交集,并保存_HIS/////////////////////////////////////
Dim pInputTable As ITable
Set pInputTable = m_pActiveLayer
Dim pOverlayTable As ITable
Dim pTmpMapLyr As IFeatureLayer
Set pTmpMapLyr = New FeatureLayer
Dim wks As IFeatureWorkspace
Dim wsf As IWorkspaceFactory
Set wsf = New ShapefileWorkspaceFactory
'wsf.Create "c:", "tmp", Nothing, 0
Set wks = wsf.OpenFromFile("c:tmp", 0)
Set pTmpMapLyr.SpatialReference = m_pMap.SpatialReference
Set pTmpMapLyr.FeatureClass = wks.OpenFeatureClass("test1")
' Dim lyfc As IFeatureClass
' Set lyfc = pTmpMapLyr.FeatureClass
Dim pOutBuff As IFeatureBuffer
Dim pOutFeat As IFeature
' pOutFeat.FeatureType = esriFTSimple
Dim pOutCurTmp As IFeatureCursor
Set pOutCurTmp = pTmpMapLyr.FeatureClass.Insert(True)
Set pOutBuff = pTmpMapLyr.FeatureClass.CreateFeatureBuffer
Set pOutFeat = pOutBuff
Set pOutFeat.Shape = pGeometry
pOutCurTmp.InsertFeature pOutBuff
' Set pOverlayTable = pGeometry
Set pOverlayTable = pTmpMapLyr
' Error checking
If pInputTable Is Nothing Then
MsgBox "Table QI failed"
Exit Sub
End If
If pOverlayTable Is Nothing Then
MsgBox "Table QI failed"
Exit Sub
End If
Set wks = Nothing
' Set output location and feature class name
Dim pNewWSName As IWorkspaceName
Set pNewWSName = New WorkspaceName
pNewWSName.WorkspaceFactoryProgID = "esriCore.ShapeFileWorkspaceFactory.1"
' Define the output feature class name
Dim pFeatClassName As IFeatureClassName
Set pFeatClassName = New FeatureClassName
Dim pDatasetName As IDatasetName
Set pDatasetName = pFeatClassName
'输出Dataset的Name为Union_result
pDatasetName.Name = m_pActiveLayer.Name & "_his"
Set pDatasetName.WorkspaceName = pNewWSName
' Set the tolerance. Passing 0.0 causes the default tolerance to be used.
' The default tolerance is 1/10,000 of the extent of the data frame's spatial domain
Dim pTol As Double
pTol = 0#
' Perform the union
Dim pBGP As IBasicGeoprocessor
Set pBGP = New BasicGeoprocessor
Dim pOutputFeatClass As IFeatureClass
pBGP.SpatialReference = m_pMap.SpatialReference
Set pOutputFeatClass = pBGP.Clip(pInputTable, False, pOverlayTable, False, pTol, pFeatClassName)
'/////////////////2.做差异分析,保留差异,全落在内的部分删除.////////////////
'找出被面覆盖的要素用于更新,删除或改变
Dim pOutCur As IFeatureCursor
Set pOutCur = pOutputFeatClass.Search(Nothing, False) '相交的要素
Dim pfeatureOut As IFeature
Dim pfeatureIn As IFeature
Set pfeatureOut = pOutCur.NextFeature
''''''''''''''''''''''''
Dim pFInClass As IFeatureClass
Dim pInCursor As IFeatureCursor
Dim pFilter As ISpatialFilter
' Dim pFeatureLayer As IGeoFeatureLayer
Set pFeatureLayer = m_pActiveLayer
Set pFInClass = pFeatureLayer.FeatureClass
Dim pTopo As ITopologicalOperator
Dim pGeom As IGeometry
Dim mArea As IArea
Do While Not pfeatureOut Is Nothing
Set pFilter = New SpatialFilter
With pFilter
Set .Geometry = pfeatureOut.Shape
.GeometryField = "SHAPE"
.SpatialRel = esriSpatialRelIntersects
Set .OutputSpatialReference(.GeometryField) = m_pMap.SpatialReference
End With
''''''''''''''''''''''''''''''''''''''''''''
Set pInCursor = pFInClass.Update(pFilter, False) '叠加分析出的图元,准备更新的
Set pfeatureIn = pInCursor.NextFeature
Do While Not pfeatureIn Is Nothing
'如果完全在多边形内
Set pTopo = pfeatureIn.Shape
Set pGeom = pfeatureIn.Shape
Set pGeom = pTopo.Difference(pfeatureOut.Shape) '叠加的差异部分
If Not pGeom Is Nothing And pGeom.IsEmpty = False Then
' Set mArea = pGeom
' If (mArea.Area > 100) Then
Set pfeatureIn.Shape = pGeom
pfeatureIn.Store
Else
' pfeatureIn.Delete
pInCursor.DeleteFeature
End If
Set pfeatureIn = pInCursor.NextFeature
Loop
'''''''''''''''''''''''''''''''''''''''''''''
Set pfeatureOut = pOutCur.NextFeature
Loop
|
