library(readxl)
library(reshape2)
library(ggplot2)
library(lubridate)
library(patchwork)
hchobg=structure(list(datetime = structure(c(1499558400, 1499562000,
1499565600, 1499569200, 1499572800, 1499576400, 1499580000, 1499583600,
1499587200, 1499590800, 1499594400, 1499598000, 1499601600, 1499605200,
1499608800, 1499612400, 1499616000, 1499619600, 1499623200, 1499626800,
1499630400, 1499634000, 1499637600, 1499641200), class = c("POSIXct",
"POSIXt"), tzone = "UTC"), L1 = c(0, 0, 0, 0, 0, 0, -0.144, -0.639,
-1.314, -1.395, -3.096, -3.42, -2.16, -1.746, -2.358, -2.367,
-1.152, -0.324, -0.081, 0, 0, 0, 0, 0), L2 = c(-0.01, -0.01,
-0.01, -0.01, -0.01, -0.072, -0.189, -0.576, -0.99, -0.864, -1.476,
-1.179, -0.558, -0.504, -0.729, -0.72, -0.333, -0.126, -0.045,
-0.009, -0.009, 0, 0, 0), P1 = c(0.02, 0.02, 0.02, 0.02, 0.02,
0.45, 1.269, 2.232, 2.916, 2.718, 4.23, 4.059, 2.493, 2.169,
2.826, 2.961, 1.728, 0.927, 0.405, 0.063, 0.018, 0.018, 0.018,
0.018), P2 = c(0.03, 0.02, 0.02, 0.02, 0.02, 0.261, 0.936, 1.881,
2.646, 2.547, 3.492, 3.141, 1.863, 1.656, 2.142, 2.025, 1.017,
0.558, 0.225, 0.045, 0.027, 0.018, 0.018, 0.018), P3 = c(0, 0,
0, 0, 0, 0.009, 0.054, 0.099, 0.117, 0.099, 0.153, 0.135, 0.072,
0.063, 0.09, 0.099, 0.054, 0.027, 0.009, 0, 0, 0, 0, 0), P4 = c(0.18,
0.15, 0.15, 0.15, 0.13, 0.162, 0.18, 0.243, 0.333, 0.414, 0.54,
0.711, 0.837, 0.729, 0.666, 0.594, 0.513, 0.414, 0.333, 0.252,
0.171, 0.135, 0.117, 0.117), P5 = c(0.03, 0.03, 0.03, 0.03, 0.02,
0.198, 0.342, 0.729, 1.044, 0.963, 1.521, 1.494, 0.945, 0.828,
1.143, 1.179, 0.684, 0.369, 0.162, 0.045, 0.036, 0.036, 0.027,
0.027)), row.names = 385:408, class = "data.frame")
#观测和净生成组成数据框
#计算净生成速率
hcho_netp=data.frame(datetime=hchobg[,1], netrate=rowSums(hchobg[,-1]))
#hcho_netp$datetime<-as.POSIXlt(floor(as.double(hcho_netp$datetime) / 60) * 60, origin = '1970-01-01',tz='GMT')
hcho_rate=data.frame(datetime=c(hcho_netp$datetime),rate=c(hcho_netp[,2]),gtype=c(rep("Net production", length(hcho_netp$datetime))))
#速率转化成长表格式
hchobg_melt=melt(hchobg,id.vars="datetime")
#将反应途径编号转化为因子型数据
hchobg_melt$variable=factor(hchobg_melt$variable,levels=c('L1', 'L2', 'P1', 'P2', 'P3', 'P4', 'P5'))
#先画出草图
p=ggplot() +
geom_area(data = hchobg_melt, aes(x=datetime, y=value, fill=variable),stat = "identity")+
geom_line(data=hcho_rate, aes(x=datetime, y=rate, linetype=gtype),size=1.2)
#x轴间隔和格式
p=p+scale_x_datetime(breaks = "3 hour", date_labels = "%H",minor_breaks="1 hours",expand=c(0,0))
#改写x和y轴名，移除legend标题
p=p+labs(x="2017-07-09", y="Rate (ppbv/h)", fill=NULL, linetype=NULL)
#HCHO颜色和图例
hcho_pal=scale_fill_manual(
breaks = c('P1', 'P2', 'P3', 'P4', 'P5', 'L1', 'L2'),
labels = c(
'L1'=bquote("HCHO photoysis"),
'L2'=bquote("HCHO + OH"),
'P1'=bquote(O[2]~"+"~CH[3]*O),
'P2'=bquote(O[2]~"+"~RO~"(from alkenes)"),
'P3'=bquote(O[2]~"+"~RO~"(from other VOCs)"),
'P4'=bquote(O[3]~"+"~VOC),
'P5'=bquote("Other reactions")
),values = c(
'L1'="#9254de",
'L2'="#2f54eb",
'P1'="#fa541c",
'P2'="#faad14",
'P3'="#a0d911",
'P4'="#52c41a",
'P5'="#fadb14"##13c2c2
))
p=p+hcho_pal
#主题风格调整
p=p+theme_bw()
x1=p+annotate("text", x = as.POSIXct("2017-07-09 23:00:00",tz="GMT"), y = Inf, hjust = 2, vjust = 2, label = "(a)")
x2=p+annotate("text", x = as.POSIXct("2017-07-09 23:00:00",tz="GMT"), y = Inf, hjust = 2, vjust = 2, label = "(b)")
x1/x2+plot_layout(guides = 'collect',ncol = 2)&theme(legend.margin=unit(0, 'cm'))
