install.packages("maptools")#画地图的包
install.packages("spdep")#空间统计，moran'I
install.packages("tripack")
install.packages("RANN")
library("maptools")
library("spdep")
library("tripack")
library("RANN")

#读空间数据
rdata <- readShapePoly("D:/data/cairo.shp")#读取shp面数据
names(rdata)#显示数据字段名字
head(rdata@data)#显示属性表

#按照字段两幅画图
spplot(rdata[c("TFR96_03","TFR86_03") ],
       main = "spatial distribute of TFR96_03",#图名
       xlab = "X Coords",#横坐标名字
       ylab = "Y Coords",#纵坐标名字
       cut = 30#分段
       )

#按照字段一幅画图
spplot(rdata["TFR96_03"],
       main = "spatial distribute of TFR96_03",#图名
       xlab = "X Coords",#横坐标名字
       ylab = "Y Coords",#纵坐标名字
       cut = 30#分段
       )

#按边邻接角邻接生成邻居（方式一）
queen_nb <- poly2nb(rdata, queen = TRUE)#有8个
rook_nb <- poly2nb(rdata, queen = FALSE)#有4个

#获取中心点坐标编号
coords <- coordinates(rdata)

#地图数据转数据框-->生成ID
IDs <- row.names(as.data.frame(rdata))

#设置画参数
oopar <- par(mfrow = c(1,2),
             mar = c(3,3,1,1)+0.1)

#画边界
plot(rdata, border = "grey", main = "Queen-Style")

#add修改上一个图
plot(queen_nb, coords, col = "dodgerblue",
     add = TRUE, pch = 19, cex = 0.5)

#生成邻接关系（方式二）knearneigh()定义K-near
k4_nb <- knn2nb(knearneigh(coords, k = 4), row.names = IDs) 


#识别邻接关系是否对称
is.symmetric.nb(k4_nb, verbose = FALSE, force = TRUE)
#补全邻接关系 
n.comp.nb(k4_nb)$nc k4_w <- nb2listw(k4_nb)

#计算moran'I 
moran.test(rdata$TFR96_03, listw = k4_w) 
#计算moran'I (蒙特卡洛方法)
moran.mc(rdata$TFR96_03, listw = k4_w, nsim = 999)