白话空间统计之二十五：空间权重矩阵（三）解构空间权重矩阵

中国古代就有“天圆地方”一说，所谓的“方”就是所谓的矩阵……

在军事上面，最为讲究就是团体的力量，需要集中多人之力，也为了方便指挥，所以就需要针对不同的情况，来组织阵型。

而最古老，但是一直沿用到今天的阵型，也是所谓的“方阵”。（凡阵有十：有方阵，有圆阵，有疏阵，有数阵，有锥行之阵，有雁行之阵，有钩行之阵，有玄襄之阵，有火阵，有水阵。此皆有所利。方阵者，所以剸（tuan，音：团，截断的意思）也。——孙膑兵法 十阵）

那么我们要将数据组织为最方便易读易用的形式，最简单的也是组织成矩阵（比如编程模型里面的数组，或者R语言里面的数据框（data frame），都是以矩阵的方式来进行数据组织的。

那么从空间权重矩阵这个名字来看，信息的组织形式自然也是通过矩阵来进行存储。

最简单的空间权重矩阵，示意如下：

但是在这种矩阵里面，我们发现了两个问题：

1、首先它是一个无向图（特殊矩阵），所以理论上只需要存储上（下）三角就可以了，全矩阵的存储，浪费了超过50%的存储空间。

2、里面如果有大量的0的话，实际上是可以不需要存储的（稀疏矩阵）。

稀疏矩阵是一种用三元组表达矩阵存储的方式，用于非0值的数量与上三角总元素数量比小于0.05时候，用来进行存储该矩阵的一种数据结构。

结构描述为：

(i,j,aij)

其中，i和j分别是行列号，aij是该行列号所代表元素的值。

那么上面的全矩阵，存储为稀疏矩阵的话，结构如下：

(1,2,1):A和B 为1，A是第一行，B是第二列，值为1。
(1,4,1):A和D为1，A是第一行，D是第四列，值为1。
(2,3,1):B和C为1，B是第二行，C是第三列，值为1。
(2,5,1):B和E为1，B是第二行，E是第五列，值为1。
(3,6,1):C和F为1，C是第三行，F是第六列，值为1。
(4,5,1):D和E为1，D是第四行，E是第五列，值为1。
(4,7,1):D和G为1，D是第四行，G是第七列，值为1。
(5,8,1):E和H为1，E是第五行，H是第八列，值为1。
(6,8,1):F和H为1，F是第六行，H是第八列，值为1。

解构完成。

可以看见，只需要8行，就可以完整的讲整个8*8的矩阵存储完成，节约了很多的存储空间，而且检索的时候，速度会快很多。

但是问题都是有两面性的，我们的矩阵是一个有向图（非特殊矩阵），而且还有可能组成了网状结构——所有的要素都与其他要素都有临近关系（矩阵的稠密度很高），那么解构成三元组的话，就会非常庞大了：

像这样的全网状结构，也就是矩阵中，0元素不存在……这种情况如果用三元组来表达，就会出现很庞大的情况。

ArcGIS的空间权重矩阵文件，是一个以swm为后缀名的二进制的文件——这也是ArcGIS空间统计模块一直为人诟病的一个历史遗留问题——无论是R语言的空间权重矩阵还是geoda的空间权重矩阵，都是ASCII明码表示的，我们可以打开它来直接进行查看或者进行手动修改。而ArcGIS的空间权重矩阵采用了二进制文件，给人感觉的就是一种浓浓的C语言时代的struct风格。

ArcGIS官方没有给出swm的结构，但是给出了读写的Python代码，所以可以猜测的结构如下：

第一行是swm的描述信息，如下所示：

'VERSION@10.1;UNIQUEID@UID;SPATIALREFNAME@GCS_WGS_1984;
INPUTFC@E:\\PyTestData\\blog\\python\\SWMDemo\\data\\world_chinese.shp;
WTYPE@5;DISTANCEMETHOD@EUCLIDEAN;EXPONENT@#;THRESHOLD@#;NUMNEIGHS@0;INPUTTABLE@#;
TIMEFIELD@#;TIMETYPE@#;TIMEVALUE@#;INPUTNET@#;IMPEDANCEFIELD@#;BARRIERFC@#;
UTURNPOLICY@#;RESTRICTIONS@#;USEHIERARCHY@#;SEARCHTOLERANCE@#;ADDCONCEPT@#;FIXEDWEIGHTS@True\n'

内容如下：

"VERSION@" 版本（ArcGIS以10.1为分界点，10.1之后，这里的版本信息全部默认为10.1）
"UNIQUEID@" 唯一值字段
"SPATIALREFNAME@" 空间参考
"INPUTFC@" 输入要素类（构建这个SWM的要素类）
"WTYPE@" 空间关系概念（一共有12种关系，如下：
wTypeDispatch = {0: 'INVERSE_DISTANCE',
                 1: 'FIXED_DISTANCE',
                 2: 'K_NEAREST_NEIGHBORS',
                 3: 'DELAUNAY_TRIANGULATION',
                 4: 'CONTIGUITY_EDGES_ONLY',
                 5: 'CONTIGUITY_EDGES_CORNERS',
                 6: 'CONVERT_TABLE',
                 7: 'ZONE_OF_INDIFFERENCE',
                 8: 'GET_SPATIAL_WEIGHTS_FROM_FILE',
                 9: 'SPACE_TIME_WINDOW',
                 10: 'NETWORK',
                 -1: 'UNKNOWN'}

）,
"DISTANCEMETHOD@" 距离方法（一般都是欧式距离）,
"EXPONENT@" 是否采用了距离指数（有否强调距离的作用）,
"THRESHOLD@" 是否采用了阈值,
"NUMNEIGHS@" 临近点的数量控制,
"INPUTTABLE@" 是否采用输入表生成的空间权重矩阵,
"TIMEFIELD@" 时间字段（用来做时空权重矩阵的，以下雷同）
"TIMETYPE@" 时间类型
"TIMEVALUE@" 时间的值,
"INPUTNET@" 输入的网络数据集（是否通过网络数据集构造的空间权重矩阵）,  
"IMPEDANCEFIELD@" 距离阻抗字段（网络数据集用的，以下雷同）,
"BARRIERFC@" 输入障碍要素,
"UTURNPOLICY@" 转弯策略（是否可以左右转弯）,
"RESTRICTIONS@" 限行策略,
"USEHIERARCHY@" 是否采用分级策略,
"SEARCHTOLERANCE@" 搜索容差,
"ADDCONCEPT@" 是否采用自顶向下的算法
"FIXEDWEIGHTS@" 是否采用了固定搜索方法。

看起来好像很复杂（好吧，实际上也很复杂），但是我们基本上不用去管，因为这些信息都是你创建空间权重矩阵时候输入的信息。而且在实际使用的时候，也没多大用……反正我没用过。

之后就是各种空间权重的信息了，基本结构如下：

Master ID : 主要素的ID。

Near ID : 与之临近的要素的ID。（可以有多个）

Weight : 权重。（与Near ID的数目相同）

示例如下：

35
12,22,32,46
0.25,0.25,0.25,0.25

如果让我们通过二进制文件的方法来读取，是一件很痛苦的事情，幸好ArcGIS并没有把“恶事”做绝……ArcGIS提供了一个官方方法，可以把swm文件转换为表格：

Spatial Statistics Tools —— Utilities ——Convert Spatial Weights Matrix to Table

下面我们来看看这个工具怎么用：

比起ArcGIS其他那些一堆参数的工具，这个工具极之简单……就一个输入和一个输出

1、输入你要转换的SWM文件。

2、选择你要输出的表格。

设置完成，点击OK之后，输出结果表格结果如下：

得到这样一个表格之后，你就可以在表格上进行修改，或者进行可视化了。

连线可视化的方法如下：

要连线，就必须有起点和终点的坐标，这些数据可以从原始数据里面去获取：

原始数据如下：包括了X,Y和UID三个关键字段。其中XY是坐标，UID是唯一标识码

将原始数据导出一份来（因为我们要同时设置UID和NID两个，而ArcGIS里面，不支持将同一份数据连续两次进行Join，所以需要导出一份来备用。

导出一份数据备用：

1、选择导出

2、选择新的名字：

3、保存：

4、在表格上点击右键，选择Joins and Relates  —— Join…

5、选择要关联的数据和字段。首先是要获得UID的数据，所以先输入的是UID的关联。

6、再次选择join，此次要关联的是NID，以获得NID的坐标：使用我们导出的备用数据：

挂接完成之后，将多余的字段勾掉，方便查看：

将图表导出来一份：

结果如下：

红色框里面是UID的坐标，蓝色框里面是NID的坐标，然后我们就可以进行连线了：

连线的工具如下：

在ArcToolbox 里面的 Data management Tools —— Features —— XY To Line

打开工具，设置好各种参数：

点击OK，生成连线：

到此可视化完成。

当然，这个过程也可以做成modelbuilder或者写成Python工具，这个代码以后有空在放出。

待续未完。