基于R的生态系统服务簇解释

最新推荐文章于 2024-12-04 17:36:23 发布

HHXX好好学习

最新推荐文章于 2024-12-04 17:36:23 发布

阅读量1.8k

点赞数 9

文章标签： r语言

本文链接：https://blog.youkuaiyun.com/xzy123138/article/details/129906253

版权

自组织映射（Self-Organizing Maps, SOM）是一种神经网络算法，可以用于聚类分析，由芬兰学者Kohonen提出，在R语言中对应的工具包是kohonen。正文如下：

#MGWR_用生态系统服务簇解释------------
setwd("E:/GIS data/GEP_YREB/R")
install.packages("kohonen")
install.packages("kohonen")
library(kohonen)
library(readxl)
library(ggplot2)
#watershed scale input data
data=read_xls("MGWR_CU.xls",sheet="1")
#data pre-processing
str(data)
X<-scale(data[,-1]) #由于各生态系统服务量纲不同，需进行标准化，-1表示除去第一项
#数据处理
WY=X[,2:5]
SC=X[,c(3:5,7)]
CF=X[,c(1,3:5)]
# WY--------------
set.seed(210000000)
g<-somgrid(xdim=4,ydim=1,topo="rectangular")
#xdim和ydim确定最终聚类数量，例如此处确定最终聚类数为6。由于是非监督的聚类，建议通过多次实验确定聚类数量使得：1）各类别间差异明显；2）各类别具有可解释性
map<-som(WY,grid=g,radius=1)
#X为标准化后的输入数据，grid为网格设置，radius为邻接半径，各参数含义和设置参考“kohonen”包说明
plot(map, type="counts")
#左图展示每一类别样本数量
plot(map, type="codes")#SOM聚类结果。
#右图花瓣越长代表某一特征值越高。
#output results
map$codes #展示各类别特征
map$unit.classif #展示各所属生态系统服务簇
output1<-cbind(data,map$unit.classif)#连接开始
write.table(output1,file="WY_CU13.csv",sep = ",",row.names=FALSE)
# SC--------------
set.seed(210000000)
g2<-somgrid(xdim=4,ydim=1,topo="rectangular")
#xdim和ydim确定最终聚类数量，例如此处确定最终聚类数为6。由于是非监督的聚类，建议通过多次实验确定聚类数量使得：1）各类别间差异明显；2）各类别具有可解释性
map2<-som(SC,grid=g2,radius=1)
#X为标准化后的输入数据，grid为网格设置，radius为邻接半径，各参数含义和设置参考“kohonen”包说明
plot(map2, type="counts")
#左图展示每一类别样本数量
plot(map2, type="codes")#SOM聚类结果。
#右图花瓣越长代表某一特征值越高。
#output results
map2$codes #展示各类别特征
map2$unit.classif #展示各所属生态系统服务簇
output2<-cbind(data,map2$unit.classif)
write.table(output2,file="SC_CU13.csv",sep = ",",row.names=FALSE)
# CF--------------
set.seed(210000000)
g3<-somgrid(xdim=4,ydim=1,topo="rectangular")
#xdim和ydim确定最终聚类数量，例如此处确定最终聚类数为6。由于是非监督的聚类，建议通过多次实验确定聚类数量使得：1）各类别间差异明显；2）各类别具有可解释性
map3<-(CF,grid=g3,radius=1)
#X为标准化后的输入数据，grid为网格设置，radius为邻接半径，各参数含义和设置参考“kohonen”包说明
plot(map3, type="counts")
#左图展示每一类别样本数量
plot(map3, type="codes")#SOM聚类结果。颜色设置：NL（#fdf06f）、TDI(#c5beaa)、NDVI(#17c37b)、HAI(#e57066)
#右图花瓣越长代表某一特征值越高。
#output results
map3$codes #展示各类别特征
map3$unit.classif #展示各所属生态系统服务簇
output3<-cbind(data,map3$unit.classif)
write.table(output3,file="CF_CU13.csv",sep = ",",row.names=FALSE)
#输出所属生态系统服务簇

pdf (file = "yebMGWR_CU10.pdf", width = 13, height =4)
CU=cowplot::plot_grid( wycu, sccu, cfcu,ncol = 3)
#svg (file = "yebtd_DEM2.svg", width = 17.1, height =4)
CU
dev.off() #must close device to close file!

#2 绘图--------------
#2.1 柱状图----------------
cu=na.omit(read_xls("MGWR_CU.xls",sheet="66"))#数据清洗
cu2=read_xls("MGWR_CU.xls",sheet="66")
#data pre-processing
cols=c("#4d73b1","#4daa92","#fedd90","#de3026")
#2.1.1 WY---------
wycu=ggplot(cu2, aes(fill=cu2[[2]], y=cu2[[3]], x=cu2[[1]])) +
geom_bar(position="fill", stat="identity") + scale_fill_manual(values = cols)+ theme(axis.line = element_line(size = 0.5,
linetype = "solid"), axis.ticks = element_line(colour = "black"),
panel.grid.major = element_line(linetype = "blank"),
panel.grid.minor = element_line(colour = "gray72",
linetype = "dashed"), axis.text = element_text(family = "serif",
size = 16, colour = "gray0"), axis.text.x = element_text(family = "serif"),
axis.text.y = element_text(family = "serif"),
plot.title = element_text(family = "serif"),
panel.background = element_rect(fill = "gray100",
colour = "aquamarine4"), plot.background = element_rect(fill = "white",
colour = NA, linetype = "solid"),
legend.position = "none") +labs(x = NULL, y = NULL, fill = NULL)+coord_fixed(ratio =7 / 1)
#2.1.2 sc---------
sccu=ggplot(cu, aes(fill=cu[[5]], y=cu[[6]], x=cu[[4]])) +
geom_bar(position="fill", stat="identity") + scale_fill_manual(values = cols)+ theme(axis.line = element_line(size = 0.5,
linetype = "solid"), axis.ticks = element_line(colour = "black"),
panel.grid.major = element_line(linetype = "blank"),
panel.grid.minor = element_line(colour = "gray72",
linetype = "dashed"), axis.text = element_text(family = "serif",
size = 16, colour = "gray0"), axis.text.x = element_text(family = "serif"),
axis.text.y = element_text(family = "serif"),
plot.title = element_text(family = "serif"),
panel.background = element_rect(fill = "gray100",
colour = "aquamarine4"), plot.background = element_rect(fill = "white",
colour = NA, linetype = "solid"),
legend.position = "none") +labs(x = NULL, y = NULL, fill = NULL)+coord_fixed(ratio =4 / 1)

#2.1.3 cf---------
cfcu=ggplot(cu, aes(fill=cu[[8]], y=cu[[9]], x=cu[[7]])) +
geom_bar(position="fill", stat="identity") + scale_fill_manual(values = cols)+ theme(axis.line = element_line(size = 0.5,
linetype = "solid"), axis.ticks = element_line(colour = "black"),
panel.grid.major = element_line(linetype = "blank"),
panel.grid.minor = element_line(colour = "gray72",
linetype = "dashed"), axis.text = element_text(family = "serif",
size = 16, colour = "gray0"), axis.text.x = element_text(family = "serif"),
axis.text.y = element_text(family = "serif"),
plot.title = element_text(family = "serif"),
panel.background = element_rect(fill = "gray100",
colour = "aquamarine4"), plot.background = element_rect(fill = "white",
colour = NA, linetype = "solid"),
legend.position = "none") +labs(x = NULL, y = NULL, fill = NULL)+coord_fixed(ratio =4 / 1)

#2.2 饼图--------------------
bt=na.omit(read_xls("MGWR_CU.xls",sheet="77"))#数据清洗
# Create data
library(ggplot2)
ggplot(bt, aes(x="", y=bt[[2]], fill=bt[[1]])) +
geom_bar(stat="count", width=1) +
coord_polar(theta="y")+
theme_void() + # 去除背景和坐标轴
ggtitle("扇形图") # 添加标题

ggplot(bt, aes(class))+
geom_bar(aes(fill=bt[[1]]))+
coord_polar(theta = "x")

# create color palette:
library(RColorBrewer)
coul <- brewer.pal(3, "Pastel2")

# Transform this data in %
data_percentage <- apply(data, 2, function(x){x*100/sum(x,na.rm=T)})

# Make a stacked barplot--> it will be in %!
barplot(data_percentage, col=coul , border="white", xlab="group")