R语言曲面拟合代码详细分析（1）

#代码参考《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