R语言计算moran‘I

R语言计算moran‘I

R语言计算moran‘I

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)

 

05-11 18:19