R语言学习笔记——空间自相关：全局Moran’I/局部Moran’I /Geary’c/Moran散点图

data:

ChinaRD2<-readRDS("gadm36_CHN_1_sp.rds")
# download the file: https://gadm.org/download_country_v3.html
names(ChinaRD2)[4]<-c("province_name")
popDataSP<-merge(ChinaRD2, data_sp, by=c("province_name"),all.x=TRUE)
popDataSP<-popDataSP[!is.na(popDataSP$zao_va),]
summary(popDataSP)
coordinates(popDataSP)->coord
ID<-popDataSP@data$province_name
library(maptools)
library(spdep)
sh1<-knn2nb(knearneigh(coord,k=2),row.names=ID)# method 1
sh1_nb<-nb2listw(sh1)# method 1
soco_nbq <- poly2nb(popDataSP)# method 2
soco_nbq_w <- nb2listw(soco_nbq)#method 2
wm <- nb2mat(soco_nbq, style='B')#method 2
(moran(popDataSP$zao_va,sh1_nb,length(sh1),Szero(sh1_nb)))
moran.test(popDataSP$zao_va,sh1_nb,randomisation=F)
geary.test(popDataSP$zao_va, sh1_nb, randomisation=TRUE, zero.policy=NULL,
           alternative="greater", spChk=NULL, adjust.n=TRUE)
oid<-order(ID)
resl<-localmoran(popDataSP$zao_va,sh1_nb,p.adjust.method="bonferroni")
printCoefmat(data.frame(resl[oid,],row.names =popDataSP$province_name),cheek.names=FALSE)
moran.plot (popDataSP$zao_va,sh1_nb,labels =(popDataSP$province_name),pch=19)# method 1
moran.plot (popDataSP$zao_va,soco_nbq_w,labels =(popDataSP$province_name),pch=19,xlab = "First Trimester",ylab = "Vitamin A") #method 2
moran.mc(popDataSP$zao_va, soco_nbq_w, nsim=99)           

可以画出左上角的图