R语言：如何在一张图上画多个填色等值线图

　　在R语言：填色等值线图及其色标（color bar）设置中我们介绍了filled.contour函数的用法，它可以很方便的绘制带色标的填色等值线图。但是我们平时可能更多的需要将多个填色图放在同一张图上（如下图所示），这种图该如何绘制？

　　相比R语言：填色等值线图及其色标（color bar）设置介绍的图形，这里的难点主要有两个，一是怎样进行分面，二是怎样在地图上绘制陆地轮廓，下面我们将依次解决这两个难题。

一、一页多图

　　绘图时进行分面操作常用的函数有par函数和layout函数，其中使用par函数分面主要通过调整mfrow参数实现，例如par(mfrow= c(2, 3))即是把当前绘图区域等分为2行3列，其缺点是只能对页面进行等份，而我们绘图的页面布局如下，显然每个绘制区域并不相同，par(mfrow)方法并不适用。

　　layout函数则可以处理这类问题，简单介绍一下layout函数：

　　layout(mat,
　　　　　 widths =rep.int(1, ncol(mat)),
      　　　 heights = rep.int(1, nrow(mat)),
　　　　　 respect = FALSE)

　　其中参数mat是一个由整数1-N构成的矩阵，这里假定我们要在同一个页面上绘制N个图形，数字1-N即是这N个图形绘制的顺序，而数字1-N在矩阵中的相对位置则表示图像在页面上出现的相对位置。在这里我们先绘制填色图，然后绘制色标，最后绘制标题，页面布局矩阵如下：

　　参数widths和heights是两个向量分别表示矩阵中每一列的宽度和每一行的高度，这里如果widths和heights给出的值是数字的话，则表示相对比例，如果是lcm函数的话则是以cm计量的绝对距离。在这里我们的页面布局中，宽度是等分的，因此widths可以设置为c(12, 12, 12)；在高度上标题和色标的高度要比主体的填色图小一点，所以heights设置为c(2, 4, 4, 4, 4, 2)。

　　参数respect如果为TRUE，表示绘图的宽高比严格按照上面widths和heights给定的比例。

mat <-rbind(14, matrix(1:12, 4, 3), 13)
widths <-c(12, 12, 12)
heights <-c(2, 4, 4, 4, 4, 2)
layout(mat,widths = widths, heights = heights)
nf <-layout(mat, widths = widths, heights = heights)
layout.show(nf)

　　上面的代码中，我们使用layout.show函数将页面布局画出（如下图所示），与我们之前设计的一样。

　　但是，即使是使用layout函数进行页面划分之后，使用filled.contour函数绘制填色等值线图，依旧不会出现分面的效果，

mat<-rbind(14, matrix(1:12, 4, 3), 13)
widths <-c(12, 12, 12)
heights <-c(2, 4, 4, 4, 4, 2)
layout(mat,widths = widths, heights = heights)
filled.contour(x= x, y = y, z = z, levels = levels,
               las = 1,
               plot.title = title(main ="Jan", cex.main = 2),
               plot.axes = list(axis(1,seq(100, 160, by = 20), c('100E', NA, '140E', NA)),
                                axis(1, 180, 180),
                                axis(1,seq(200, 280, by = 20), c(NA, '140W', NA, '100W', NA)),
                                axis(2,seq(-20, 20, by = 10), c('20S', '10S', 'Eq', '10N', '20N'))),
               key.title = title(main ='degC'))

　　这是因为filled.contour函数自身已经调用过一次layout函数，我们直接在R中输入filled.contour，可以看到该函数的代码，原来filled.contour绘制的填色图和色标的分面已经使用过layout函数。但是不要着急，我们继续阅读代码，可以发现其实filled.contour函数中绘制填色图的部分使用了一个内置函数.filled.contour(x, y, z, levels, col)，它只有4个参数，用法与filled.contour函数相同，我们可以使用.filled.contour绘制我们需要的图形。

# 分面
mat<-rbind(14, matrix(1:12, 4, 3), 13)
widths <-c(12, 12, 12)
heights <-c(2, 4, 4, 4, 4, 2)
layout(mat,widths = widths, heights = heights)
# 页面设置
mar <- c(2.5,4, 2, 1)
par(mar = mar,las = 1)
# 经纬度标签
lon_num <-seq(120, 270, 30)
lat_num <-seq(-10, 10, 10)
lon_lab <-c("120E", "150E", "180",  "150W", "120W","90W")
lat_lab <-c("10S", "EQ", "10N")
# 绘填色图
plot( xlim,ylim,
      xlim = range(xlim), ylim = range(ylim),
      xaxs = 'i', yaxs = 'i',
     xaxt = 'n', yaxt = 'n', type = 'n', ann=F)
 
.filled.contour(x,y, z, levels = levels, col = col)
# 绘陆地轮廓
map('world2Hires',xlim = range(x), ylim = range(y), add = T)
# 补齐边框
axis(1, x,labels = NA, tcl = 0)
axis(3, x,labels = NA, tcl = 0)
axis(2, y,labels = NA, tcl = 0)
axis(4, y,labels = NA, tcl = 0)
# 经纬度标签
axis(1, lon_num,lon_lab,
     mgp = c(4, 1.2, 0),
     cex.axis = 1.3)
axis(2, lat_num,lat_lab,
     mgp =c(4, 0.8, 0),
     cex.axis = 1.3)
# 月份标签
axis(3, wz ,labels = 'Jan',
     mgp = c(3, 0.4, 0),
     tcl = 0,
     cex.axis = 1.6)

　　成功啦！

　　之后可以写一个循环，将12个月的数据都绘制出来，最后绘制色标和标题，代码如下

# 色标
mar <- c(3,4, 1.5, 6)
unit <-expression(paste(''^'o','C', sep = ''))
nlev <-length(levels)
par(mar = mar,las = 1)
plot(x = 1:nlev,y = rep(1, nlev),
     xaxs = 'i', yaxs = 'i',
     xaxt = 'n', yaxt = 'n', type = 'n', ann =F)
rect(xleft =1:(nlev - 1),
     ybottom = rep(0, (nlev - 1)),
     xright = 2:nlev,
     ytop = rep(2, (nlev - 1)),
     col = col, border = 1)
par(new=T)
plot(x = 1:nlev,y = rep(1, nlev),
     xaxs = 'i', yaxs = 'i',
     xaxt = 'n', yaxt = 'n', type = 'n', ann =F)
axis(1, at =(2:(nlev - 1)),
     label = round(levels[2:(nlev - 1)],2),
     tcl = 0, cex.axis = 1.8)
axis(4, at = 1,
     label = unit,
     tcl = 0, cex.axis = 2)
 
# 标题
title <- 'TheDevelopment of SSTA in 1997'
size_t <- 3
par(mar = c(0,0, 0, 0),las = 1)
plot(c(0, 4),c(0, 4),
     xlim = c(0, 4), ylim = c(0, 4),
     xaxs = 'i', yaxs = 'i',
     xaxt = 'n', yaxt = 'n', type = 'n', ann =F, axes = F)
text(x = 2, y =2, labels = title, cex = 3)

二、绘制陆地轮廓

　　陆地和国家的轮廓信息我们可以调用maps包或者mapdata包，其中mapdata里的地图信息更加丰富，我们这里使用mapdata包的world2Hires数据集。首先我们将要使用的地图信息绘制出来

x <- seq(from= 98.5, to = 288.5, by = 1)
y <- seq(from= -20.5, to = 20.5, by = 1)
map('world2Hires',xlim = range(x), ylim = range(y))

　　map函数中的add参数如果设置为TRUE，则会在现有的图像上叠加地图，例如：

plot( range(x),range(y),
      xlim = range(x), ylim = range(y),
      xaxs = 'i', yaxs = 'i',
      xaxt = 'n', yaxt = 'n', type = 'n',ann=F)
map('world2Hires',xlim = range(x), ylim = range(y), add = T)

拓展阅读：

数据来源可参考R语言处理气象数据：NetCDF格式数据的读写

filled.contour函数的用法可参考R语言：填色等值线图及其色标（color bar）设置