#代码参考《R语言数据可视化之美-专业图表绘制指南》例4.2.1,主要代码链接参见
#https://github.com/EasyChart/Beautiful-Visualization-with-R,
#在添加了详细的注释基础上,对一处BUG代码进行了勘误
rm(list = ls())#清空工作区,可在括号内加装变量等定向清除某个变量等
library(lattice)
library(gridExtra)
library(reshape2)
library(RColorBrewer)
colormap<-colorRampPalette(rev(brewer.pal(11,'RdYlGn')))(100)
path<-readline()#括号不写内容,默认读取下一行
D:\Program Files\RStudio\Beautiful-Visualization-with-R-master\Beautiful-Visualization-with-R-master\第4章 数据关系型图表
#设定工作目录,输出图片等需要到的目录
Mypath<- gsub("\\\\", "/",path)
setwd(Mypath)
getwd()
#--------------------------------------------------澶氶」寮忔嫙鍚?------------------------------------------
mydata <- read.csv("Surface_Data.csv", sep= ",", header=T)
#澶氶」寮忔嫙鍚坺=f(x,y)=a+bx+cy+dxx+eyy
x <- mydata$x
y <- mydata$y
z <- mydata$z
x2<-x*x
y2<-y*y
poly_z <- lm(z ~ x + y +x2+y2)
print(poly_z)
#随便建立的模型,实际中应该是自己建立合适的模型
N<-30
xmar <- seq(min(x),max(x),(max(x)-min(x))/N)#max(x)=25,min(x)=8,(max(x)-min(x))/N=0.566667
ymar <- seq(min(y),max(y),(max(y)-min(y))/N)#max(y)=36,min(x)=10,(max(x)-min(x))/N=0.566667
#该指令的意思是最大值到最小值之间等分成N份
Grid_xy<-expand.grid(list(x=xmar,y=ymar))#按x,y的等分,建立起基础网格点
Grid_xy$x2<-Grid_xy$x*Grid_xy$x
Grid_xy$y2<-Grid_xy$y*Grid_xy$y#基础网格点坐标再平方
Grid_z <- predict.lm(poly_z, newdata=Grid_xy)#利用已经建立的模型,输出其他点的预测值
df<-data.frame(matrix(Grid_z, length(xmar), length(ymar)))#length(xmar)=length(ymar)=N+1
#将预测值与(N+1)*(N+1)的方阵相互联系,组成三维矩阵;其x,y为1:N+1以1为单位递增的数列
colnames(df)<-xmar#通过改变列名称的方法让矩阵的x重定义到xmar数列上
df$x<-ymar#通赋值的方法让矩阵增加一列x,其值为ymar
melt_df<-melt(df,id.vars='x', variable.name ="y",value.name = "z")#把三维矩阵的三个维度信息用x,y,z表示,并整理为三列
#df为31*32=992个数据,如果刨除x列,则有961个(三维,其行列信息也是)数据;melt_df有961*3个数据,
#显然是把每个数据对应的行列信息都拓展出来了
melt_df$y<-as.numeric(melt_df$y)#这里melt_df$y从因子变量变为数值变量,但重新失去了原有ymar数列信息
#trellis.par.set("axis.line",list(col=NA,lty=1,lwd=1)) # Removes the border of the plot if you want
melt_df$y<-(melt_df$y-1)*(max(x)-min(x))/N+min(x)#勘误,这样才保留原始的y轴信息
surface_plot1 <- wireframe(z ~ y*x, data=melt_df,
xlab = "as.numeric转换后校正",
ylab = "as.numeric转换后校正",
zlab="Power (KW)",
zlim=c(20,180),
drape = TRUE,
colorkey = TRUE,
scales = list(arrows=FALSE),
light.source = c(10,0,10),
col.regions = colorRampPalette(rev(brewer.pal(11,'RdYlGn')))(100),
screen = list(z = -60, x = -60)
)
surface_plot1