用R语言进行莫兰指数分析并生成高清莫兰指数散点图

geoda上进行莫兰指数分析生成的图可编辑性较低，同时生成的图分辨率太低，于是想到了R语言。

R语言进行分析需要spdep包，自行安装。
代码如下：

 #设置工作所在路径
setwd("工作路径")

#读取相关包
library(sf)
library(spdep)
library(grid)
library(gridExtra)
library(ggplot2)

#读取数据
data <- st_read("file.shp")
 #计算权重表
weigth <- poly2nb(data) 
 #转换为权重矩阵
weigth_mattrix <- nb2listw(weigth)
#莫兰指数计算结果
moran_result <- moran.test(data$bianliang, listw = weigth_mattrix) 
#标准化变量，Geoda生成的散点图进行了标准化处理，如果不需要标准化，则跳过这一步，包括后面的代码进行修改
sc <- as.vector(scale(data$bianliang))

#莫兰指数散点图仅需变量和变量空间滞后值
#莫兰散点图和数据
dat <- moran.plot(sc, listw_1) 
View(dat)  #其中x为横坐标数据(变量)，wx为纵坐标数据(变量空间滞后值)

#用ggplot2画图
#注意：geom_hline(yintercept=0, lty=2) + geom_vline(xintercept=0, lty=2) 中的参考线如果是非标准化数据用均值代替，利用mean()函数

p1 <- ggplot(data = dat,aes(x=x,y=wx)) +
  #点属性
  geom_point(shape=1,colour="red",size=2,stroke = 0.7) +
 #拟合曲线属性
  geom_smooth(formula=y ~ x, method="lm", colour='black', se=FALSE, linewidth=0.5) +
  #x=0,y=0虚线
  geom_hline(yintercept=0, lty=2) + 
  geom_vline(xintercept=0, lty=2) + 
  #x轴刻度及范围
  scale_x_continuous(
    limits=c(-max(abs(floor(min(dat$x))), abs(ceiling(max(dat$x)))),max(abs(floor(min(dat$x))), abs(ceiling(max(dat$x))))),
    breaks=seq(-max(abs(floor(min(dat$x))), abs(ceiling(max(dat$x)))),max(abs(floor(min(dat$x))), abs(ceiling(max(dat$x)))), 1)) +
  #y轴刻度及范围
  scale_y_continuous(
    limits=c(-max(abs(floor(min(dat$x))), abs(ceiling(max(dat$x)))),max(abs(floor(min(dat$x))), abs(ceiling(max(dat$x))))),
    breaks=seq(-max(abs(floor(min(dat$x))), abs(ceiling(max(dat$x)))),max(abs(floor(min(dat$x))), abs(ceiling(max(dat$x)))), 1)) +
  #图上添加莫兰指数，p值，z值
  annotate('text',x=0.2,y=max(abs(floor(min(dat$x))), abs(ceiling(max(dat$x))))-0.5, 
           label= paste0("Moran's I:",round(moran_result_1$estimate[1],3),"\n","P:",
                         ifelse(round(moran_result_1$p.value,3)==0,0.001,round(moran_result_1$p.value,3)),
                         "\n","Z:",round(moran_result_1$statistic,3)),
           size=4,color='black',family = "serif",hjust=0) +
  #画多图时用于添加每个图的标签
  #annotate('text',x=-max(abs(floor(min(dat_1$x))), abs(ceiling(max(dat_1$x)))),
  #         y=max(abs(floor(min(dat_1$x))), abs(ceiling(max(dat_1$x)))), 
  #         label= "a", 
  #         size=5,color='black',family = "serif",hjust=0) +
  theme_bw() + 
  theme_classic() +
  theme(panel.grid = element_blank(),
        axis.ticks.length = unit(-0.1,"cm"),
        axis.text = element_text(family = "serif"),
        title = element_blank())

#多图排版
#merged_plot <- grid.arrange(grobs=list(p1,p2,p3,p4),nrow=2,ncol=2, 
#                            left=textGrob(paste0("滞后值","\n"," lag value"),
#                                          gp=gpar(fontsize = 15,fontfamily=c("SimSun","serif")),rot = 90), 
#                            bottom=textGrob(paste0("变量","\n","variable"),
#                                            gp=gpar(fontsize = 15,fontfamily=c("SimSun","serif"))))

#多图排版保存把plot=p1更换为plot=merged_plot
#保存图片
ggsave(filename = "moran.tiff", plot=p1, width=4, height=4, dpi=300)