ArcEngine-空间分析解析课件.ppt

1、1第11讲 空间分析 国土信息工程系复习一下专题地图的几种样式 重点掌握几个对象 国土信息工程系11空间分析 空间分析是在对地理空间中的目标进行形态结构定义与分类的基础上，对目标的空间关系和空间行为进行捆，为目标的空间查询和空间相关分析提供参考，进一步为时间性服务的功能体系 通过空间分析不仅可以知道数据库的数据，而且可以通过这些数据云提示更深刻、更内在的规律和特征。衡量一个GIS系统的好坏，关键要看其空间分析功能是否强大，实用与灵活性 形态分析、空间查询、空间关系描述，空间相关分析和空间决策 国土信息工程系11.1空间查询与分析的基础 11.1.1基础概念 空间查询与分析的基础的基础是对空间图

2、形的运算操作 如：交、并、差、异或、缓冲，裁减，谁计算凸壳、切分，简化，移动，缩放，长度和面积 国土信息工程系11.1.2使用ITopologicalOperator接口来操作图形 ITopologicalOperator接口是AE当中非常重要的接口，实现它的对象包括：GeometryBag,Multipoint,Point,polyline,polygon.只要把运算的图形转换为五个当中的一个，就可以用它来运算 国土信息工程系组件对象模型图 国土信息工程系 Boundary:只读，获得图形的边框。面的国边框为polyline,线的边框为Multipoint,点或Multipoint的边框为空

3、或者Multipoint Buffer:缓冲区操作 Clip:裁剪操作 ConstructUnion:两个以上多边形的合并 Difference:差运算 Intersect:交运算 国土信息工程系 IsKnownSimple:Boolean,指示图形是否为简化的，拓朴是否正确 IsSimple:Boolean,指示是否使用了Simplify方法对图形进行了简化操作 Simplify:简化对象 SymmetricDifference:Xor，异或操作 Union:并操作 国土信息工程系Simplify方法 国土信息工程系simplify 国土信息工程系11.1.2.1Union 国土信息工程系U

4、nion代码(参看VB与VB2005的代码)并操作Public Function Union(g1 As IGeometry,g2 As IGeometry)As IGeometry Dim pTopo As ITopologicalOperator Dim pTopo2 As ITopologicalOperator If Not TypeOf g1 Is ITopologicalOperator Or Not TypeOf g2 Is ITopologicalOperator Then MsgBox 图形不能是裁切矩形 Union=g1 Exit Function End If pTopo

5、=g1 使用ITopologicalOperator接口 pTopo.Simplify 简化图形pTopo2=g2 pTopo2.SimplifyUnion=pTopo.Union(g2)合并pTopo=NothingpTopo2=NothingEnd Function 国土信息工程系11.1.2.2InterSect（交）操作 国土信息工程系InterSect(参看VB与VB2005的代码)Dim pTopo As ITopologicalOperatorDim pTopo2 As ITopologicalOperatorDim iDimension As esriGeometryDimen

6、sionIf Not TypeOf g1 Is ITopologicalOperator Or Not TypeOf g2 Is ITopologicalOperator ThenMsgBox(图形不能是裁切矩形)InterSect=g1Exit FunctionEnd IfpTopo=g1pTopo.Simplify()pTopo2=g2pTopo2.Simplify()iDimension=IIf(g1.Dimension g2.Dimension,g1.Dimension,g2.Dimension)InterSect=pTopo.InterSect(g2,iDimension)pTopo

7、=NothingpTopo2=Nothing 国土信息工程系11.1.2.3Difference(差)操作 国土信息工程系Difference(参看VB与VB2005的代码)Dim pTopo As ITopologicalOperatorDim pTopo2 As ITopologicalOperato If Not TypeOf g1 Is ITopologicalOperator Or Not TypeOf g2 Is ITopologicalOperator Then MsgBox(图形不能是裁切矩形)Difference=g1 Exit Function End IfIf g1.Di

8、mension g2.Dimension ThenMsgBox(两个图形的维数(Dimension)必需相同)Exit FunctionEnd IfpTopo=g1pTopo.Simplify()pTopo2=g2pTopo2.Simplify()Difference=pTopo.Difference(g2)pTopo=Nothing 国土信息工程系11.1.2.4SymmetricDifference（异或）操作 国土信息工程系SymmetricDifference(参看VB与VB2005的代码)Dim pTopo As ITopologicalOperatorDim pTopo2 As I

9、TopologicalOperatorIf Not TypeOf g1 Is ITopologicalOperator Or Not TypeOf g2 Is ITopologicalOperator ThenMsgBox(图形不能是裁切矩形)SymmetricDifference=g1Exit FunctionEnd IfIf g1.Dimension g2.Dimension ThenMsgBox(两个图形的维数(Dimension)必需相同)Exit FunctionEnd IfpTopo=g1pTopo.Simplify()pTopo2=g2pTopo2.Simplify()Symme

10、tricDifference=pTopo.SymmetricDifference(g2)pTopo=NothingpTopo2=Nothing 国土信息工程系11.1.2.5Clip（剪切）操作 国土信息工程系裁剪(参看VB与VB2005的代码)Public Function Clip(g1 As IGeometry,g2 As IGeometry)As IGeometry Dim pTopo As ITopologicalOperator If Not TypeOf g1 Is ITopologicalOperator Then MsgBox 图形1不能是裁切矩形 Set Clip=g1 E

11、xit Function End If If Not TypeOf g2 Is IEnvelope ThenMsgBox 绘制的第二个图形必需是裁剪矩形，请刷屏然后重新绘制！Set Clip=g1 Exit Function End If Set pTopo=g1 pTopo.Simplify pTopo.Clip g2Clip=pTopopTopo=NothingEnd Function 国土信息工程系11.1.2.6Cut（切分）操作 国土信息工程系11.1.2.6 Cut(切分)操作代码Public Function Cut(g1 As IGeometry,g2 As IGeometry

12、)As IGeometry Dim pTopo As ITopologicalOperator Dim pTopo2 As ITopologicalOperator Dim pL As IGeometry Dim pR As IGeometry If Not TypeOf g1 Is ITopologicalOperator Then MsgBox 图形1不能是裁切矩形Cut=g1 Exit Function End If If Not TypeOf g2 Is IPolyline Then MsgBox 绘制的第二个图形必需是Polyline，请刷屏然后重新绘制！Cut=g1 Exit Fu

13、nction End If MsgBox 屏幕上显示的是切分出来的左多边形！pTopo=g1 pTopo.SimplifypTopo2=g2 pTopo2.Simplify On Error GoTo Err pTopo.Cut g2,pL,pRCut=pLpTopo=NothingpTopo2=Nothing Exit FunctionErr:MsgBox 切分线可能无法完整的切分图形！End Function 国土信息工程系11.1.2.7 Buffer(缓冲区)操作 国土信息工程系11.1.2.7 Buffer（缓冲区）操作代码Public Function Buffer(g1 As I

14、Geometry,g2 As IGeometry)As IGeometry Dim pTopo As ITopologicalOperator Dim pTopo2 As ITopologicalOperator Dim pG1 As IGeometry Dim pG2 As IGeometry If Not TypeOf g1 Is ITopologicalOperator Then MsgBox 图形1不能是裁切矩形Buffer=g1 Exit Function End If If g1.Dimension g2.Dimension Then MsgBox 本示例代码要求图形1和图形2的维

15、数相同,以便合并操作！Buffer=g1 Exit Function End IfpTopo=g1 pTopo.SimplifypG1=pTopo.Buffer(50）pTopo2=g2 pTopo2.SimplifypG2=pTopo2.Buffer(50）Buffer=Union(pG1,pG2)pTopo=NothingpTopo2=NothingEnd Function 国土信息工程系11.1.2.8 Boundary(提取边界)操作 国土信息工程系Boundary(提取边界)操作代码Public Function Boundary(g1 As IGeometry,g2 As IGeo

16、metry)As IGeometry Dim pTopo As ITopologicalOperator Dim pTopo2 As ITopologicalOperator Dim pG1 As IGeometry Dim pG2 As IGeometry If Not TypeOf g1 Is ITopologicalOperator Then MsgBox 图形1不能是裁切矩形Boundary=g1 Exit Function End If If g1.Dimension g2.Dimension Then MsgBox 本示例代码要求图形1和图形2的维数相同,以便合并操作！Bounda

17、ry=g1 Exit Function End IfpTopo=g1 pTopo.SimplifypG1=pTopo.Boundary pTopo2=g2 pTopo2.SimplifypG2=pTopo2.Boundary Boundary=Union(pG1,pG2)pTopo=NothingpTopo2=NothingEnd Function 国土信息工程系ITopologicalOperator注意事项 1）两个Geometry不仅要有坐标系统，而且必须是相同的坐标系统。如果不是，把一个Geometry投影到另一个的坐标系统中去。先设置IGeometry:SpatialReferenc

18、e，然后调用IGeometry:Project并把坐标系统作为参数。有时没有必要设置坐标系统，只要调用Project就行了 国土信息工程系 2）Snap both geometries to the spatial reference.This can be done by calling IGeometry:SnapToSpatialReference.捕捉Geometry到坐标系统。可以通过调用IGeometry:SnapToSpatialReference来实现。国土信息工程系 3）Simplify both geometries.This can be done by a QI to

19、ITopologicalOperator2 and setting IsKnownSimple to False and then calling Simplify.保证两个Geometry都是拓扑简单的。先通过QI查询ITopologicalOperator2接口并设置IsKnownSimple为False，然后调用Simplify。国土信息工程系11.2空间查询 空间查询是使用一定的条件在空间数据中搜索相应的结果集，和一般的查询不同的是，这种查询是对空间对象的搜索，而不是对文字信息的搜索。基于空间属性的查询 基本空间位置的查询 国土信息工程系11.2.1QuerFilter对象的重要属性方

20、法（参看程序）Addfield:向QuerFilter中添加一个字段到Subfields用来查询 Outputspatialreference:以字段设置或者读取空间参考 Where clause:SQL语句（不含select*from查询）FeatureCursor是一个可以包含多个Feature对象，他实现了I Feature接口与ICursor接口，通过：IFeatureCursor nextfeature方法可以遍历所有Feature对象，通过ICursor:nextRow，可以所有属性记录 国土信息工程系Geodatabase类库当中 国土信息工程系11.2.2基于空间位置的查询 I

21、spationFiter 国土信息工程系 Geometryfield:图形字段的字段名称,即shape字段 Spatialrel:查询的空间关系。相交，相接，穿越，包含等 Spataildescription:自定义的查询的空间关系 Geometry:用来查询的图形，比如拉框查询的“框”国土信息工程系SpatilRel属性的类型esrispatiolrelenum,定义如下 国土信息工程系ISpatialFilter.SpatialRelDescription 国土信息工程系点查询 Dim pPt As esriGeometry.Point Dim pGeo As IGeometry Dim

22、pSpatialFilter As ISpatialFilter Dim pFeaturelayer As IFeatureLayerDim pTopo As ITopologicalOperator 创建点查询的“点”,注意在点查询中，真正的点对于点、线对象的相交是很难的 因此我们设置一个容差，用一个很小的面来代替点是很好的做法 Set pPt=MapControl.ToMapPoint(x,y)Set pTopo=pPt Set pGeo=pTopo.Buffer(0.01)创建SpatialFilter Set pSpatialFilter=New SpatialFilter With

23、pSpatialFilter Set.Geometry=pGeo .SpatialRel=esriSpatialRelIntersects End With 国土信息工程系执行查询 pSpatialFilter.GeometryField=pFeaturelayer.FeatureClass.ShapeFieldName Set pFeatureCursor=pFeaturelayer.Search(pSpatialFilter,False)Set pFeature=pFeatureCursor.NextFeature 国土信息工程系多边形查询 Set pPolygon=MapControl.

24、TrackPolygon 创建SpatialFilter Set pSpatialFilter=New SpatialFilter With pSpatialFilter Set.Geometry=pPolygon .SpatialRel=esriSpatialRelIntersects End With 国土信息工程系线查询 创建点查询的“线”Set pPolyline=MapControl.TrackLine 创建SpatialFilter Set pSpatialFilter=New SpatialFilter With pSpatialFilter Set.Geometry=pPoly

25、line .SpatialRel=esriSpatialRelIntersects End With 国土信息工程系11.3联合查询 Set pPolygon=MapControl.TrackPolygon 创建SpatialFilter Set pSpatialFilter=New SpatialFilter With pSpatialFilter Set.Geometry=pPolygon .SpatialRel=esriSpatialRelIntersects .WhereClause=pLabel.Caption End With 国土信息工程系IRubberBand 经常需要在项目中

26、画地理要素（Feature）时或者画元素（Element），其实IRUbberBand接口就实现了绘制几何形体（Geometry）的方法TrackNew，以及移动一个一个几何形体的方法TrackExisting。IRUbberBand接口有两个方法：1，TrackExisting方法；2，TrackNew方法；国土信息工程系11.4地图代数功能1:求面积 Dim pActView As IActiveView=AxMapControl1.Map Dim pScreen As IScreenDisplay=pActView.ScreenDisplay Dim pRubber As IRubber

27、Band=New RubberPolygonClass()Dim pPolygon As IPolygon=TryCast(pRubber.TrackNew(pScreen,Nothing),IPolygon)Dim pArea As IArea=TryCast(pPolygon,IArea)Dim dArea As Double=pArea.Area MessageBox.Show(面积为：&dArea)国土信息工程系地图代数功能2:求多边形线的长度 Dim pActView As IActiveView=AxMapControl1.Map Dim pScreen As IScreenDis

28、play=pActView.ScreenDisplay Dim pRubber As IRubberBand=New RubberLine Dim pPol As IPolyline=TryCast(pRubber.TrackNew(pScreen,Nothing),IPolyline)Dim dlen As Double=pPol.Length MessageBox.Show(长度为：&dlen)国土信息工程系以下是VBA的一些功能代码 使用方法1打开属性表,选择计算的字段,右点选择Calculate Values；2.选择“是”，进入Field Calculator；2选择Advance选

29、项；3 在Pre-Logic VBA Script Code编辑框中输入VBA代码；4在下面编辑框中输入赋值部分.国土信息工程系1 多边形面积VBA部分:Dim pGeo As IGeometrySet pGeo=ShapeDim pPolygon As IPolygonSet pPolygon=pGeoDim pArea As IAreaSet pArea=pPolygon赋值部分:pArea.Area 国土信息工程系2-多边形周长 VBA部分:Dim pGeo As IGeometrySet pGeo=ShapeDim pPolygon As IPolygonSet pPolygon=pG

30、eo赋值部分:pPolygon.Length 国土信息工程系3-多边形重心X VBA部分:Dim pGeo As IGeometrySet pGeo=ShapeDim pPolygon As IPolygonSet pPolygon=pGeoDim pArea As IAreaSet pArea=pPolygonDim pPoint As IPointSet pPoint=pArea.Centroid赋值部分:pPoint.X 国土信息工程系5-多边形重心YVBA部分:同上赋值部分:pPoint.Y 国土信息工程系6-Polyline长度 VBA部分:Dim pGeo As IGeometry

31、Set pGeo=ShapeDim pPolyline As IPolylineSet pPolyline=pGeoDim pCurve As IPolycurveSet pCurve=pPolyline赋值部分:pCurve.Length 国土信息工程系7-点坐标XVBA部分:Dim pGeo As IGeometrySet pGeo=ShapeDim pPoint As IPointSet pPoint=pGeo赋值部分:pPoint.X 国土信息工程系8-点坐标Y VBA部分:同上赋值部分:pPoint.Y 国土信息工程系 9-表示点坐标XVBA部分:Dim pDoc As IMxDoc

32、umentSet pDoc=ThisDocumentDim pSpRef As ISpatialReferenceSet pSpRef=pDoc.FocusMap.SpatialReferenceDim pClone As ICloneSet pClone=ShapeDim pGeo As IGeometrySet pGeo=pClone.CloneDim pPoint as IPointSet pPoint=pGeopGeo.Project pSpRef赋值部分:pPoint.X 国土信息工程系10-表示点坐标YVBA部分:同上赋值部分:pPoint.Y坐标值为On the Fly显示的坐标

33、,不是文件存储的固有坐标 国土信息工程系11-连续编号VBA部分:Static lCount as longlCount=lCount+1赋值部分:lCount(从1开始)lCount-1(从0开始)国土信息工程系11.5缓冲区分析 地理目标的延展或扩展，在实际问题中可以理解为地理目标的影响范围，如污染源的影响范围。分析步骤：获得要被进行缓冲运算的图形（Igeometry）;第二步对获得的图形进行缓冲运算生成一个新的图形（Igeometry）；第三步：再利用这个新的图形进行复杂的空间查询等操作。国土信息工程系 Dim pPt As IPoint Dim pGeo As IGeometry Di

34、m pMapCtrl As MapControl Dim pFeatLyr As IFeatureLayer Dim pFSelection As IFeatureSelection Dim pFilter As ISpatialFilter Dim pTopo As ITopologicalOperator Dim i As Integer,Count As Integer Set pMapCtrl=m_pHook.hook Set pPt=pMapCtrl.ToMapPoint(x,y)Set pTopo=pPt Count=pMapCtrl.LayerCount 国土信息工程系 Set

35、pFilter=New SpatialFilter With pFilter Set.Geometry=pTopo.Buffer(0.05).SpatialRel=esriSpatialRelIntersects End With For i=0 To Count-1 Set pFSelection=pMapCtrl.Layer(i)pFSelection.SelectFeatures pFilter,esriSelectionResultNew,False Next pMapCtrl.Refresh esriViewGeoSelection 国土信息工程系11.6叠加分析 叠置分析（over

36、lay analysis）是地理信息系统提取空间隐含信息的手段之一。它是将一个区域两组或两组以上的地物要素进行叠加，从而产生新特征的分析方法。矢量叠置分析与叠置栅格分析 国土信息工程系11.6.1基于矢量的叠置分析 点与多边形的叠置分析：比如分析城市中每个学校的数量 线与多边形的叠置分析：如：科用一个河流要素层，可以找出哪条河流在哪个宗地内以及这条河流有多长一段落在这块宗地内等 多边形与多边形的叠置分析：如将宗地地块要素层与污染区域要素层进行叠置分析（如把宗地要素层和水系要素层进行叠置，来获得河流在宗地上的边界信息，擦除叠置）国土信息工程系 AE 开发中，矢量图层叠加分析需要用到的主要类为 B

37、asicGeoprocessor，其主要接口为 IBasicGeoprocessor。IBasicGeoprocessor接口提供了基本的空间数据处理的方法和属性，其中包括叠加求交（Interset）和叠加求和（Union)。国土信息工程系IBasicGeoprocessor Interface MembersDescription CancelTracker The cancel tracker.ClipClips features.Dissolve Dissolves features.Intersect Intersects features.MergeMerges features.S

38、patialReferenceThe output spatial reference.UnionCreates a union of features.BasicGeoprocessor：CoClasses that implement IBasicGeoprocessor 国土信息工程系下面提供一个叠加求交的开发实例：Private Sub M_OverLayer_Click(ByVal sender As Object,ByVal e As System.EventArgs)Try 分析层 Dim pLayer As ILayer=Me.axMapControl1.get_Layer(0

39、)Dim pInputFeatLayer As IFeatureLayer=TryCast(pLayer,IFeatureLayer)Dim pInputTable As ITable=TryCast(pLayer,ITable)Dim pInputFeatClass As IFeatureClass=pInputFeatLayer.FeatureClass 叠加表 pLayer=Me.axMapControl1.get_Layer(1)Dim pOverlayTable As ITable=TryCast(pLayer,ITable)国土信息工程系 叠加分析表 Dim pFeatClassN

40、ame As IFeatureClassName=New FeatureClassNameClass()pFeatClassName.FeatureType=esriFeatureType.esriFTSimple pFeatClassName.ShapeFieldName=shape pFeatClassName.ShapeType=pInputFeatClass.ShapeType 工作空间名称 Dim pNewWSName As IWorkspaceName=New WorkspaceNameClass()pNewWSName.WorkspaceFactoryProgID=esriDat

41、aSourcesFile.ShapefileWorkspaceFactory pNewWSName.PathName=C:temp 国土信息工程系 数据集名称 Dim pDatasetName As IDatasetName=TryCast(pFeatClassName,IDatasetName)pDatasetName.Name=ss pDatasetName.WorkspaceName=pNewWSName 几何处理 Dim pBGP As IBasicGeoprocessor=New BasicGeoprocessorClass()Dim pOutputFeatClass As IFea

42、tureClass=pBGP.Intersect(pInputTable,False,pOverlayTable,False,0.01,pFeatClassName)国土信息工程系 输出要素层设置 Dim pOutputFeatLayer As IFeatureLayer=New FeatureLayerClass()pOutputFeatLayer.FeatureClass=pOutputFeatClass pOutputFeatLayer.Name=pOutputFeatClass.AliasName Me.axMapControl1.AddLayer(DirectCast(pOutput

43、FeatClass,ILayer),0)axMapControl1.Update()Catch ex As Exception MessageBox.Show(ex.Message)End Try End Sub 国土信息工程系11.6.2基于栅格的数据分析 栅格数据具有空间信息隐含但属性信息明显的特点，利用多层栅格数据进行各种各样的叠加，可以发现大量的空间现或空间过程，比如将同一地区不同时期的影像叠加，可以获得城市扩张，森林砍伐，洪水淹没，山体滑坡等信息。国土信息工程系RasterMathOps组件组件 The RasterMathOpsmethods provide access to a

44、 full suite of mathematical operators and functions.These operators and functions enable the values in multiple rasters to be combined arithmeticallyaddition,subtraction,multiplication,and division,for example.In addition,the mathematical manipulation of the values in a single input raster(sine,expo

45、nent,power,and so on),the evaluation of multiple input rasters(Boolean And,Greater Than,combinations,and so forth),and the evaluation and manipulation of the values in the binary format(Bitwise And and Bitwise Left Shift,for example)are also supported by these methods.国土信息工程系esriSpatialAnalyst.raste

46、rmathops 参看SpatialAnalyst模型图 1）基于栅格的代数运算 2）基于数学变换的数学变换（指数，对数，三角变换）等 3）基于多上栅格要素层的运算，包括代数运算（+-/*）和逻辑运算等 国土信息工程系Imathop Abs：用于计算栅格像无的绝对值并写入新数据集。Divide:用于把两输入栅格像元值并写入新数据集 Exp:用于栅格像元以e为底的次方并写入新数据集 Exp10 Exp2 Float Int 国土信息工程系 Log10 Log2 Minus:两像元相减 Mod Negatge Plus Power Square Times 国土信息工程系 ILogialOp接口：

47、逻辑运算方法 ITrigOP接口：提供了栅格三角运算方法。国土信息工程系 Dim pRas01 As IRasterSet pRas01=OpenDefaultIRaster(C:Data,InRas01)Dim pRas02 As IRasterSet pRas02=OpenDefaultIRaster(C:Data,InRas02)Dim pEnv As IRasterAnalysisEnvironmentSet pEnv=pMapAlgebraOpSet pEnv.OutWorkspace=OpenIRasterWorkspace(C:Workspace)Dim pRasOut As I

48、Raster 国土信息工程系 Set pMapAlgebraOp=New RasterMapAlgebraOppMapAlgebraOp.BindRaster pRas01,Ras01pMapAlgebraOp.BindRaster pRas02,Ras02Set pRasOut=pMapAlgebraOp.Execute(Ras01+Ras02)Save Output RasterDim sRasNmOut As StringsRasNmOut=outRas01RasterSaveOnDisk(pRasOut,C:Workspace,sRasNmOut,GRID)在帮助中找到了以上的代码，能不能Set pRasOut=pMapAlgebraOp.Execute(Ras01+Ras02)这个代码多加几个数据集就可以了