带南海九段线分位数地图可视化（R语言版）

今天带来一篇承诺虾神的可视化博客。内容是使用R语言进行带南海九段线分位数地图可视化。虾神的原博文地址如下（Python版）。

Python实现带南海九段线分位数地图完整可视化版本（附代码及数据）

1999-2017年中国各省旅游外汇收入分析及可视化（附代码及数据）

数据及代码github地址

1 数据下载

虾神把代码和数据放在了github上，没接触过github的人可能对如何下载不太熟悉，这里也简单介绍下两种方式。如果想进一步了解git这种神器的可以安装git，然后按第一种方式下载（以下介绍默认已安装git bash）。而第二种方式则完全不用了解git方面的内容。

1 git clone

首先新建一个你放数据和分析的文件夹。然后先点击浏览虾神提供的github地址，在页面中点击clone or download，即跳出如下的页面，复制https的地址。

然后右击打开git bash。敲入命令行。

git clone https://github.com/allenlu2008/PythonDemo.git

然后敲击回车，即可发现开始下载。当然由于虾神这个仓库内容比较多，是一个比较漫长的下载过程。

2 直接下载（Download ZIP）

第二种方式依旧是点击clone or download。这回是点击Download ZIP。

接下来就进入传统的浏览器点击下载存储文件的范畴了。后续就不阐述了。

2 R语言可视化

1 所需包的安装

由于R语言需要一些包的支持才能实现读取空间数据和可视化。具体这块前面有博客介绍过。这里不详述。

R语言读取空间数据以及ArcGIS中OLS工具回归结果可视化R语言版

这里主要提一下需要本次可视化需要安装的几个包：rgdal，ClassInt以及sp包。

2 绘图

说到绘图的话，其实R语言中有多套绘图系统（我称之为系统是指他们的语法各有特色，使用起来有所差异），比如最基础的plot系统，lattice系统，闻名遐迩的ggplot2系统。本篇博客最早是打算使用基于lattice的spplot绘图来进行。后面发现利用spplot绘制这个图并不简单，但我也会先介绍使用spplot的画法。最后我采用的是plot系统来绘制。本次绘图函数样例选择2008年外汇收入数据来进行。此外后续出现的whsr表示外汇收入表格，chinau表示使用的中国省份地图数据（关联有外汇收入数据），jdx表示九段线地图。

由于本次是绘制分位数地图，因此在读取数据的基础上需要先计算分位数。从虾神Python代码里也发现有计算分位数进行分类的代码。这里也在R里面写了一下。这里有两种方式，一种是采用R包ClassInt来做分类。就如代码的前两行，一行代码解决。此外也可以使用这个包做熟悉的自然间断点法分类。不过似乎这个分类与虾神的分位数地图结果有点不一致，所以我也写了一个自己手动计算的代码。quantinle是R的分位数函数，输入参数是给定一组向量（数据）以及需要的分位（例如25分位数就是0.25），返回值对应分位数的值。具体的计算步骤见代码。不详述。

##use ClassInt
q5 <- classIntervals(chinau$y2008, 5, style = 'quantile')

##self calculation
q25 <- quantile(whsr$y2008, probs = c(0.25))
q50 <- quantile(whsr$y2008, probs = c(0.5))
q75 <- quantile(whsr$y2008, probs = c(0.75))
iqr <- q75 - q25
u <- q75 + 1.5 * iqr
d <- q25 - 1.5 * iqr
u <- ifelse(u > max(whsr$y2008), max(whsr$y2008), u)
d <- ifelse(d < min(whsr$y2008), min(whsr$y2008), d)
box <- c(d, q25, q50, q75, u)

获取到box分位数分界线后，在R里可以通过ifelse给不同省份赋值（或者叫分类，就是按照分位数范围将外汇收入划分为不同类别），这里将分类设置为plot字段（值为1到6）。代码如下。

chinau$plot <- ifelse(chinau$y2008 < box[1], 1, 
                      ifelse(chinau$y2008 < box[2], 2, 
                             ifelse(chinau$y2008 < box[3], 3, 
                                    ifelse(chinau$y2008 < box[4], 4, 
                                           ifelse(chinau$y2008 < box[5], 5, 6)))))

接着可以设置色带。这个颜色设置看个人喜好，我是随机挑了一种配色方案，使用colorRamPalette。

colspic <- colorRampPalette(c("#f38181", "#fce38a", "#eaffd0", "#95e1d3"))(length(box)+1)

以上就是主要的准备工作，接下来就进行可视化语句就行了。不过这是基于spplot系统的绘制地图准备工作，如果是基于plot还会多一步，后面会说。这里先关注如何用spplot来绘制地图。先放代码和结果再来解释。

spplot(chinau, zcol = "plot", scales = list(draw = T), 
       main = "2008年中国各省国际旅游外汇收入分位数专题图\n数据来源:国家统计局", col.regions = colspic,
       at = seq(1, 7, 1))

等一等，为什么有X。

其实我后面画的时候发现这个系统画这个图感觉乱七八糟很麻烦，连画九段线也有很多毛病，所以就只是先展示下spplot怎么绘图，后面有空再来水一篇博——啊呸，说错了，再来写一篇博客。前面那个结果图是有问题的，仅仅是教学示意，一定要加上九段线，因为我们中国一点都不能少！！！

先来讲spplot的语法，首先第一个参数是绘图的矢量或栅格数据（这里是chinau），zcol表示要用来映射颜色的字段（选用了之前的plot字段，即参照分位数的分类），scales = list(draw = T)意思是要把经纬度标出来，如果是F就不标。main是填标题，at是映射字段的绘制数值。其他参数详情参见文档或者等以后我来水——呸，我来介绍。

这里也展示下ClassInt分类的结果。当然也是教学用，没有画九段线。

接下来进入plot系统的正经实现效果。首先就如虾神博客效果，我们知道本次实现的可视化本身就是一个拼图操作，所以先需要用par函数来分割绘图位置。par函数的参数十分丰富。这次主要用两个，一个是fig，一个是new。代码如下。

par(fig = c(0, 0.3, 0, 1), new = T)

fig的参数是一个向量，形式是(x1,x2,y1,y2)，这个向量是以画布左下角为原点，其实表示的就是接下来绘制的图片占据图片位置和大小。按照这里的代码，接下来绘制的图是占画布横轴的前面30%，纵轴的全部。new是表示允许叠加新的图，这样才能把两个图画在一个画布上。解释得有点绕口，但是如果你自己尝试画个图就能理解，所以我就不继续展开了。

接下来先将plot绘制地图的关键做一些简要介绍，这里需要先在plot函数可视化前，做一个准备工作。即新建一个字段col，把颜色映射到这个字段上（字段值是颜色的十六进制），也是用ifelse的方法。做完之后就可以开始做可视化。

由于虾神给的九段线数据是只有九段线那块区域，数据经纬度范围与省份数据不相同，所以可视化的第一步首先要计算一个能包含这两个数据的数据经纬度范围。这里就需要bbox了。

你说的是这个bbox？

还是这个bbox？

或者是这个？

哦，对不起，串戏了。是这个。

这个函数可以把空间数据的坐标范围返回给我们。这样子就可以通过简单的计算得到包含两个数据的经纬度范围。当然我把所有的代码都连在了一行，所以看起来比较复杂。这个函数首先先绘制九段线，同时，绘制的数据范围是包含两个数据的。xlim和ylim是表示绘制的数据范围，需要的是一个二元向量（如c(20,30），即表示绘制经（纬）度的范围是20到30）。main跟spplot一样是标题。axes表示不需要画任何坐标轴。这样子结果如图。

plot(jdx, 
     xlim = c(min(c(bbox(chinau)[1,1], bbox(jdx)[1,1])), max(c(bbox(chinau)[1,2], bbox(jdx)[1,2]))),
     ylim = c(min(c(bbox(chinau)[2,1], bbox(jdx)[2,1])), max(c(bbox(chinau)[2,2], bbox(jdx)[2,2]))), 
     main = '2008年中国各省国际旅游外汇收入分位数专题图\n数据来源:国家统计局',
     axes = F)

现在先把省份的数据绘制上去就能得到基本的地图。col是颜色映射的字段，add表示可以把新的图绘制到原来的图上，相当于PS的叠加一个新图层。这样我们就看到了基本的图形。

plot(chinau, col = chinau$col, add = T)

上面的图很单调，需要格网做点缀。这个时候就需要另外一个函数axis了。axis的第一个参数是side，取值为1（下）,2（左）,3（上）,4（右），意思就是画上下左右的四条坐标轴。at表示要绘制的刻度值范围，tck是表示刻度值的长度，0的话是不标刻度，1的话是网格线，默认为-0.01，col是颜色。这里再次用到了bbox，这是为了划分网格，这里采用的是10经度为划分标准。当然我还加上了四条轴，具体的代码后面会分享。这里就不详述了，主要阐述函数用法。结果如图，越来越接近虾神效果图了。

axis(1, at = seq(round(min(c(bbox(chinau)[1,1], bbox(jdx)[1,1]))-10), 
                 round(max(c(bbox(chinau)[1,2], bbox(jdx)[1,2]))+10),
                 10), 
     tck = 1, col = 'grey')

地图目前相较于虾神的图（标注不考虑），最关键的差距就是图例了。这里使用polygon和text来绘制图例。polygon是绘制多边形，首先需要给出组成x和y的多边形坐标，border是线的颜色，col是多边形填充颜色。text是根据坐标标注文字。x,y是点坐标，字符串是标注内容。当然生成多边形坐标以及标注文字坐标依旧使用了bbox，这里可以根据具体地图替换主要改变前面的系数:1.1还有后面的加数灵活绘制。效果如图。

polygon(x = c(rep(round(1.1*min(c(bbox(chinau)[1,1], bbox(jdx)[1,1]))), 7), 
              seq(round(1.1*min(c(bbox(chinau)[1,1], bbox(jdx)[1,1]))), 
                  round(1.1*min(c(bbox(chinau)[1,1], bbox(jdx)[1,1])))+6, 
                  1),
              rep(round(1.1*min(c(bbox(chinau)[1,1], bbox(jdx)[1,1])))+6, 7),
              seq(round(1.1*min(c(bbox(chinau)[1,1], bbox(jdx)[1,1])))+5, 
                  round(1.1*min(c(bbox(chinau)[1,1], bbox(jdx)[1,1]))), 
                  -1)),
        y = c(seq(round(1.1*min(c(bbox(chinau)[2,1], bbox(jdx)[2,1]))), 
                  round(1.1*min(c(bbox(chinau)[2,1], bbox(jdx)[2,1])))+3, 
                  0.5), 
              rep(round(1.1*min(c(bbox(chinau)[2,1], bbox(jdx)[2,1])))+3, 7),
              seq(round(1.1*min(c(bbox(chinau)[2,1], bbox(jdx)[2,1])))+2.5, 
                  round(1.1*min(c(bbox(chinau)[2,1], bbox(jdx)[2,1]))), 
                  -0.5),
              rep(round(1.1*min(c(bbox(chinau)[2,1], bbox(jdx)[2,1]))), 7)), 
        border = 'black',
        col = 'white')
text(x = round(1.1*min(c(bbox(chinau)[1,1], bbox(jdx)[1,1])))+10, 
     y = (round(1.1*min(c(bbox(chinau)[2,1], bbox(jdx)[2,1]))) + 
                  round(1.1*min(c(bbox(chinau)[2,1], bbox(jdx)[2,1]))) + 3)/2, "无数据")

后续就一步步添加图例，再加上箱线图就可以得到虾神的效果图了。箱线图绘制是R语言中比较简单的，使用boxplot即可。

当然我们也可以做一版有小窗的，原理同上面一样，也是设置par，把九段线部分重绘制一遍。结果如图。

最后虾神要求的，箱线图+点绘图。这里也描述下。由于只针对单年数据所有省份箱线图，先设置一个类别字段让所有省份数据都相同，设置为1，然后以这个字段来绘制boxplot箱线图。加上点的思路就是以设置的类别字段作为x，外汇收入为y，颜色与地图相同。代码如下。详细意义就不阐述了。详情可参见boxplot的文档。效果最后如图。

chinau$boxu <- 1

boxplot(y2008 ~ boxu, data = as.data.frame(chinau))
points(factor(chinau$boxu), chinau$y2008, col = chinau$col, pch = 16)

打完收工，至于动图，由于R语言做动图会牵扯到另外一个软件，我们就等下一篇博客再水——呸，再介绍。所有的代码后面会公开给大家。请稍安勿躁。本文所有使用数据来自于虾神Github，解释权由他说了算。最后一句——bbox牛逼。