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)

这个我看过，我是想按年份分割后，直接生成各个截面的morani.不过还是谢了