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源文件

"""
多边形的邻接矩阵生成
pnm.py - by toolbar

update time: 2006.12.19
"""
#引入操作系统库
import os
#引入TK库
from Tkinter import *
#引入消息对话框
from tkMessageBox import *
#引入文件对话框
from tkFileDialog import *
#引入_GDAL.dll
import _gdal
#引入OGR库
from gdal import ogr

#主程序类
class Main:
    #字体定义
    _font=('宋体',9)
    #初始化
    def __init__(self,master):
        self.master=master
        Frame(master).pack(side=BOTTOM,expand=YES,fill=BOTH)
        master.title('多边形邻接计算'.decode('mbcs'))
        master.geometry('300x200+350+250')
        #第一行框架
        frame1=Frame(master)
        frame1.pack(side=TOP,expand=YES)
        #第二行框架
        frame2=Frame(master)
        frame2.pack(side=TOP,expand=YES)
        #第一行控件（文件输入）
        Label(frame1,text='输入文件'.decode('mbcs'),font=self._font).pack(side=LEFT,fill=BOTH)
        fin=StringVar()
        Entry(frame1,relief=SUNKEN,textvariable=fin,width=36).pack(side=LEFT)
        list=Listbox(frame2,height=3,width=37)
        Button(frame1,text='...',width=2,command=lambda w=fin,l=list,s=self:w.set(s.filein(l))).pack(side=LEFT)
        #第二行控件（文件输入）
        Label(frame2,text='选择变量'.decode('mbcs'),font=self._font).pack(side=LEFT,fill=BOTH)
        fieldin=StringVar()
        #list=Listbox(frame,height=3,width=37)
        scroll=Scrollbar(frame2,command=list.yview)
        list.config(yscrollcommand=scroll.set)
        list.pack(side=LEFT)
        scroll.pack(side=LEFT,fill=Y)
        #第三行框架
        frame3=Frame(master)
        frame3.pack(side=TOP,expand=YES)
        #第三行控件（文件输出）
        Label(frame3,text='输出文件'.decode('mbcs'),font=self._font).pack(side=LEFT,fill=BOTH)
        fout=StringVar()
        Entry(frame3,relief=SUNKEN,textvariable=fout,width=36).pack(side=LEFT)
        Button(frame3,text='...',width=2,command=lambda w=fout,s=self:w.set(s.fileout())).pack(side=LEFT)
        #第四行框架
        frame4=Frame(master)
        frame4.pack(side=TOP,expand=YES)
        #第四行控件（执行按钮）
        Button(frame4,text='确 认'.decode('mbcs'),font=self._font,
               command=lambda w1=fin,w2=fout,l=list,s=self:s.adjcent(w1.get(),w2.get(),l.selection_get())).pack(side=LEFT)
        Button(frame4,text='清 空'.decode('mbcs'),font=self._font,
               command=lambda w1=fin,w2=fout:(w1.set(''),w2.set(''))).pack(side=RIGHT)
        #GUI显示循环
        master.mainloop()        

    #文件输入
    def filein(self,list):
        a=askopenfilename(title='打开文件'.decode('mbcs'),filetypes=[('Shape file','*.shp'),('All files','*.*')])
        if a:
            shapefile=ogr.Open(a)
            shapelayer=shapefile.GetLayer()
            layerinfo=shapelayer.GetLayerDefn()
            for i in range(0,layerinfo.GetFieldCount()-1):
                list.insert(END,layerinfo.GetFieldDefn(i).GetName())
            return os.path.abspath(a)
        else:
            return ''

    #文件输出
    def fileout(self):
        a=asksaveasfilename(title='保存文件'.decode('mbcs'),filetypes=[('Neighbor file','*.txt'),('All files','*.*')])
        if a:
            return os.path.abspath(a)
        else:
            return ''

    #执行邻接函数
    def adjcent(self,FileIn,FileOut,list):        
        textfile=open(FileOut,'w')
        textfile.writelines(list+'\n')
        shapefile=ogr.Open(FileIn)
        shapelayer=shapefile.GetLayer()
        #showinfo(title='',message=str(fcount))
        shapefeature=shapelayer.GetNextFeature()
        while shapefeature:
            shapegeometry=shapefeature.GetGeometryRef()
            for i in range(0,shapelayer.GetFeatureCount()-1):
                shapefeature1=shapelayer.GetFeature(i)
                if shapefeature1.GetFID()==shapefeature.GetFID():continue
                shapegeometry1=shapefeature1.GetGeometryRef()
                if shapegeometry1.Intersect(shapegeometry):
                    #s=str(shapefeature.GetFID())+' '+str(shapefeature1.GetFID())+' 1\n'
                    t=shapefeature.GetFieldIndex(list)
                    s=str(shapefeature.GetField(t))+' '+str(shapefeature1.GetField(t))+' 1\n'
                    textfile.writelines(s)
            self.master.title(str(shapefeature.GetFID()))
            shapefeature=shapelayer.GetNextFeature()
        textfile.close()
        showinfo(message='Done!')

if __name__=='__main__':
    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~