简介
在查阅文献的过程中，看到了几幅非常不错的出版图，今天就跟着小编一起学习下，他们是怎么使用 R 绘制出来的。
今天主要介绍 第四幅图（D) —— 实现双 Y 轴，并且添加坐标轴的微小刻度线。这个图在科研绘图中较为常用，例如：将算法的收敛情况和计算所耗时间同时绘制。
前三幅图的详细代码介绍可见：基于 R 语言的科研论文绘图技巧详解（3）基于 R 语言的科研论文绘图技巧详解（2）基于 R 语言的科研论文绘图技巧详解（1）。后面几幅图会一一介绍，读者在学习过程中，可以将内部学到的知识点应用到自己的图形绘制中。推文已经将主要知识点进行罗列，更有利于读者学习和查阅。
那我们来看看，作者是怎么实现这个功能的吧，本文知识点较多，大家耐心学习，建议自己实践。对应代码、数据可在 GitHub - marco-meer/scifig_plot_examples_R: Scientific publication figure plotting examples with R[1] 中找到。
主要知识点
- 实现双 Y 轴；
- 学会修改坐标轴为对数尺度；
- 添加坐标轴的微小刻度线。
绘图
加载包
首先加载一些需要使用到的包。
library(ggplot2) # Grammar of graphics
设置主题
接下来，为了方便起见，作者在绘图前设置好了主题，并将该函数命名为 my_theme。
这一部分在第一篇推文 基于 R 语言的科研论文绘图技巧详解（1）给出，代码将在文末中完整代码给出。
手动修改大部分面板，具体可以参考本篇文章[2]。或者观看我在 B 站发布的《R 语言可视化教程》，里面也有一些简单主题设置介绍。
导入数据
首先使用 read.csv() 导入数据，其中一个数据前几行如下所示。
data_D1 = read.csv("./data_D1.csv") data_D2 = read.csv("./data_D2.csv") # format data to ggplot's liking data_D = data.frame("width"=c(data_D1$width,data_D2$width), "unit"=c(rep("shear_stress",nrow(data_D1)), rep("velocity",nrow(data_D2))), "value"=c(data_D1$shear_stress,data_D2$velocity) ) head(data_D) # width unit value # 1 20.0 shear_stress 0.00174000 # 2 6.0 shear_stress 0.01622000 # 3 2.0 shear_stress 0.16065000 # 4 1.0 shear_stress 0.65696000 # 5 4.0 shear_stress 0.02312000 # 6 0.6 velocity 0.01726544
这里得到的数据，一共有三列，两个数据集的值在 value 中，width 放了两个数据集各自的width， unit 为离散数据。
绘图步骤详解
关键在于如何构建双 Y 轴，下面来看看作者是怎么设置的吧。
绘制单轴
首先，处理下第一个线性图所需要的数据，一共是两列。
curve_D1 = data.frame(width=data_D1$width, shear_stress=33.28/(pi*18*data_D1$width^2)) curve_D1 # width shear_stress # 1 20 0.001471299 # 2 6 0.016347767 # 3 2 0.147129903 # 4 1 0.588519612 # 5 4 0.036782476
这里绘制，小编带大家一步步解释，尤其注意作者的思想。
顺便提一下，很多公众号只是给了搬运的代码，但是并没有解释其中的奥秘。这对于初学者而言，就很难理解了。小编这个系列就是带大家一起学习作者画图思路。学会融会贯通，用到自己的科研绘图中。
先简单绘制出线性图，可以看到：在 x 轴附近， y 轴下降的很快。
ggplot(data=data_D1, aes(x=width,y=shear_stress)) + geom_point(fill="red",size=3,pch=22) + geom_line(data=curve_D1)
我们可以对其做下对数变换，使得数据呈线性形式。使用 scale_x_log10() 和 scale_y_log10() 对刻度进行对数变换。内部参数这里不做解释，大家看着修改，就知道内部含义了。
scale_x_log10(expand=c(0,0), breaks=c(0.5,1,2,5,10,20,50), labels=c(0.5,1,2,5,10,20,50), limits=c(0.5,50)) + scale_y_log10(expand = c(0, 0), labels = trans_format('log10', math_format(10^.x)), breaks=c(0.001,0.01,0.1,1), limits = c(0.001,1))
之后，我们对主题进行调整。尤其注意这里的 annotation_logticks(sides = "l")，这个代码可以增加对数坐标轴的标记（左边位置）。
annotation_logticks(sides = "l") + theme_void()+ theme( line = element_blank(), # exclude everything outside axes bc it messes with positioning of grob in panel_D text = element_blank(), title = element_blank(), axis.line.y = element_line(colour = "black") ) + ylab("shear stress (Pa)")
添加到另一张图形中
之后，将前面的图添加到另一张线性图中。首先把另一张图绘制出来：
ggplot(data=data_D2, aes(x=width,y=velocity))+ geom_point(fill="blue",size=3,pch=21)
之后使用 annotation_custom(ggplotGrob(panel_D1)) 将前面那幅图添加到该图中。此时结果如下：
注意:
annotation_custom()是一个特殊的集合对象，用于静态注释。注释不会影响缩放。
这时，恭喜你两幅图已经合并啦！但是存在几个问题：
- 两幅图的 Y 轴重复了。这时候使用
scale_y_continuous()将原图的 Y 轴位置往右放置（position = "right"）。 - 但是变换完之后，左边标签没有，而左边的 Y 轴就是第一幅图得到的结果，我们需要添加缺失的标签。处理方式为：使用
sec.axis = sec_axis(~., name = "shear stress (Pa)",breaks=rescale(c(-3,-2,-1,0), to = c(0,1)),labels = c(expression("10"^"-3","10"^"-2", "10"^"-1", "10"^"0"))))手动添加坐标标签和数值。 - 两幅图的 x 轴不一致，使用
scale_x_log10()修改结果。 - 使用
annotation_logticks(sides = "b")添加 x 轴的 ticks。
所以，如果读者想使用这里的代码，还是需要一定能力看懂这些代码的～ 小编已经教会你们啦，以后要学会根据自己的结果更改代码噢～
scale_y_continuous(expand = c(0,0), breaks = seq(0,1,0.1), limits = c(0,1), position = "right", sec.axis = sec_axis(~., name = "shear stress (Pa)",breaks=rescale(c(-3,-2,-1,0), to = c(0,1)),labels = c(expression("10"^"-3","10"^"-2", "10"^"-1", "10"^"0")))) + scale_x_log10(expand=c(0,0), breaks=c(0.5,1,2,5,10,20,50), labels=c(0.5,1,2,5,10,20,50), limits=c(0.5,50)) + annotation_logticks(sides = "b")
后面就是添加一些线段，文字以及修改主题啦～前面几个推送都已经介绍过了，这里不做解释了。
my_theme() + geom_line(color="blue") + geom_vline(xintercept = 1.1, linetype="dashed") + geom_hline(yintercept = 0.9, linetype="dashed") + theme( plot.margin = unit(c(0.1,0,0,0.5), "cm"), # to match other panels axis.title.y = element_text(margin = margin(r=1)), axis.text.y = element_text(margin = margin(r=6)), axis.text.y.right = element_text(margin = margin(l=7)), axis.title.y.right = element_text(angle = 90) ) + xlab(expression(lumen~width~(mu*m))) + ylab("relative flow velocity") + annotate(geom = "text",x =6 ,y =0.85 ,label = "tau == 0.5~Pa",parse=T) + annotate(geom = "text",x =1.4 ,y =0.4 ,label = "b == 1.1*mu*m",parse=T,angle=90) + annotate(geom = "text",x =5.6 ,y =0.6 ,label = "tau == frac(4*mu*Q,pi*a*b^2)",parse=T)
完整代码
library(ggplot2) base_size = 12 my_theme <- function() { theme( aspect.ratio = 1, axis.line =element_line(colour = "black"), # shift axis text closer to axis bc ticks are facing inwards axis.text.x = element_text(size = base_size*0.8, color = "black", lineheight = 0.9, margin=unit(c(0.3,0.3,0.3,0.3), "cm")), axis.text.y = element_text(size = base_size*0.8, color = "black", lineheight = 0.9, margin=unit(c(0.3,0.3,0.3,0.3), "cm")), axis.ticks = element_line(color = "black", size = 0.2), axis.title.x = element_text(size = base_size, color = "black", margin = margin(t = -5)), # t (top), r (right), b (bottom), l (left) axis.title.y = element_text(size = base_size, color = "black", angle = 90, margin = margin(r = -5)), axis.ticks.length = unit(-0.3, "lines"), legend.background = element_rect(color = NA, fill = NA), legend.key = element_rect(color = "black", fill = "white"), legend.key.size = unit(0.5, "lines"), legend.key.height =NULL, legend.key.width = NULL, legend.text = element_text(size = 0.6*base_size, color = "black"), legend.title = element_text(size = 0.6*base_size, face = "bold", hjust = 0, color = "black"), legend.text.align = NULL, legend.title.align = NULL, legend.direction = "vertical", legend.box = NULL, panel.background = element_rect(fill = "white", color = NA), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(size = base_size, color = "black"), ) } data_D1 = read.csv("./data_D1.csv") data_D2 = read.csv("./data_D2.csv") # format data to ggplot's liking data_D = data.frame("width"=c(data_D1$width,data_D2$width), "unit"=c(rep("shear_stress",nrow(data_D1)), rep("velocity",nrow(data_D2))), "value"=c(data_D1$shear_stress,data_D2$velocity) ) head(data_D) curve_D1 = data.frame(width=data_D1$width, shear_stress=33.28/(pi*18*data_D1$width^2)) curve_D1 panel_D1 <- ggplot(data=data_D1, aes(x=width, y=shear_stress)) + geom_point(fill="red", size=3, pch=22) + geom_line(data=curve_D1) + scale_x_log10(expand=c(0,0), # prevent gap between origin and first tick breaks=c(0.5,1,2,5,10,20,50), labels=c(0.5,1,2,5,10,20,50), limits=c(0.5,50)) + scale_y_log10( expand = c(0, 0), # using trans_format from the scales package, but one can also use expressions labels = trans_format('log10', math_format(10^.x)), breaks=c(0.001,0.01,0.1,1), limits = c(0.001,1) ) + annotation_logticks(sides = "l") + theme_void()+ theme( line = element_blank(), # exclude everything outside axes bc it messes with positioning of grob in panel_D text = element_blank(), title = element_blank(), axis.line.y = element_line(colour = "black") ) + ylab("shear stress (Pa)") panel_D1 panel_D <- ggplot(data=data_D2, aes(x=width, y=velocity))+ # add plot of first dataset as grob as a trick to introduce two y-axes with different scalings geom_point(fill="blue", size=3, pch=21) + annotation_custom(ggplotGrob(panel_D1)) + scale_y_continuous(expand = c(0,0), breaks = seq(0,1,0.1), limits = c(0,1), # putting the y axis of the second plot to the right position = "right", # now the secondary axis becomes the left axis # we need the axis text+title for panel_D1 # They were excluded in panel_D1 bc they were messing with the positioning sec.axis = sec_axis(~., name = "shear stress (Pa)", # rescale breaks bc sec_axis inherits scale from primary y axis breaks=rescale(c(-3,-2,-1,0), to = c(0,1)), labels = c(expression("10"^"-3", "10"^"-2", "10"^"-1", "10"^"0"))) ) + scale_x_log10(expand=c(0,0), breaks=c(0.5,1,2,5,10,20,50), labels=c(0.5,1,2,5,10,20,50), limits=c(0.5,50)) + annotation_logticks(sides = "b") + my_theme() + geom_line(color="blue") + geom_vline(xintercept = 1.1, linetype="dashed") + geom_hline(yintercept = 0.9, linetype="dashed") + theme( plot.margin = unit(c(0.1,0,0,0.5), "cm"), # to match other panels axis.title.y = element_text(margin = margin(r=1)), axis.text.y = element_text(margin = margin(r=6)), axis.text.y.right = element_text(margin = margin(l=7)), axis.title.y.right = element_text(angle = 90) ) + xlab(expression(lumen~width~(mu*m))) + ylab("relative flow velocity") + annotate(geom = "text",x =6 ,y =0.85 ,label = "tau == 0.5~Pa",parse=T) + annotate(geom = "text",x =1.4 ,y =0.4 ,label = "b == 1.1*mu*m",parse=T,angle=90) + annotate(geom = "text",x =5.6 ,y =0.6 ,label = "tau == frac(4*mu*Q,pi*a*b^2)",parse=T) panel_D
小编有话说
本文主要学到的知识点如下：
- 使用
annotation_custom(ggplotGrob())图中添加其他图形； - 使用
scale_x_log10()和scale_y_log10()对刻度进行对数变换； - 使用
annotation_logticks(sides = "b")添加 x 轴的 ticks； - 使用
scale_y_continuous(position = "right")改变 Y 轴位置。
