带南海九段线分位数地图可视化（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，把九段线部分重绘制一遍。结果如图。

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