Data input

library(readxl)
data <- read_excel("D:/plots for PGG/20220307/data.xlsx")
head(data)
## # A tibble: 6 x 8
##    year `N_(mg/L)`   N_SD `P_(mg/L)`     P_SD `Chla_(ug/L)` Chla_SD    EI
##   <dbl>      <dbl>  <dbl>      <dbl>    <dbl>         <dbl>   <dbl> <dbl>
## 1  1999      0.045 0.0225    0.0208  0.0104           NA     NA     0.125
## 2  2000      0.154 0.0771    0.00488 0.00244          NA     NA     0.05 
## 3  2001      0.175 0.0877    0.0141  0.00705          NA     NA     0.05 
## 4  2002      0.361 0.180     0.00868 0.00434          NA     NA     1.36 
## 5  2003      0.232 0.118     0.0122  0.00199           1.37   0.681 0.2  
## 6  2004      0.208 0.156     0.00418 0.000884          2.57   0.398 0.170

Plot ! The key is geom_ribbon !

library(ggplot2)
ggplot(data, aes(x = year))+
    geom_line(aes(y = 80*`P_(mg/L)`, color = "PO4 (mg/L)"), size = 1)+
    geom_ribbon(aes(x = year, ymin=(`P_(mg/L)`-P_SD)*80, ymax=(`P_(mg/L)`+P_SD)*80),fill='orange',alpha = 0.15)+
    geom_line(aes(y = `Chla_(ug/L)`, color = "Chl a (ug/L)"), size = 1 )+
    geom_ribbon(aes(x = year, ymin=`Chla_(ug/L)`-Chla_SD, ymax=`Chla_(ug/L)`+Chla_SD),fill='red',alpha = 0.1)+
    geom_line(aes(y = EI, color = "EI"), size = 1)+
    geom_line(aes(y = `N_(mg/L)`*15, color = "DIN (mg/L)"), size = 1)+
    geom_ribbon(aes(x = year, ymin=(`N_(mg/L)`-N_SD*0.5)*15, ymax=(`N_(mg/L)`+N_SD*0.5)*15),fill='black',alpha = 0.1)+theme_bw()+
    scale_y_continuous(sec.axis = sec_axis(~./15, name = "DIN (mg/L)", breaks = c(0, 0.1, 0.2, 0.3, 0.4, 0.5)))+
    #~./15表示次轴与主轴的关系是:次轴等于主轴除以15
    #name表示次y轴的名称
    scale_x_continuous(breaks = data$year, labels = c('1983-1984', '1990', '1996-1997', '1999', 2003:2018))+
    scale_color_manual(values = c("black", "red", "blue", "orange"), limits = c("DIN (mg/L)", "Chl a (ug/L)", "EI", "PO4 (mg/L)"))+
    theme(axis.title = element_text(size = 18, face = 'bold'), axis.text = element_text(size = 16, colour = 'black'), legend.position = 'left', legend.text = element_text(size = 16, colour = 'black'), axis.title.y.right = element_text(angle = 90, hjust =0.5), axis.text.x = element_text(size = 12, angle = 320, vjust = 0.5, hjust = 0.2), legend.title = element_blank(), axis.ticks = element_line(size = 0.6), axis.ticks.length = unit(1.5, 'mm')) + xlab('Year') + ylab('Chl a (ug/L)') + theme(panel.border = element_rect(size = 1.2))
## Warning: Removed 7 row(s) containing missing values (geom_path).