R语言ggplot2 绘图，曲线拟合，双坐标轴展示

在进行数据处理时，经常会有一些数据丢失，我们可以使用不同得拟合方法将数据进行拟合，本文以地表温度、降雨数据为例，用ggplot2进行绘图并进行数据的拟合。用双坐标轴对结果进行展示。

首先是设置工作空间及数据的导入

rm(list = ls())
setwd("E:/test")
library(ggplot2)
library(dplyr)
library(cowplot)

ls1=read.csv("54511.csv",header = T)
ls2=read.csv("Beijing.csv",header = T)

因为后期需要用到拟合，所以设置一个Π值为3.14，建一个有12个值得数组用于存放每月得气温。

用for循环依次取出月这一列数值等于1-12的所有行，进行均值计算，均值取出来存到事先建好的数组里。再把温度数据框中加一个X值表示1-12个月。


Pi=3.14

E_tem=array(1:12)

for( i in 1:12)
{
  Mdata = ls1[ls1$Mon==i,]
  E_tem[i] = mean(Mdata$TEM)
}

E_tem=data.frame(E_tem) #均气温
E_tem=data.frame(x=c(1:12),y=E_tem$E_tem)

把地表温度的值取出，然后进行拟合，此次选用正弦+余弦拟合，y=sin(Ωpix)+cos(Ωpix),其中Ω的度数需要根据具体数据进行调整。

e_TEM=ls2$RASTERVALU  #地表温度

e_TEM=data.frame(x=c(1:12),y=e_TEM)

lis1=lm(data=E_tem,formula ='y~(sin(1/6*pi*x)+cos(1/6*Pi*x))')  #气温拟合
lis2=lm(data=e_TEM,formula ='y~(sin(1/6*pi*x)+cos(1/6*Pi*x))')  #地表温度拟合

拟合完成后也可用summary函数查看其各个拟合值得参数，这里是拼接了一下拟合后得公式

MODEL_Ta=summary(lis1)
names(MODEL_Ta)
coefficients_R2 = MODEL_Ta $r.squared
coefficients_Ta = MODEL_Ta $coef[,"Estimate"]
coefficients_R2
coefficients_Ta
names(coefficients_Ta)
Intercept_Ta = coefficients_Ta[[1]]
print(Intercept_Ta)

Sin_Ta= coefficients_Ta[[2]]
print(Sin_Ta)

Cos_Ta=coefficients_Ta[[3]]
print(Cos_Ta)

formula_Ta=paste("y=",Sin_Ta,"*sin(1/6*pi*x)",Cos_Ta,"*cos(1/6*Pi*x)","+",Intercept_Ta)
print(formula_Ta)

整理数据进行绘图，sec.axis 双坐标轴绘制，关于图例须注意的是如果写在括号里，只写能表明循序的值即可，不用标注颜色，在使用scale_color_manual方法时，会根据你标注的ABC顺序设置颜色，顺序的依据为英文字母前后顺序。


df1=data.frame(x=c(1:12),y=lis2$fitted.values-lis1$fitted.values) #这里时取了两个拟合的差值

df=data.frame(x=c(1:12),a=E_tem$y,b=e_TEM$y,c=df1$y)

tiff(file="Surface Temperature and air Temperature.tif",height = 400*8,width = 1000*8,res = 100*8,units = "px",compression = "lzw")                             ##输出格式

f1=ggplot()+
  
  geom_point(size=1.5,data=df,mapping = aes(x=x,y=a,color="A"))+
  geom_line(data=df,mapping = aes(x=x,y=a,color="A"))+
  geom_smooth(data = df,mapping = aes(x=x,y=a),method = "lm",formula ='y~(sin(1/6*pi*x)+cos(1/6*Pi*x))',color="red")+          #画均温及拟合
  
  geom_point(size=1.5,data=df,mapping = aes(x=x,y=b,color="B"))+
  geom_line(data=df,mapping = aes(x=x,y=b,color="B"))+
  geom_smooth(data = df,mapping = aes(x=x,y=b),method = "lm",formula ='y~(sin(1/6*pi*x)+cos(1/6*Pi*x))',color="orange")+       #地表温度及拟合
  
  geom_line(data=df,mapping = aes(x=x,y=c,color="C"))+                                                                         #ΔTs-a
  
  scale_x_continuous(name='Month',breaks = c(1:12),labels=c("1","2","3","4","5","6","7","8","9","10","11","12"))+
  scale_y_continuous(
    name='Temperature（℃）',
    sec.axis = sec_axis(~.*1,breaks = seq(0,30,10),
                        name = "ΔTs-a（℃）",
                        labels=c("0","1","2","3"))      
  )+                                                                                                                            #双坐标轴
  
  geom_hline(yintercept = 0,color="gray")+                                                                                      #垂线
  labs(title="Beijing")+                                                                                                        #图名
  scale_color_manual(name="",values = c("green","blue","yellow"),labels=c("Ta","Ts","ΔTs-a"))+                                  #图例
  theme_bw()

print(f1)
dev.off()

完整代码展示：对于拼接，我选择的是cowplot包比较实用一句代码就能搞定，拼接在文末

rm(list = ls())
setwd("E:/test")
library(ggplot2)
library(dplyr)
library(cowplot)

ls1=read.csv("54511.csv",header = T)
ls2=read.csv("Beijing.csv",header = T)


Pi=3.14

E_tem=array(1:12)

for( i in 1:12)
{
  Mdata = ls1[ls1$Mon==i,]
  E_tem[i] = mean(Mdata$TEM)
}

E_tem=data.frame(E_tem) #均气温
E_tem=data.frame(x=c(1:12),y=E_tem$E_tem)

e_TEM=ls2$RASTERVALU  #地表温度

e_TEM=data.frame(x=c(1:12),y=e_TEM)

lis1=lm(data=E_tem,formula ='y~(sin(1/6*pi*x)+cos(1/6*Pi*x))')  #气温拟合
lis2=lm(data=e_TEM,formula ='y~(sin(1/6*pi*x)+cos(1/6*Pi*x))')  #地表温度拟合

#公式拼接

MODEL_Ta=summary(lis1)
names(MODEL_Ta)
coefficients_R2 = MODEL_Ta $r.squared
coefficients_Ta = MODEL_Ta $coef[,"Estimate"]
coefficients_R2
coefficients_Ta
names(coefficients_Ta)
Intercept_Ta = coefficients_Ta[[1]]
print(Intercept_Ta)

Sin_Ta= coefficients_Ta[[2]]
print(Sin_Ta)

Cos_Ta=coefficients_Ta[[3]]
print(Cos_Ta)

formula_Ta=paste("y=",Sin_Ta,"*sin(1/6*pi*x)",Cos_Ta,"*cos(1/6*Pi*x)","+",Intercept_Ta)
print(formula_Ta)
maxPoint_Ta=(4*Sin_Ta*sin(1/6*pi)*Intercept_Ta-(Cos_Ta*cos(1/6*pi))*Cos_Ta*cos(1/6*pi))/4*Sin_Ta*sin(1/6*pi)
print(maxPoint_Ta)

MODEL_Ts=summary(lis2)
names(MODEL_Ts)
coefficients_R2 = MODEL_Ts $r.squared
coefficients_Ts = MODEL_Ts $coef[,"Estimate"]
coefficients_R2
coefficients_Ts
names(coefficients_Ts)
Intercept_Ts = coefficients_Ts[[1]]
print(Intercept_Ts)

Sin_Ts= coefficients_Ts[[2]]
print(Sin_Ts)

Cos_Ts=coefficients_Ts[[3]]
print(Cos_Ts)

formula_Ts=paste("y=",Sin_Ts,"*sin(1/6*pi*x)",Cos_Ts,"*cos(1/6*Pi*x)","+",Intercept_Ts)
print(formula_Ts)

df1=data.frame(x=c(1:12),y=lis2$fitted.values-lis1$fitted.values)

df=data.frame(x=c(1:12),a=E_tem$y,b=e_TEM$y,c=df1$y)

#y= -7.1359 * sin(1/6 * pi * x) -14.2869 * cos(1/6 * Pi * x)  + 13.6185 

tiff(file="Surface Temperature and air Temperature.tif",height = 400*8,width = 1000*8,res = 100*8,units = "px",compression = "lzw")                             ##输出格式

f1=ggplot()+
  
  geom_point(size=1.5,data=df,mapping = aes(x=x,y=a,color="A"))+
  geom_line(data=df,mapping = aes(x=x,y=a,color="A"))+
  geom_smooth(data = df,mapping = aes(x=x,y=a),method = "lm",formula ='y~(sin(1/6*pi*x)+cos(1/6*Pi*x))',color="red")+          #画均温及拟合
  
  geom_point(size=1.5,data=df,mapping = aes(x=x,y=b,color="B"))+
  geom_line(data=df,mapping = aes(x=x,y=b,color="B"))+
  geom_smooth(data = df,mapping = aes(x=x,y=b),method = "lm",formula ='y~(sin(1/6*pi*x)+cos(1/6*Pi*x))',color="orange")+       #地表温度及拟合
  
  geom_line(data=df,mapping = aes(x=x,y=c,color="C"))+                                                                         #ΔTs-a
  
  scale_x_continuous(name='Month',breaks = c(1:12),labels=c("1","2","3","4","5","6","7","8","9","10","11","12"))+
  scale_y_continuous(
    name='Temperature（℃）',
    sec.axis = sec_axis(~.*1,breaks = seq(0,30,10),
                        name = "ΔTs-a（℃）",
                        labels=c("0","1","2","3"))      
  )+                                                                                                                            #双坐标轴
  
  geom_hline(yintercept = 0,color="gray")+                                                                                      #垂线
  labs(title="Beijing")+                                                                                                        #图名
  scale_color_manual(name="",values = c("green","blue","yellow"),labels=c("Ta","Ts","ΔTs-a"))+                                  #图例
  theme_bw()

print(f1)
dev.off()

#以下是绘制日均温代码

tiff(file="Daily mean temperature and Precipitation.tif",height = 400*8,width = 1000*8,res = 100*8,units = "px",compression = "lzw") 

day=array(data =c(31,28,31,30,31,30,31,31,30,31,30,31))
z=0
m_tem=array(1:365)
m_rain=array(1:365)
for (i in 1:12)
{
  df2=ls1[ls1$Mon==i,]
  q=day[i]
  for (j in 1:q) 
  {
    
    df3=df2[df2$Day==j,]
    z=z+1
    
    m_tem[z]=mean(df3$TEM)
    
  }
  
}

c=0

for (a in 1:12)
{
  d1=ls1[ls1$Mon==a,]
  p=day[a]
  for (b in 1:p) 
  {
    
    d2=d1[d1$Day==b,]
    c=c+1
    m_rain[c]=mean(d2$PRE_1h)
    
  }
  
}



m_tem=data.frame(x=c(1:365),y=m_tem)
m_rain=data.frame(x=1:365,y=m_rain)

#fft_rain=fourierFit(m_rain,3)  #傅里叶变化，没有成功

f2=ggplot()+
  
  geom_point(data = m_tem,mapping = aes(x=x,y=y,color="A"))+
  geom_line(data = m_tem,mapping = aes(x=x,y=y,color="A"))+
  
  
  geom_point(data = m_rain,mapping = aes(x=x,y=y*10,color="B"))+
  geom_line(data = m_rain,mapping = aes(x=x,y=y*10,color="B"))+
  
  
  scale_x_continuous(name="day")+
  scale_y_continuous(name='Temperature(℃)',sec.axis = sec_axis(~.*1,breaks = seq(0,30,10),
                                                               name = "pre(mm)",
                                                               labels=c("0","1","2","3")))+#双坐标轴
  scale_color_manual(name="",values = c("green","blue","yellow"),labels=c("Ta","Precipitation"))+ 
  labs(title="Beijing")+     
  
  theme_bw()
print(f2)
dev.off()


tiff(file="Beijing.tif",height = 400*8,width = 1000*8,res = 100*8,units = "px",compression = "lzw")
p5 <- cowplot::plot_grid(f1,f2,nrow = 1, labels = LETTERS[1:2])                            #将p1-p4四幅图组合成一幅图，按照两行两列排列，
print(p5)                                                                           #标签分别为A、B、C、D。（LETTERS[1:4] 意为提取26个大写英文字母的前四个:A、B、C、D）

dev.off()

下面展示几张结果图

这张图中可以明显看出九月份数据是缺失的，但是通过拟合自动补全了。

这是对点数据的拟合。

cowplot拼接图

原文链接：https://blog.csdn.net/qq_42940285/article/details/121014965

你可能也喜欢