此文也发表在我的公众号，欢迎关注：magicaircode
01 前言
绘图应该是大家在科研和工作中最常用的R语言功能之一了。
鉴于此，我想开通一个系列《复现》。
该系列旨在通过R语言代码复现一些好看或者常用的图表，给你一些绘图的灵感、思路和方法。
图表来自大气化学或空气质量相关的科研文章或报告。
如果你有喜欢的和想要复现的图表，也欢迎后台私信我，也许我就帮你复现出来了。
今天的《复现》第一期的图来自我近期发表的文章中的甲醛速率日变化的图。
02 图片概况
这是一幅组合图。代表同一天两个地点的甲醛化学速率的日变化图。为便于讲解，我让左右图数值一致了。
图例中色块代表了甲醛生成和去除了几种途径，黑色实线表示净生成速率。图中y轴表示反应速率（生成为正，去除为负），x轴表示时间。
(Chen et al., 2022)
03 绘图思路
初步分析可知，绘图步骤是：绘制面积堆积图层（生成和去除速率），绘制线形图层（净生成速率），两个图层叠加，得到第一个图。再绘制另一个图，组合两个图。
04 数据概况
下方显示了绘图示例数据前五行。第一列表示时间，第二至三列表示甲醛的去除速率，第四至八列表示甲醛的去除速率。我们将这个表命名为hchobg。
datetime L1 L2 P1 P2 P3 P4 P5
1 2017-07-09 00:00:00 0 -0.01 0.02 0.03 0.00 0.18 0.03
2 2017-07-09 01:00:00 0 -0.01 0.02 0.02 0.00 0.15 0.03
3 2017-07-09 02:00:00 0 -0.01 0.02 0.02 0.00 0.15 0.03
4 2017-07-09 03:00:00 0 -0.01 0.02 0.02 0.00 0.15 0.03
5 2017-07-09 04:00:00 0 -0.01 0.02 0.02 0.00 0.13 0.0205 数据预处理
先计算好净生成速率，即每一行求和（生成速率和去除速率求和）。净生成速率的表命名为hcho_netp。
hcho_netp=data.frame(datetime=hchobg[,1],
netrate=rowSums(hchobg[,-1]))现在我们有两个表：hchobg和hcho_netp。
绘图需要三列表的格式（第一列时间，第二列变量数值，第三列变量类别）。在本例中，第二列即生成和去除速率的数值，第三列即生成和去除途径的代号。因此，我们需要先对数据进行格式转换。
这里使用了reshape2包的melt函数，一键转化生成和去除速率表hchobg。转化后的命名为hchobg_melt：
library(reshape2)
hchobg_melt=melt(hchobg,id.vars="datetime")此外，因为绘制面积堆积图时，各个类别（第三列）需要是因子型数据，所以要对第三列从文本型转化为因子型：
hchobg_melt$variable=factor(hchobg_melt$variable,
levels=c('L1', 'L2', 'P1', 'P2', 'P3', 'P4', 'P5'))对于净生成速率表hcho_netp，本身只有两列，因此我们重新制作一下，加上第三列表示类别即可，新表命名为hcho_rate：
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和hcho_rate。他们的前五行分别如下：
hchobg_melt（生成和去除速率的时间序列表）：
datetime variable value
1 2017-07-09 00:00:00 L1 0
2 2017-07-09 01:00:00 L1 0
3 2017-07-09 02:00:00 L1 0
4 2017-07-09 03:00:00 L1 0
5 2017-07-09 04:00:00 L1 0
6 2017-07-09 05:00:00 L1 0hcho_rate（净生成速率的时间序列表）：
datetime rate gtype
1 2017-07-09 00:00:00 0.25 Net production
2 2017-07-09 01:00:00 0.21 Net production
3 2017-07-09 02:00:00 0.21 Net production
4 2017-07-09 03:00:00 0.21 Net production
5 2017-07-09 04:00:00 0.18 Net production
6 2017-07-09 05:00:00 1.12 Net production06 绘图
先画出草图（别忘了载入ggplot2包）。用ggplot()创建空白图，然后采用了geom_area()和geom_line()两个函数分别绘制面积堆积图和实线图，然后两个图层叠加。
library(ggplot2)
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)
有点粗糙吧。别急。我们先整理一下y轴。改一下时间的显示方式。
p=p+scale_x_datetime(breaks = "3 hour", date_labels = "%H",minor_breaks="1 hours",expand=c(0,0))
是不是舒服了一些。我们再改写x轴和y轴的轴名，移除图例的标题。
p=p+labs(x="2017-07-09", y="Rate (ppbv/h)", fill=NULL, linetype=NULL)
图例的文字需要改写成具体内容而不是代号。然后我们还需要换一套更加有区分度的舒适的颜色。
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。
x1=p+annotate("text", x = as.POSIXct("2017-07-09 23:00:00",tz="GMT"), y = Inf, hjust = 2, vjust = 2, label = "(a)")
然后接下来为了演示两个图的组合，我们依样画葫芦，再用同一个图，加另一个子图编号，生成另一个图。生成的图命名为x2。
x2=p+annotate("text", x = as.POSIXct("2017-07-09 23:00:00",tz="GMT"), y = Inf, hjust = 2, vjust = 2, label = "(b)")
然后就是合并两个图了。依旧是采用我特别喜欢的patchwork包。通过这个包，我们只要输入两个图名称，并采用符号连接，例如x1/x2，这两个图就会左右分布了。
然后我们在后面加上这两句代码，以合并图例，并缩小空白。
library(patchwork)
x1/x2+plot_layout(guides = 'collect',ncol = 2)&theme(legend.margin=unit(0, 'cm'))
这样，我们就复现出了本文开始的图了。
07 尾声
如果你想试试，可以后台私信索取代码。这个系列刚刚开始，欢迎提供意见和建议。也欢迎提供提供你感兴趣的或者想复现的图表。撰文不易，欢迎点击“赞”、“在看”和“分享”，非常感谢。
08 本期图表来源
此文也发表在我的公众号，欢迎关注：magicaircode