《ArcGIS点数据统计分析工具开发.docx》由会员分享，可在线阅读，更多相关《ArcGIS点数据统计分析工具开发.docx（17页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、ArcGIS点数据统计分析工具开发摘 要：空间点模式的分析方法主要包括样方分析、核密度估计和最近邻距离分析。本文主 要介绍了基于Python脚本语言的ArcGIS点模式分析工具箱的具体实现过程。关键词：空间点模式；Python； ArcGIS1引言任何的事件对象都可以抽象为空间中的一个位置点，地理事件点的分布模式包括随机分 布、均匀分布和聚集分布，针对点状事件的空间点模式分析应用广泛。文献1基于空间点模式分析了地震的空间分布趋势；文献2使用空间点模式分析确定 了大兴安岭地区雷击火的空间分格局；文献3使用空间点模式分析确定了保山市隆阳区地 质灾害的空间集聚特征；文献2基于空间点模式分析方法研究了
2、就业中心对不同产业的集 聚程度。在探究空间点事件的分布模式时，我们可以基于样方计数法、核密度估计法进行一阶效 应分析，同时可以基于最近邻距离进行二阶效应分析。本文主要介绍了样方分析、核密度估 计和最近邻距离分析的原理，基于Python语言开发相应工具的算法流程，并给出了相应的 主要代码，最后介绍了脚本工具界面中相应参数的设置。2点模式分析空间数据分析与常规的统计数据分析不同，在经典统计学中，研究样本是相互独立的, 而空间对象在空间位置上都具有一定的相关性，他们在空间上的变化不是随机的，且空间数 据的研究结果可能会因研究尺度和研究单元的组合方式而发生改变。探索性数据分析与常规 的数据分析模式也有
3、所不同，它防止了传统数据分析对数据模式、趋势和关系的假设。探索 性数据分析在没有先验假设的前提下，通过对数据的深入分析，最大限度的洞察数据，揭示 潜在结构，提取重要变量，探测奇异值，从而构成有检验价值的假设。从20世纪90年代开 始，广大地学研究者逐步把探索性数据分析方法应用到地理数据的分析中，来完善和开展空 间分析技术的理论和方法，从而形成了探索性数据分析(Exploratory Data Analysis)这一领 域。探索性空间统计分析可以用来概括空间数据的性质，探索空间数据的模式，产生和地 理数据相关的假设，并在地图上识别异常数据的分布位置，已发现是否存在热点区域等。空间点模式分析是研究
4、点状的地理实体或者事件的方法。点模式分析是探索性空间分析 的重要内容，因为任何的事件对象都可以抽象为空间中的一个位置点。地理事件点的分布模 式包括随机分布、均匀分布和聚集分布。对于空间点模式的一阶效应分析和二阶效应分析, 分别主要基于点密度和点距离的方法。当空间上过程的强度变化时，一阶效应显著，事件间 的绝对位置有决定作用，单位面积的事件数量在空间上有比拟清楚的变化，强调用观测事件 的空间密度估计，主要的分析方法有样方分析和核密度估计。当二阶效应强烈时，位置之间 存在相互作用，事件间的相对位置和距离具有决定作用，这时增加和减少相邻事件之间距离 的效果明显，强调用事件之间的距离估计，主要的分析方
5、法有K函数和基于最近邻距离的G4. 3最近邻距离工具最近邻距离分析原理最近邻距离指的是该模式中某一事件与其最近的事件之间的距离。位于Pi和Pj的两个 事件之间的表示为距离d（pg）,其计算使用欧基里德距离：d（Pi，Pj）; J（Xi Xj）2 + （ yj）2公式中，点pi和Pj的位置坐标分别为（Xi,yi）, （Xj,yj）o4. 3. 2最近邻距离分析流程最近邻距离的处理思路为使用SearchCursor方法遍历各要素点，获取其坐标位置，从而 遍历出与其最近的要素点及对应的距离。图3是本文基于点数据进行最近邻距离分析的算法 流程图。图7最近邻距离分析流程图4. 3. 3最近邻距离分析关键
6、步骤源代码导入模块，并设置相关参数：import arcpy,os,sysfc=arcpy.GetParameter(O) 使用For循环遍历图层中各点，获取各点x坐标和y坐标的数组:x = #x坐标数组y = #y坐标数组desc=arcpy.Describe(fc)shapefieldname = desc.ShapeFieldNamerows=arcpy.SearchCursor(fc)rows.reset()for row in rows:feat=row.getValue(shapefieldname)pnt=feat.getPart(O)x.append(pnt.X)y.appen
7、d(pnt.Y)使用While循环遍历各要素点，搜索与其距离最近的要素点及其对应距离# 定义计算两点之间距离的函数def distancedl,i2,j 1J2):d=math.sqrt(i 1 -i2)*(i 1 -+(j 11 -j2)return d# 使用While循环遍历各要素点，对各要素点的最近邻点及其距离的初始值进行赋值 n=len(x)dmin = 11 near = i = 0 while i n:if i = 0:dmin.append(distance(xi,xf 1 ,yi,y 1) near.append(l) i=i+l continuedmin.append(di
8、stance(xi,xO,yi,y 0) near.append(O) i=i+l# 使用While循环遍历各要素点，搜索各要素点的最近邻点及其距离 d 二口i = 0j = 1while i n: while j n: * d.append(O) ifli二j:d.append(distance(xi,xj,yi,yj)if = 2: del dj-2i=i+l(4)在点要素属性表中添加最近邻点字段和最近邻距离字段，使用For循环遍历图层中 各点，为最近邻字段和最短距离字段赋值：#属性表中添加最近邻点、最短距离字段field 1 _name=nnear_pointnfield2_name=H
9、near_disnarcpy.AddField_management(fc,fieldl_name,nSHORTn)arcpy.AddField_management(fc,field2_name,DOUBLE)#使用For循环遍历图层中各点，为最近邻字段和最短距离字段赋相应的值cursor=arcpy.da.UpdateCursor(fc, field l_name,field2_name)1 = 0for row in cursor:row0 = nearirowl = dminicursor. updateRow(row) i=i+l4. 3. 4自定义最近邻距离分析工具箱在点模式分析工
10、具中添加样方分析的脚本工具(Nearest Distance)。在Script File设置 关联的最近邻距离分析的脚本文件；然后设置样方分析的参数，包括输入型要素图层数据 Point layer,参数设置界面如图8所示。最近邻距离分析脚本工具相应参数设置完毕后，运 行脚本工具，界面如图9所示。Add Script图8最近邻距离分析参数设置回区三 Naerest Distance图9最近邻距离分析工具界面4应用实例在ArcMap的ArcToolbox中添加点模式分析工具箱(Point Analysis.tbx),在ArcMap 中即可直接调用相应工具进行分析。4. 3样方分析工具应用使用大新庄
11、东碑村地块数据进行样方分析。设置样方分析范围及单元格个数，工具界面 如图10所示，运行工具，结果统计了每个单元格所包含的事件点的数量，结果如图11所示。Quadrat AnalysisPoint layer |大新庄乡东碑村地块Boundary layer 庆新庄乡东碑村地块Extentlayer extentGri d (n*n)rnCancelEnvironments. Hide HelpGrid (nF)No description availablenrTool Help图10样方分析工具运行界面Quadrat AnalysisCompletedClose this dialog wh
12、en completed successfullyExecuting: QuadratAnalysis 大新庄乡东碑村 地块 大新庄乡东碑村地块layer extent 5 一 Start Time: Wed Apr 27 21:23:14 2016 Running script Quadra七Analysis 大新庄乡东碑村地块Shapeextent X: 1grid_0_0: 1, grid_0_l: 0, grid_0_ 21: 0r,grid_0_3,: 2,-1grid_0_41: 0, 1grid_3_4,:grid_3_31: 0f-grid_3_21: 17 grid_3_l1
13、: 27 7grid_3_01: 07 三 1grid_2_4,:grid_2_21: 0r-grid_2_31: 07 grid_2_0: 07 grid_2_l1: 07 1grid_4_4,: 0r-fgrid_4_01: 0r-grid_4_ I1: 17 7grid_4_2: 17 7grid_4_3: 0二一 1grid_l_l: 0,1grid_l_01: 0r-grid_l_ 31: 07 grid_l_2: 07 grid_l_41: 0? Completed script Quadra七Analysis Succeeded at Wed Apr 27 21:23:14 20
14、16 (Elapsed Time: 0.00 seconds)图H样方分析工具运行结果4. 3核密度估计工具应用使用某地气温数据进行核密度估计。设置栅格单元格大小及核密度分析的核带宽，工具 运行界面如图12所示，运行工具，分析结果如图13所示。口回 汉Search radiusNo description availableSearch radius200060000OK JCancelEnvironments. VBScript、JavaScript Jscritp 和Perl在内的脚本语言。与系统语言C+和.NET等相比，脚本语言更加简洁。VBA曾是ArcGIS 应用最普遍的编程工具，但A
15、rcGIS 10后，VBA不再作为ArcGIS的默认安装程序，Python 作为VBA的替代品被广泛使用。Python作为ArcGIS的脚本语言，具有以下优点：与C+、Visual Basic等高度结构化 的语言相比，语法简单，容易掌握；Python是一款开源软件；Python是一种解释性语言， 不需要编译可以直接运行，这使Python具有较强的可移植性；可伸缩性强，大型工程或小 型的一次性程序都可使用。ArcGIS 10引入了 ArcPy站点包，用来访问ArcGIS的地理处理功能。ArcPy是在 ArcGISScripting的基础上开展来的，继承了 ArcGISScripting的功能。通
16、过Python调用ArcPy 站点包可以进行地理数据分析、数据转换等。在使用ArcPy之前，首先要先导入ArcPy站 点包。导入ArcPy站点包后，可以调用ArcGIS地理处理工具箱中的所有工具网。ArcGIS 9.x版本已经引入Python工具箱，可实现自定义工具箱的功能，又便于创立新 的脚本工具。在ArcCatalog中添加一个新的工具箱(New Toolbox),在Toolbox中添加脚 本工具(Add Script)；将实现相应功能Python文件与该脚本工具进行关联，通过Script File 选取相应py文件；分析工具与Python文件关联后，设置相应参数，包括参数的数据格式、 输
17、入输出类型等。工具界面设置完成后，就可在ArcToolbox中添加工具箱(Add Toolbox), 在文件夹中选取相应的工具箱文件即可。4点模式分析工具箱开发样方分析工具4.1.1 样方分析原理样方分析方法是在一组固定大小的单元或样方内，记录各研究单元内的事件数量。样方 分析有两种计数方法，一种是对完全覆盖该研究区域的不重叠样方进行枚举检查，另一种是 在研究区域内随机设置样方并计算其事件数量。此处，我们开发的是基于枚举检查的样方分 析工具。4.1.2 样方分析流程样方分析的处理思路为根据网格行列数对网格进行编号，根据网格范围、网格数判断数 据点所在的网格，从而统计每个内的事件数量。图1是本文
18、基于点数据进行样方分析的算法 流程图。图1样方分析流程图样方分析关键步骤源代码导入模块，并设置相关参数：import arcpy, os, sys#点要素图层#边界图层#选取范围#格网范围layer = arcpy.GetParameter(O)layer2 = arcpy.GetParameterAsText( 1)extentType = arcpy.GetParameterAsText(2) nGrid = arcpy.GetParameter(3)使用If选择语句设置网格的范围：if(extentT ype=nLayer exrentn): ext=layer.getExtent()e
19、lse:mxd = arcpy.mapping.MapDocument(,CURRENTn)df = arcpy.mapping.ListDataFrames(mxd)Oext = df.extent 根据网格范围和网格数目计算单元格的宽度和高度：dx 二(ext.XMax - ext.XMin)/nGrid#单元格宽度dy = (ext.YMax - ext.YMin)/nGrid#单元格高度使用SearchCursor方法获取各要素点，使用For循环遍历图层中各要素点，根据单 元格的宽度、高度，判断其所在的网格编号，从而统计各网格内的要素点数目：# 使用For循环遍历各网格，将各网格内点要
20、素的数目的初始值设置为零grids = for i in range(nGrid):for j in range(nGrid):gridsngrid_n+str(i)+n_n+str(j)=O# 使用For循环遍历各点，根据单元格的宽度、高度，判断其所在的网格的编号，以此 统计各网格内的要素点数目desc = arcpy.Describe(layer)shapefieldname = desc.ShapeFieldNamerows = arcpy.SearchCursor(layer)rows.reset()for row in rows:feat 二 row.getValue(shapefi
21、eldname)pnt = feat.getPart(O)xidx = int (pnt.X-ext.XMin)Zdx) #计算行的编号yidx = int (pnt.Y-ext.YMin)/dy)#计算列的编号# 根据要素点的行列编号，使用If选择语句，判断其所在的网格的编号if xidx=nGrid:continueif yidx=nGrid:continuegridsngrid_n4-str(xidx)+n_n+str(yidx) = gridsngrid_H+str(xidx)+n_，+str(yidx) +1 del rows,rowarcpy.AddMessage(grids)自定
22、义样方分析工具箱在ArcCatalog中添加一个点模式分析的工具箱(Point Analysis, tbx),在点模式分 析工具中添加样方分析的脚本工具(Quadrat Analysis) o在Script File设置关联的样方 分析脚本文件；然后设置样方分析的参数，包括输入型要素图层数据Point layer、输入型要 素图层数据Boundary layer、输入型字符串数据Extent、输入型长整型数据Grid (n*n),参数 设置界面如图2所示。样方分析脚本工具相应参数设置完毕后，运行脚本工具，界面如图3 所示。图2样方分析参数设置图3样方分析工具界面4. 2核密度估计工具核密度估计
23、原理核密度估计是在样方分析的基础上开展而来的，与样方分析相比，核密度估计的优点在 于，它考虑了研究区域中各点之间的关系。研究范围内的任一点都有一个密度，不只是在有 事件发生的位置上，核密度通过计数一个区域中的，或中心在待估计位置的核中的事件数量 进行估计。，核密度估计的最简单的方法是在估计的位置上使用一个圆，然后计数圆内的事 件点的数目，然后用事件点的数目除以圆的面积。于是在点P有了一个核密度强度估计。_#(SeC(Plr) 入p-其中，c(pj)是中心在感兴趣点P，半径为r的圆。在核密度估计中，核带宽的选择具 有重要影响。如果带宽选择过大，估计的密度就会处处相似，且接近于整个模式的平均密度。
24、 但是如果带宽太小，估计结果就会强调单个事件，在远离事件的位置上，核密度为零。4. 2. 2核密度估计流程核密度估计的处理流程为根据点数据范围生成网格，判断以各网格为圆心，以带宽为半 径的圆内的事件的数量，然后除以圆的半径，从而求取每个单元格的核密度，进而获得整个 外表的核密度。图2是本文基于点数据进行核密度估计的算法流程图。图4核密度估计流程图在属性表中添加核密度强 度字段，并更新相应值CreateFishnet 工具 生成矢量化网格网格数据遍 历 各 点 的 x,y 坐 标属性表更新后的网格数据单元格的核密度强度原始点数据中各点的x, y坐标的数组FeatureToRaster 工具将矢量
25、数据转化为栅格数据栅格格式的核密度 估计数据遍 历 各 点 的 x,y 坐 标遍历单元格中心点数 据，判断原始点数据 中落在搜索半径内的 点，计算单元格强度2. 3核密度估计关键步骤源代码导入模块，并设置相关参数: import arcpy, os, sys#输入点要素图层#输出栅格图层#栅格图层的单元格大小#核带宽#输入点要素图层#输出栅格图层#栅格图层的单元格大小#核带宽layer = arcpy.GetParameter(O) outraster = arcpy.GetParameter( 1) cellsize = arcpy.GetParameter(2) r = arcpy. Ge
26、tParameter(3)使用CreateFishnet工具生(2)依据原始点要素图层范围和设置的输入单元格大小参数, 成矢量化的网格数据和各单元的中心点数据：#设置网格生成范围mxd=arcpy.mapping.MapDocument(HCURRENTn)df=arcpy.mapping.ListDataFrames(mxd)Oext=df.extent# 设置CreateFishnet方法的相应参数outfeatureclass=nE:/pythoexc/wwq.gdb/fishnetnorigincoordinate=arcpy.Point(ext.XMin,ext.YMin)yaxis
27、coordinate=arcpy.Point(ext.XMin,ext.YMax)cellsizewidth = cellsizecellsizeheight 二 cellsizenumrows = 0numcolum = 0oppositecoomor=arcpy.Point(ext.XMax,ext.YMax) labels = true, templateExtent = *#geometryType = POLYGON1# 调用 CreateFishnet 方法arcpy.CreateFishnet_management(outfeatureclass,origincoordinate
28、,yaxiscoordinate,cellsizew idth,cellsizeheight,numrows,numcolum,oppositecoornor,labels,templateExtent,geometryType)使用For循环分别遍历原始点数据中的各要素点和各单元格的中心点，分别获取原始 点数据和单元格中心点数据中各点的x坐标和y坐标的数组：# 定义获取图层中各点x坐标的数组的函数def getx(layer):desc=arcpy.Describe(layer)shapefieldname = desc.ShapeFieldNamex二口rows=arcpy.SearchC
29、ursor(layer)rows.reset()for row in rows:feat=row.getValue(shapefieldname)pnt=feat.getPart(O) x.append(pnt.X) return x# 定义获取图层中各点y坐标的数组的函数def gety(layer):desc=arcpy.Describe(layer)shapefieldname = desc.ShapeFieldName y =rows=arcpy.SearchCursor(layer)rows.reset()for row in rows:feat=row.getValue(shape
30、fieldname)pnt=fcat.gctPart(O) y.append(pnt.Y) return y#获取各单元格中心点中各点的x坐标和y坐标的数组以及单元格中心的个数 layer 1 =6shabeLshpx_fish=getx(layer 1)y_fish=gety(layer 1)n_fish=len(x_fish)#获取原始点数据中各点各点的x坐标和y坐标的数组以及单元格中心的个数x = getx(layer) y = gety(layer) n = len(x)使用While循环遍历各单元格中心点，内嵌While循环遍历原始点数据中的各点，根 据设置的核带宽r,判断原始点数据
31、中落在以各个单元格中心点为圆心，核带宽r为半径的 圆内的点的数目，并依据核密度估计公式：入=空誓 计算各单元格的核密度强度：p nr2# 计算每单元格的核密度强度 i = 0while i n_fish:#遍历单元格中心点j 二。k = 0while j n:#遍历原始数据各点dis=distance(x_fishi,xj,y_fishi,yj) if dis 输出型栅格图层Output raster、输入型双精度数据Cellsize、输入型双精度数据Searchradius,参数设置界面如图5所示。最近邻距离分析脚本工具相应参数设置完毕后，运行脚 本工具，界面如图6所示。Add Script
32、Click any parameter above to see its properties below.Di splay NameData TypePoint layerFeature LayerOutput rasterRaster LayerCellsi zeDoubleSearch radiusDoubleParameter PropertiesPropertyValueTypeRequiredDirecti on MultiValue Def aul t EnvironmentInput NoFilterNoneTo add a new parameter, type the name into an empty row in the name column, click in the Data Type column to choose a data type, then edit the Parameter Properties.上一步d) Finish取消图5核密度估计参数设置口回 汉KernelDensity3： Kernel Density Point layer圄 Output raster回 Cellsize Search radiusrrrCancelEnvironments. Hide HelpTool Help图6核密度估计工具界面