博文
shp2mat：基于shape文件生成空间权重矩阵
|||
【导读】 |
程序shp2mat基于shape文件生成各种形式的空间权重矩阵，一步到位。 |
一、函数动机 |
无论是用GeoDa、ArcGIS或是R语言spdep包等，生成的空间权重矩阵都存在如下问题：
（1）有可能存在有单元无邻居情况（GeoDa、spdep包相关函数），或是不能直接基于经纬度计算弧度距离（ArcGIS）；
（2）不是n*n矩阵格式，无法方便地用之于matlab空间计量经济分析函数。
针对以上这样那样的问题，这里推出基于R语言spdep包相关函数编写的函数shp2mat，直接从shape文件得到所需的各种n*n形式空间权重矩阵，使用非常方便。具体说明如下：
（1）元素主要采取两种取值方式：一是取值为0或1（binary），二是取值为单元间的距离（distance）。
（2）权重生成方法为邻接关系、k个最近邻居和阈值距离。空间单元为面和点。对于邻接关系，如是面，则根据queen规则赋值；如为点，则基于泰森多边形判断邻接关系。
（3）对于缺失邻居的单元i，如queen规则下的“孤岛”、距离阈值过小等，直接用其距离最近的邻居j填补之。同时在单元j增加邻居i，以保持其对称关系。
摒弃一些不必要的细枝末节，如邻接关系判断规则（queen或rook差别不大）、邻居阶数（2阶以上很少考虑）等，你只需设定函数取值形式、最近邻居数k或阈值距离td，就能得到理想且适用的空间权重矩阵。
二、R语言程序shp2mat |
shp2mat()基于shape文件生成n*n形式的空间权重矩阵，其参数设定如下：
shp2mat(shp, type = "b", td = NULL, k = NULL)
参数说明：
shp为shape文件路径，后缀名非必需；shape文件空间坐标用经纬度表示。
type取"b"-binary(0-1值元素)或"d"-distance(距离值元素)；
td为阈值距离,单位为km; k为最近邻居个数；
td和k取值只能二者居其一。
为不致混淆，空间权重矩阵单元排列顺序依照shape文件属性表顺序从1开始依次编号。
现取中国大陆31个省（面）和31个省会城市（点）为例说明其用法。
代码 |
> setwd("…\\shp2mat") > library(maptools) #查看shape文件数据情况 > provinces <- readShapeSpatial("province31_shp\\province31") > capitals <- readShapeSpatial("capitals31_shp\\capitals") > plot(provinces); plot(capitals, pch = 19, cex = 0.6, add = T) > provinces@data > capitals@data
#1 元素取值为0-1 #1.1 省份面 #1.1.1 基于queen邻接关系设定邻居,赋值为1 > pmat1 <- shp2mat("province31_shp\\province31") #1.1.2 将距离最近的3个单元设为邻居,赋值为1 > pmat2 <- shp2mat("province31_shp\\province31", k=3) #1.1.3 将500公里以内的单元设为邻居,赋值为1 > pmat3 <- shp2mat("province31_shp\\province31", td=500) #1.2 省会点 #1.2.1 基于泰森多边形邻接关系设定邻居,赋值为1 > cmat1 <- shp2mat("capitals31_shp\\capitals") #1.2.2 将距离最近的3个单元设为邻居,赋值为1 > cmat2 <- shp2mat("capitals31_shp\\capitals", k=3) #1.2.3 将500公里以内的单元设为邻居,赋值为1 > cmat3 <- shp2mat("capitals31_shp\\capitals", td=500)
#2 元素取值为距离,这里基于经纬度计算弧度距离,以km为单位 #2.1 省份面 #2.1.1 给出所有单元间距离 > pmat1 <- shp2mat("province31_shp\\province31", type="d") #2.1.2 给出最近三个邻居的距离 > pmat2 <- shp2mat("province31_shp\\province31", type="d", k=3) #2.1.3 给出500公里范围内的邻居距离,如无邻居则用距离最近的单元替代 > pmat3 <- shp2mat("province31_shp\\province31", type="d", td=500) #2.2 省会点 #2.2.1 给出所有点对间距离 > cmat1 <- shp2mat("capitals31_shp\\capitals", type="d") #2.2.2 给出最近三个点的距离 > cmat2 <- shp2mat("capitals31_shp\\capitals", type="d", k=3) #2.2.2 给出500公里范围内的邻居距离,如无邻居则用距离最近的点替代 > cmat3 <- shp2mat("capitals31_shp\\capitals", type="d", td=500)
#给距离矩阵行列赋名并写出权重矩阵 > unit.names <- provinces@data$NAME > dimnames(pmat1) <- list(unit.names, unit.names) > write.csv(pmat1, "province31_distance.csv") > unit.names <- capitals@data$NAME > dimnames(cmat1) <- list(unit.names, unit.names) > write.csv(cmat1, "capital31_distance.csv") |
数据和代码下载请关注微信公众号：“思达区域经济研究方法”，SDAR-workshop
交流请加微信群：“思达区域经济研究方法”，需邀请加入，请注明“单位名称+姓名”
扫码或长按，关注该微信号 |
版权申明 本号所有图文资料，除特别说明外，其版权归浙江工业大学王庆喜领衔的“思达工作室”所有。转发、引用请注明原始出处。 |
网络链接 1、科学网博客：http://blog.sciencenet.cn/u/workshopofsdar 2、网易博客：http://wqx1976.blog.163.com/ 3、人大经济论坛账号：R语言区域经济 4、知乎账号：sdar |
公众号功能 1、致力于打造中国R语言空间数据分析第壹号； 2、探讨R语言、ArcGIS、matlab、GeoDa、Stata、SPSS等软件在空间数据分析中的应用； 3、讨论产业集聚、空间溢出、区域创新网络、城市群发展等热门研究主题。 |
https://blog.sciencenet.cn/blog-3376208-1118132.html
上一篇：swm2mat：将swm权重文件转成n-by-n矩阵
下一篇：一招搞定空间横截面数据回归模型
