arcgis中邻接矩阵文件的生成

首先我要说明的是，arcgis 9.2中已经实现基于基本一阶邻接关系的空间自相关计算，也
就是说arcgis系统本身已通过一定算法解决了本文以及之前blog中涉及的相关问题。

以下内容仅限于算法设计的兴趣及对arcgis相关功能实现的探讨。

arcgis 9.2中对于空间关联矩阵文件的构建做了拓展，你可以DBF属性中任意值为Unique的
字段来建立矩阵。原先矩阵的关联字段只能使用Object ID构建，而现在你可以使用其它唯
一值的字段了，甚至是中文的地名字段（当然需要注意的是，字段值中不能有空格！而且
似乎用户自定义的文件中若使用OID构建运行时会报错）

下面是之前的面向9.0版本开发的VBA代码，现增添了计算MoranI值的程序，一并打包（生
成的邻接文件可用于9.0/9.1版本的相关计算，也可以用于VBA自带的MoranI计算，但不能
直接用于9.2中）：
http://lilybbs.net/file/T/toolbar/VBA.rar

为了满足9.2中的邻接矩阵构建，鄙人近日对算法进行了扩充，并基于GDAL/OGR运算包在P
ython平台下进行了算法实现。特别要指出的是，由于GDAL相关函数的限制，算法中相应的
相邻多边形计算未能完全实现，仅仅搜索出了外接矩形（Envelop）相邻的多边形。因此大
家会发现生成的邻接矩阵行数会比原来的多一些。（相信随着GDAL的完善，最终可以在Py
thon平台上脱离ArcGIS实现相应功能）

附python源文件

  1 """
  2 多边形的邻接矩阵生成
  3 pnm.py - by toolbar
  4 
  5 update time: 2006.12.19
  6 """
  7 #引入操作系统库
  8 import os
  9 #引入TK库
 10 from Tkinter import *
 11 #引入消息对话框
 12 from tkMessageBox import *
 13 #引入文件对话框
 14 from tkFileDialog import *
 15 #引入_GDAL.dll
 16 import _gdal
 17 #引入OGR库
 18 from gdal import ogr
 19 
 20 #主程序类
 21 class Main:
 22     #字体定义
 23     _font=('宋体',9)
 24     #初始化
 25     def __init__(self,master):
 26         self.master=master
 27         Frame(master).pack(side=BOTTOM,expand=YES,fill=BOTH)
 28         master.title('多边形邻接计算'.decode('mbcs'))
 29         master.geometry('300x200+350+250')
 30         #第一行框架
 31         frame1=Frame(master)
 32         frame1.pack(side=TOP,expand=YES)
 33         #第二行框架
 34         frame2=Frame(master)
 35         frame2.pack(side=TOP,expand=YES)
 36         #第一行控件（文件输入）
 37         Label(frame1,text='输入文件'.decode('mbcs'),font=self._font).pack(side=LEFT,fill=BOTH)
 38         fin=StringVar()
 39         Entry(frame1,relief=SUNKEN,textvariable=fin,width=36).pack(side=LEFT)
 40         list=Listbox(frame2,height=3,width=37)
 41         Button(frame1,text='...',width=2,command=lambda w=fin,l=list,s=self:w.set(s.filein(l))).pack(side=LEFT)
 42         #第二行控件（文件输入）
 43         Label(frame2,text='选择变量'.decode('mbcs'),font=self._font).pack(side=LEFT,fill=BOTH)
 44         fieldin=StringVar()
 45         #list=Listbox(frame,height=3,width=37)
 46         scroll=Scrollbar(frame2,command=list.yview)
 47         list.config(yscrollcommand=scroll.set)
 48         list.pack(side=LEFT)
 49         scroll.pack(side=LEFT,fill=Y)
 50         #第三行框架
 51         frame3=Frame(master)
 52         frame3.pack(side=TOP,expand=YES)
 53         #第三行控件（文件输出）
 54         Label(frame3,text='输出文件'.decode('mbcs'),font=self._font).pack(side=LEFT,fill=BOTH)
 55         fout=StringVar()
 56         Entry(frame3,relief=SUNKEN,textvariable=fout,width=36).pack(side=LEFT)
 57         Button(frame3,text='...',width=2,command=lambda w=fout,s=self:w.set(s.fileout())).pack(side=LEFT)
 58         #第四行框架
 59         frame4=Frame(master)
 60         frame4.pack(side=TOP,expand=YES)
 61         #第四行控件（执行按钮）
 62         Button(frame4,text='确 认'.decode('mbcs'),font=self._font,
 63                command=lambda w1=fin,w2=fout,l=list,s=self:s.adjcent(w1.get(),w2.get(),l.selection_get())).pack(side=LEFT)
 64         Button(frame4,text='清 空'.decode('mbcs'),font=self._font,
 65                command=lambda w1=fin,w2=fout:(w1.set(''),w2.set(''))).pack(side=RIGHT)
 66         #GUI显示循环
 67         master.mainloop()        
 68 
 69     #文件输入
 70     def filein(self,list):
 71         a=askopenfilename(title='打开文件'.decode('mbcs'),filetypes=[('Shape file','*.shp'),('All files','*.*')])
 72         if a:
 73             shapefile=ogr.Open(a)
 74             shapelayer=shapefile.GetLayer()
 75             layerinfo=shapelayer.GetLayerDefn()
 76             for i in range(0,layerinfo.GetFieldCount()-1):
 77                 list.insert(END,layerinfo.GetFieldDefn(i).GetName())
 78             return os.path.abspath(a)
 79         else:
 80             return ''
 81 
 82     #文件输出
 83     def fileout(self):
 84         a=asksaveasfilename(title='保存文件'.decode('mbcs'),filetypes=[('Neighbor file','*.txt'),('All files','*.*')])
 85         if a:
 86             return os.path.abspath(a)
 87         else:
 88             return ''
 89 
 90     #执行邻接函数
 91     def adjcent(self,FileIn,FileOut,list):        
 92         textfile=open(FileOut,'w')
 93         textfile.writelines(list+'\n')
 94         shapefile=ogr.Open(FileIn)
 95         shapelayer=shapefile.GetLayer()
 96         #showinfo(title='',message=str(fcount))
 97         shapefeature=shapelayer.GetNextFeature()
 98         while shapefeature:
 99             shapegeometry=shapefeature.GetGeometryRef()
100             for i in range(0,shapelayer.GetFeatureCount()-1):
101                 shapefeature1=shapelayer.GetFeature(i)
102                 if shapefeature1.GetFID()==shapefeature.GetFID():continue
103                 shapegeometry1=shapefeature1.GetGeometryRef()
104                 if shapegeometry1.Intersect(shapegeometry):
105                     #s=str(shapefeature.GetFID())+' '+str(shapefeature1.GetFID())+' 1\n'
106                     t=shapefeature.GetFieldIndex(list)
107                     s=str(shapefeature.GetField(t))+' '+str(shapefeature1.GetField(t))+' 1\n'
108                     textfile.writelines(s)
109             self.master.title(str(shapefeature.GetFID()))
110             shapefeature=shapelayer.GetNextFeature()
111         textfile.close()
112         showinfo(message='Done!')
113 
114 if __name__=='__main__':
115     Main(Tk())

 


（python版本2.4.1，ArcGIS 9.2已自带；GDAL版本1.3.2，需安装）

谈论完了自己的算法，也很好奇ArcGIS 9.2中是怎么实现一阶邻接对象矩阵生成的。查找
了一下toolbox下的ArcScript源代码，发现ESRI的实现主要是通过了两个步骤：（1）使用
polygon to line功能，生成的line文件将自带有弧段左右多边形的拓扑信息；（2）对生
成line文件的左右多边形字段进行Frequency统计，剔除重复的冗余，在此基础上再生成邻
接矩阵文件就不难了（一次表格遍历即可）
有兴趣的话，可以在toolbox相应功能上右键‘edit’即可查看ESRI的源代码~
enjoy~

posted @ 2007-08-02 15:14  columbus2  阅读(4472)  评论(5)  编辑  收藏  举报