厚缊

诹图——拿破仑远征图

厚缊 / 2019-05-01


一直在用R基础绘图系统作图,今天第一次回到ggplot2的怀抱,感觉十分生疏,很多东西不理解,发自内心感叹,自己知道的东西太少了。这一期的缘起,主要有两方面的原因:一方面,自从看了Murrell开发的vwline,惊叹R还可以这样玩,那时开始就一心想把这个包的主要函数封装成geom_*()函数,方便调用,但反复尝试了很多次,都没有成功;另一方面,Minard的拿破仑远征图一直是经典中的经典,重现经典也是重要的学习策略,当然,前提是我解决了如何封装vwline包的问题。这一期的内容比较多,代码差不多250行,还要说明如何利用grid绘图生态里面的grob函数来封装成ggplot2geom_*()函数,所以得保持耐心。

vwline包主要有6个grob函数,一期是肯定说不完,今天要讲的故事只包括grid.vwline()vwlineGrob()函数,顺带送个应用的拿破仑远征图。

准备工作

最新版本的vwline包依赖于gridBezier,所以要先安装gridBezier包,而且这两个包都没有放在CRAN上,需要通过GitHub进行安装。拿破仑远征图的的数据来源于HistData包,要完成全部例子,也要进行安装。

if(!require(devtools)) install.packages('devtools')
if(!require(ggplot2)) install.packages('ggplot2')
if(!require(dplyr)) install.packages('dplyr')
if(!require(concaveman)) install.packages('concaveman')
if(!require(sf)) install.packages('sf')
if(!require(HistData)) install.packages('HistData')
if(!require(mapdata)) install.packages('mapdata')
if(!require(gridBezier)) devtools::install_github('pmur002/gridbezier')
if(!require(vwline))devtools::install_github('pmur002/vwline/pkg')
library(ggplot2)
library(dplyr)
library(gridBezier)
library(vwline)
library(HistData)
library(mapdata)

自定义ggplot2扩展函数

ggplot2包提供了很好的扩展接口函数原型ggproto,能用最少的代码写geom_*()stat_*()facet_*()等扩展函数。通常来说,若是ggplot2里面已经有类似(父类)的geom_*()函数,只是没有相应的统计变换来达到自己绘图目的,只需要重写stat_*()函数,这种情况比较简单;若是ggplot2中没有类似的,需要自己重写Geom函数,这样的情况会更复杂一些。通过个小例子说说怎么写stat_*()函数,然后把主要精力放在Geom函数部分。

自定义stat_*()函数

在官方示例中,作者讲解了如何写一个stat_chull()函数画包围散点集的最小凸多边形,多边形是由点连线而成的,所以可以理解为这个函数的父类是GeomPoint或者GeomPolygon,因此只需要写统计变换函数就可以了。我这里的例子和这个类似,是画包围散点集的最小多边形,可以是凸的,也可以是非凸的。

transconcave <- function(x, 
                     y, 
                     concavity = 16,
                     length_threshold = 0,
                     buffer = FALSE,
                     dist = NULL,
                     ...){
  if(!(is.numeric(x) && is.numeric(y)))
    stop("'x' or 'y' not a numeric vector.", call. = FALSE)
  if(length(x) != length(y) )
    stop("Length of 'x' or 'y' not equal.", call. = FALSE)
  if(length(x) < 3)
    stop("Length of 'x' or 'y' should greater than 3.", call. = FALSE)
  if(is.null(dist)){
    rngx <- range(x, na.rm = TRUE)
    rngy <- range(y, na.rm = TRUE)
    dist <- min(diff(rngx), diff(rngy)) * 0.1
  }
  m <- matrix(c(x, y), ncol = 2)
  cave <- concaveman::concaveman(m, concavity = concavity, length_threshold = length_threshold)
  x <- cave[ , 1]
  y <- cave[ , 2]
  if(buffer){
    cave_bf <- sf::st_buffer(sf::st_polygon(list(cave)), dist = dist, ...)[[1]]
    x <- cave_bf[ , 1]
    y <- cave_bf[ , 2]
  }
  return(data.frame(x = x, y = y))
}

StatCave <- ggproto('StatCave', Stat, 
     required_aes = c('x', 'y'),
     compute_group = function(data, scales, params, concavity = 1, 
                              length_threshold = 0, buffer = FALSE, dist = NULL) {
                    xy <- transconcave(data$x, data$y, concavity = concavity, 
                                       length_threshold = length_threshold,
                                       buffer = buffer, dist = dist)
                    
                    xy
     }
)

stat_cave <- function(mapping = NULL, data = NULL, geom = 'polygon',
                    position = 'identity', na.rm = FALSE, show.legend = NA, 
                    inherit.aes = TRUE, concavity = 1, length_threshold = 0,
                    buffer = FALSE, dist = NULL, 
                    ...) {
  layer(
    stat = StatCave, data = data, mapping = mapping, geom = geom, 
    position = position, show.legend = show.legend, inherit.aes = inherit.aes,
    params = list(
      concavity = concavity, 
      length_threshold = length_threshold,
      buffer = buffer, 
      dist = dist, 
      na.rm = na.rm, ...)
  )
}

transconcave()函数根据散点集(X坐标和Y坐标)计算围栏多边形的坐标,返回坐标数据框。这个函数利用了concaveman包关于非凸多边形围栏计算方法,参数concavitylength_threshold 传递给concaveman()函数。参数bufferTRUE时,利用st_buffer()函数设置围栏多边形缓冲区,大小由dist来控制。其余参数均传递给st_buffer()函数。

StatCave定义继承于Stat类的统计变换函数。required_aes设置了必须的映射参数,这里是散点集的X坐标和Y坐标;compute_group()函数是对传入的数据参数做统计变换,前三个参数是datascalesparams是必须的,后面的其它参数传递给自定义的统计变换函数transconcave()函数。xy <- transconcave(...)部分对数据进行统计变换,需要解释的是引用映射参数(required_aes内的参数)需要使用data$的形式,因为数据映射之后会被整理成名为data的数据框,required_aes内的参数均独占一列。统计变换完成后返回新的数据框xy

stat_cave()函数的geom参数设置为polygon,表示非凸多边形围栏默认是画多边形,也可以设置成pointpath。从concavity参数开始的所有参数均是自定义统计变换函数参数,根据需要设置默认值。剩余其它参数保持不变。函数体内部是layer()函数,需要注意的两个地方:一是stat换成前面定义统计变换函数StatCave;二是自定义统计变换函数参数通过列表形式赋值给params

data(mpg)
ggplot(mpg, aes(displ, hwy)) + geom_point() + stat_cave(alpha = 0.6)

ggplot(mpg, aes(displ, hwy)) + geom_point() + 
  stat_cave(alpha = 0.6, fill = '#66C2A5', concavity = 5, 
            length_threshold = 1.2, buffer = TRUE, dist = 0.5)

ggplot(mpg, aes(displ, hwy)) + geom_point() + 
  stat_cave(alpha = 0.6, fill = '#66C2A5', concavity = 500, 
            buffer = TRUE, dist = 0.5)

data(mpg)
ggplot(mpg, aes(displ, hwy)) + geom_point() + 
  stat_cave(alpha = 0.6, fill = '#66C2A5', concavity = 5, 
            length_threshold = 1.2, buffer = TRUE, dist = 0.5) + 
  facet_wrap(.~year)

自定义Geom*()函数

有了stat_cave()函数的铺垫,后面的故事容易讲一点,但需要有一点grid包的基础。grid绘图系统的绘图函数有两类:一类是grid.*()函数,直接在当前绘图设备出图,对应于ggplot2中的geom_*()函数;另一类是*Grob()函数,需要调用grid.draw()才能出图,对应于ggplot2中的Geom*()函数。接下来我们先用vwline::vwlineGrob()定义新的Geom*()函数,再把这个函数封装成geom_*()函数。

先看下vwline::vwlineGrob()函数,必要参数是xyw,其它均为可选参数。xy是路径的点坐标,w是路径的宽度,在自定义Geom*()函数时必要参数必须有相应的参数传进来。

args(vwline::vwlineGrob)
## function (x, y, w, default.units = "npc", open = TRUE, linejoin = "round", 
##     lineend = "butt", mitrelimit = 4, stepWidth = FALSE, render = if (open) vwPolygon else vwPath(), 
##     gp = gpar(fill = "black"), name = NULL, debug = FALSE) 
## NULL

default_aes中定义的是默认映射参数,因为vwline::vwlineGrob()函数本质上画的是多边形,相应的映射参数主要为边框颜色(colour)、填充颜色(fill)、边框线粗细(size)、边框线类型(linetype),很多情况下还有透明度(alpha),但我觉得用透明度区分不同多边形难度比较大,就没有放在默认映射参数里面。最重要的映射参数是宽度w,默认是画宽度为1cm的线。

draw_group定义了绘图函数,前三个参数datapanel_paramscoord不用管,直接写上就行,后面的参数(除了na.rm)都是传递给vwline::vwlineGrob()的,对着原函数的默认参数一个个写上(这里多了两个参数widspunits,后面会解释)。

接下来直接跳到trans那一行,这一行很重要,是对你的数据做坐标变换,例如极坐标变换、对数变换等,变换参数都不用特殊设定。紧接着一行取了变换之后数据(若是有分面、分组就是多个数据框组成的列表)的第一行,这是因为每个多边形的边框颜色、填充色、边框线宽、线型类型均只要设置一次就可以,而不是每个点设置一次。

ggplot2:::ggname(...)是封装vwline::vwlineGrob()的部分,映射参数的参数值要从变换后的数据(trans)中取,其它跟着原函数参数一个个写下来。最后一行代码draw_key = draw_key_polygon是处理图例,vwline::vwlineGrob()是画多边形,所以选了多边形图例,ggplot2内置了多种图例,详细情况参见原文档。

最后说说多的两个参数widspunitswidsp本质上也是和w一样的宽度参数,但是我在处理的过程中发现,vwline::vwlineGrob()支持的widthSpec类宽度返回值是列表,若直接传进GeomVwlineGrob中,映射函数会把vwline::widthSpec()返回列表当成向量处理,导致数据映射出错,所以退而取其次,就新加了这个参数,当设置widthSpec类宽度时,通过widsp在映射函数aes()外设置。units参数是为了更灵活的设置线型宽度单位,并且和vwline::widthSpec()函数中的default.units参数区别开。

GeomVwlineGrob <- ggproto('GeomVwlineGrob', Geom,
    required_aes = c('x', 'y'),
    default_aes = aes(w = unit(1, 'cm'), colour = 'black', fill = 'grey40', 
                      size = 0.5, linetype = 1),
    draw_group = function(data, panel_params, coord,
                          open = TRUE,
                          linejoin = 'round',
                          lineend = 'butt', 
                          mitrelimit = 4,
                          stepWidth = FALSE,
                          render = if (open) vwPolygon else vwPath(),
                          name = NULL, 
                          widsp = NULL,
                          debug = FALSE, 
                          units = 'cm', 
                          alpha = NA, 
                          na.rm = FALSE) {
        if(is.null(widsp)){
           if(!is.unit(data$w)) {
              data$w <- unit(data$w, units)
           }
        }
        if(!is.null(widsp)){
           if(class(widsp) != 'widthSpec')
           stop("'widsp' not a 'widthSpec' class.", call. = FALSE)
        }
        n <- nrow(data)
        if (n < 2) return(zeroGrob())
        trans <- coord$transform(data, panel_params)
        first_row <- trans[1, , drop = FALSE]
                                   
        ggplot2:::ggname('geom_vwline',
        vwlineGrob(
        x = trans$x, 
        y = trans$y, 
        w = if(is.null(widsp)){
               trans$w
            }else {
               widsp
            }, 
        open=open, 
        default.units = 'native',
        linejoin=linejoin,
        lineend=lineend,
        mitrelimit=mitrelimit,
        stepWidth=stepWidth,
        render=render,
        name=name, 
        debug=debug, 
        gp = gpar(
                  col = scales::alpha(first_row$colour, alpha),
                  fill = scales::alpha(first_row$fill, alpha),
                  lwd = trans$size * .pt,
                  lty = trans$linetype
                  )))
    },
    draw_key = draw_key_polygon
)

有了前面自定义stat_cave()函数的基础,写geom_*()函数就比较简单了。一般情况下,mapping参数到inherit.aes直接写上,接下来的参数都是传递给先前定义的GeomVwlineGrob()函数,所以要对照进行设置。(注意:所有映射参数不需要在geom_*()函数中设置)。最后是layer函数,除了geom要设置为前面自定义的GeomVwlineGrob()外,其它对照着一条一条的写。

geom_vwline <- function(mapping = NULL, data = NULL, stat = 'identity', 
    position = 'identity', ..., show.legend = NA, inherit.aes = TRUE,
    open = TRUE, 
    linejoin = 'round',
    lineend = 'butt', 
    mitrelimit = 4,
    stepWidth = FALSE,
    render = if (open) vwPolygon else vwPath(),
    name = NULL,
    widsp = NULL,
    debug = FALSE, 
    units = 'cm',
    na.rm = FALSE) {
  layer(
    data = data,
    mapping = mapping,
    stat = stat,
    geom = GeomVwlineGrob,
    position = position,
    show.legend = show.legend,
    inherit.aes = inherit.aes,
    params = list(
      open=open,
      linejoin=linejoin,
      lineend=lineend, 
      mitrelimit=mitrelimit,
      stepWidth = stepWidth,
      render=render,
      name=name, 
      widsp = widsp,
      debug=debug,
      units = units,
      na.rm = na.rm,
      ...
    )
  )
}

画几个好玩的图测试下效果。在数据框外手动设置宽度的单位时,w度向量长度和数据框的行数一致。

df1 <- data.frame(x = c(.2, .2, .8, .8),
                  y = c(.2, .8, .8, .2))
w = unit(1:4/4, 'in')

ggplot(df1, aes(x, y)) + geom_vwline(w = w, open = FALSE) +
  geom_vwline(w = 0.001, color = 'white', fill = 'white', open = FALSE) +
  coord_cartesian(xlim = c(0, 1), ylim = c(0, 1))

d <- density(ToothGrowth$len)
x <- (d$x - min(d$x)) / diff(range(d$x))
w <- 0.2 * d$y / max(d$y)
df2 <- data.frame(x = x, y = 0.5, w = w)
ggplot(df2, aes(x, y)) + geom_vwline(w = w, units = 'native', fill = '#66C2A5')

ggplot(df2, aes(x, y)) + 
  geom_vwline(w = w, units = 'native', fill = '#66C2A5', open = FALSE) +
  coord_polar()

拿破仑远征图

为了出这幅图,做的准备工作有点多了,几乎把ggplot2扩展全部说了一遍。说实话,ggplot2我也不是很熟悉,在写这篇文章的过程中没有少犯低级错误,都是通过查帮助一点一点的去克服,好在作为R最流行的绘图软件包,ggplot2生态系统足够完善,各种资料应有尽有,遇到问题不妨自己好好查查,总能找到解决方案。

数据准备

拿破仑远征图的数据收录在HistData包中,分别是Minard.troopsMinard.citiesMinard.temp,对应行军路径、沿途主要城市位置和期间气温变化数据。这里我们只用了前两个。

首先生成新变量w,用每个节点剩余人数和出发时总人数的比例表示,这也是Minard定义远征图中路径宽度的变量。接下来按照行军方向(direction)和组(group)将Minard.troops数据集分成6组,其中A表示进攻,R表示撤退。之后的部分主要是我觉得这个数据画出来的图有点丑,对数据做了一些微调,不是必要的。最后两行代码是提取mapdata中的河流数据,尽量把行军范围内的主要河流还原。

minard <- Minard.troops %>% 
  mutate(w = survivors / max(survivors))
directionA1 <- filter(minard, direction == 'A', group == 1)
directionA2 <- filter(minard, direction == 'A', group == 2)
directionA3 <- filter(minard, direction == 'A', group == 3)
directionR1 <- filter(minard, direction == 'R', group == 1)
directionR2 <- filter(minard, direction == 'R', group == 2)
directionR3 <- filter(minard, direction == 'R', group == 3)
directionR1[5, 1:2] <- c(35.33, 55.23)
directionR1[6, 2] <- 55.06
directionR1[7, 2] <- 54.74
directionR1[8, 2] <- 54.62
directionR1[19, 1] <- 23.7
directionR2[2, 2] <- 54.27
directionR2[3, 2] <- 54.17
directionR2[4, 1:2] <- c(28.32, 54.23)
directionR3[1:3, 1] <- c(24.5, 24.1, 24)
directionA1[14, 1:2] <- c(35.52, 55.38)
directionA1[1, 1] <- 23.5
directionA3[1:2, 2] <- c(55.15, 55.25)

rivers <- map("rivers",plot=F,add=T)
rivers <- data.frame(long = rivers$x, lat = rivers$y)

主图

主图总共叠加了8个图层,分别是撤退路线3个、向莫斯科进攻路线3个、河流1个和城市标签1个。有几个参数需要解释下:w表示行军路径宽度,映射给剩余士兵比例,单位是英寸,也可以在aes()中指定,然后设置默认单位units = 'in'geom_vwline()函数默认stepWidth参数为FALSE,即不同宽度之间平滑递减,这里不需要平滑,所以设置为TRUElinejoin参数设置线的连接方式,按照Murrell的建议,设为mitre。最后通过xlim()ylim()限定坐标范围,目的是利用这两个函数会自动丢弃超出范围的绘图数据的机制来截取河流。

注意:我把主图赋值给了mainfig,是方便后续处理,要直接看这部分的效果可以去掉,也可以用grid.draw(mainfig)查看。

mainfig <- ggplot() + 
  geom_vwline(data = directionR1, aes(x = long, y = lat), 
              w = unit(directionR1$w, 'in'),
              stepWidth=TRUE, linejoin='mitre', 
              color = 'grey50', fill = 'grey50') +
  geom_vwline(data = directionR2, aes(x = long, y = lat), 
              w = unit(directionR2$w, 'in'),
              stepWidth=TRUE, linejoin='mitre', 
              color = 'grey50', fill = 'grey50') +
  geom_vwline(data = directionR3, aes(x = long, y = lat), 
              w = unit(directionR3$w, 'in'),
              stepWidth=TRUE, linejoin='mitre', 
              color = 'grey50', fill = 'grey50') +
  geom_vwline(data = directionA1, aes(x = long, y = lat),
              w = unit(directionA1$w, 'in'),
              stepWidth=TRUE, linejoin='mitre', 
               color = '#E5CBAA', fill = '#E5CBAA') + 
  geom_vwline(data = directionA2, aes(x = long, y = lat),
              w = unit(directionA2$w, 'in'),
              stepWidth=TRUE, linejoin='mitre', 
              color = '#E5CBAA', fill = '#E5CBAA') + 
  geom_vwline(data = directionA3, aes(x = long, y = lat),
              w = unit(directionA3$w, 'in'),
              stepWidth=TRUE, linejoin='mitre', 
              color = '#E5CBAA', fill = '#E5CBAA') +
  geom_path(data = rivers,
            aes(x = long, y = lat), color = 'grey70') +
  geom_text(data = Minard.cities, 
            aes(x = long, y = lat, label = city),
            family = 'Felipa-Regular', size = 5
            ) +
  xlim(c(23.5, 38)) + 
  ylim(c(53.7, 55.9)) +
  theme_void()

其它辅助元素

这部分添加注释文本和标题。这里写的有点复杂,思路是借助grid包,定义三个视图,分别画主图、注释文本和标题,这一块我也有点晕头转向,说不太清楚。更简洁的思路是利用R中的拼图工具,直接把上下两部分叠加在一起。

note <- "Les nombres d’hommes présents sont représenté s par les largeurs des zones colorées a raison d’un millimètre pour dix mille hommes; ils sont de\n
plus écrits en travers des zones. Le gris désigne les hommes qui entrent en Russie, le noir ceux qui en sortent. Les renseignements qui ont servi à \n 
dresser la carte ont été puisés dans les ouvrages de MM. Thiers, de Ségur, de Fezensac, de Chambray et le journal inédit de Jacob, pharmacien \n
de l’ar– mée depuis le 28 octobre. Pour mieux faire juger à l'ail diminution de lal’armée, j’ai scha et Witebsk, avaient toujours marsupposé que\n
les corps du prince Jérôme et du Maréchal Davoust qui avaient été détachés sur avec l’armée Minsk et Mobilow et ont rejoint vers Orché.\n"

grid.newpage() 
note_width <- 
pushViewport(
  viewport(
    x = unit(0, 'npc'), y = unit(0, 'npc'),
    width = unit(1, 'npc'), height = unit(3 / 5, 'npc'),
    just = c('left', 'bottom'),
    name = 'vp1'
  )
)
print(mainfig, newpage = FALSE)
popViewport()
pushViewport(
  viewport(
    x = unit(0.1, 'npc'), y = unit(3 / 5, 'npc'),
    width = unit(0.85, 'npc'), height = unit(2 / 5, 'npc'),
    just = c('left', 'bottom'),
    name = 'vp2',
    gp = gpar(fill = '#00000000')
  )
)
grid.text(label = note, 
          x = unit(0, 'npc'),
          y = unit(1, 'npc') - unit(1.6, 'inches'),
          just = c('right', 'top'),
          hjust = 0, 
          vjust = 1,
          gp = gpar(fontfamily = 'Felipa-Regular',
                    cex = 1.3,
                    col = 'grey20',
                    lineheight = 0.75
                    )
          )
popViewport()
pushViewport(
  viewport(
    x = unit(0.15, 'npc'), y = unit(1, 'npc') - unit(1.5, 'inches'),
    width = unit(0.7, 'npc'), height = unit(1, 'inches'),
    just = c('left', 'bottom'),
    name = 'vp3',
    gp = gpar(fill = '#00000000')
  )
)
grid.text(label = 'Carte figurative', 
          x = unit(0, 'npc') + unit(2, 'mm'),
          y = unit(0.9, 'inches'),
          just = c('left', 'top'),
          hjust = 0, 
          vjust = 1,
          gp = gpar(fontfamily = 'DrSugiyama',
                    fontface = 2,
                    cex = 3,
                    col = 'grey40',
                    lineheight = 0.75
          )
)

grid.text(label = 'des pertes', 
          x = unit(1.1, 'npc') - unit(2, 'mm'),
          y = unit(0.6, 'inches'),
          just = c('left', 'top'),
          hjust = 1, 
          vjust = 1,
          gp = gpar(fontfamily = 'DrSugiyama',
                    fontface = 2,
                    cex = 2,
                    col = 'grey40',
                    lineheight = 0.75
          )
)

popViewport()