R语言—热力地图复合气泡饼图

之前有整理R语言地图的绘制思路，主要包括三种方式，但是在实际应用中，因为不同的专业一般只会用到其中的一种，所以我最常用到的就是从第三方下载地理的空间多边形数据，然后使用ggplot2包做出地图即可。不过上篇写得相对理论些，我更倾向于从实践出发，以目标为导向，以实践学习为目的。
这次的图形时带有扇形图的热力地图，描述起来就是在热力地图的特定位置添加上扇形图。扇形图的半径表示数据的大小，扇形图中的每个扇叶表示每个子集的比例。很适合用于表示在不同地理空间上的分类变量的占比。他大概就长这样：

要实现这个目的，刚开始有几种想法，因为热力地图非常容易，问题就在扇形图怎么搞：一是用个for循环，生成要展示的一个个扇形图，然后根据坐标位置将这些个扇形图分别放到合适的位置；二是直接生成一个个的扇形图元素，然后手动调到合适的位置。刚开始我就是这两种想法，但是实现起来的时候错误百出。。。不是说这些想法没法实现，事实上，在我后来查找相关方法的时候也发现确实有大神这么做的，也是能做出来的，但是这些个做法严重的违背了python之禅：能用一行代码实现的，绝对不用两行。
既然我能想到用for循环去实现这种想法，那么肯定会有大神为我们造好了轮子，我就花点时间找找看，如果没有的话，嘿嘿，那我造出来轮子，我就成大神了。。。但不幸的是，我的大神之路夭折了。在ggplot2中，有geom_scatterpie()函数，就是用来实现这种图形的，那就直接拿过来用了。下面用的所有数据都可以直接下载，就算是用作训练的吧，知道数据的形式，再用自己的数据来做就方便了。

一、做热图

第一步：拿到地图数据

根据自己的需要，找到所描述地区的地理信息数据。我要画的是中国地图，那就找到中国地图的数据就行，地理信息数据有很多个，我就找到三个文件，并且我也不知道这三个文件里面到底是什么，毕竟不是专业的，没精力搞这些，反正多一个无所谓，少一个就不行，三个文件的后缀分别是.dbf、.shp、.shx（点击这里下载，提取码：lat2）。这三文件都要放在同一位置。

第二步：R语言读取地理信息数据

就是读取那个.shp后缀的文件，需要以这个包：rgdal

if (!requireNamespace('rgdal',quietly = T)){
    install.packages('rgdal')

setwd('自己工作空间的目录')

china_map=rgdal::readOGR("bou2_4p.shp")

这个时候地理信息数据就已经导入到china_map中去了，这个数据框的类型是空间多边形数据框，里面有好多个变量和维度，具体什么意思，我还没搞懂，但是不影响使用，因为后面用到的，只有其中的几个变量而已，我们把这几个变量搞懂就足够了。

第三步：画地图

地理信息数据存在了china_map中，可以直接用这个作为数据框来画图，但是之前也说了，这个东西他是空间多边形数据框，里面的信息实在是不太好理解，所以就把他转化为一般的数据框，二维的，变量加观测的那种。

china_map_frame <- fortify(china_map)

这个函数就可以了，同样，我也还没搞懂这个函数是怎么样把那个空间多边形数据变成数据框的，但是我的目的达成了，这个时候的数据框长这样：

有7个变量，91040个观测。在做地图的过程中，只需要经纬度和分组这三个变量即可，就是long、lat和group。其实看了完整的数据形式基本就很容易理解这个画图的原理的，很简单，我将做出来的地图放大，发现每个省的图形，就是一个个的闭合多边形，这些多边形就是一条条线段围成的，这些线段的端点，就是这一个个的经纬度，按照什么顺序围成的那一个个多边形，就是根据group这个分类变量，一共有925个组，也就是有925个多边形拼成了我们看到的地图，每个组内的那么多个经纬度点，就围成一圈了，如果线段足够多，趋近于无限的时候，那做出来的边界就非常的顺滑了，不过我们这9万多个点，足够用了。
理解这个原理很重要，因为现在做出来的地图仅仅是个地图的框架，要想做成热图，还需要连接上我们的数据，这个时候就能够知道如何把我们的数据跟这些地理信息数据结合再一起了。
好了，画图：

ggplot(china_map_frame,aes(x=long,y=lat,group=group))+
  geom_polygon(fill="white",colour="black")+
  coord_map("polyconic")+
  theme(
    panel.grid=element_blank(),
    panel.background=element_blank(),
    axis.text=element_blank(),
    axis.ticks=element_blank(),
    axis.title=element_blank(),
    legend.position=c(0.2,0.3)
  )

解释一下，ggplot中的三个映射参数：x，y和group就对应china_map_frame中的经纬度和分组三个变量。接着用几何对象函数geom_polygon()，就是多边形函数，将数据映射到图上，函数中的fill和colour参数就是多边形边框和填充。coord_map()函数就是将上面的坐标系转换成地理坐标系，并且使用polyconic投影方式，将根据经纬度围出来的多边形拉到二维上。再下面的theme()函数就是主题函数，其中的那几个参数的意思就是将图上的图例、背景板、坐标轴等东西全都设为无，那图上就仅剩地图了。就长这样：

第四步：连接自己的数据

要连接数据，就是将两个数据框合并起来，就需要一个公共变量作为桥梁。现在我们想做的是各个省份的热图，那么省份的名字必然就是我们要找的公共变量，当然，也可能在地理信息数据里省份的名字用拼音，不过不影响。我们在china_map，也就是最初导入的那个地理信息数据中，找到了可以用来做桥梁的变量，省份的名称，而且还是汉字。

interim <-china_map@data

就是用这个语句从china_map中提取出来的，这个interim数据框长这样：

7个变量，925个观测。很容易就想到了，他是对应上面那个china_map_frame数据框中group变量的925个分类的，不幸的是，这个interim里没有经纬度和分组数这三个用来画地图的变量，所以就要考虑两次合并了第一次是将这个有省份名字的数据框与有那三个地理信息的数据框合并，然后再根据省份名将我们的数据合并进去。思路搞好了，那就分两步走：
第一步：将有地理信息数据的数据框china_map_frame和有省份名称的数据框interim合并，但是有个问题，china_map_frame和interim中并没有公共变量作为桥梁，所以需要自己建立一个公共变量，我找了一些资料，发现新建公共变量的方法都一样，不过并没有人解释为什么，或许在地理信息数据里面，interim中各个多边形块排列的顺序就是china_map_frame里的group对应的每个多边形块排列的顺序。所以，给interim新加一个id变量，数值就和china_map_frame中的id一样，从0到924，一共925个观测。然后就可以根据id变量，将两个数据框合并起来，这里使用jion()函数，依据id变量，合并类型为full。生成一个完整可用的中国地图数据china_map_data。

interim <-data.frame(interim,id=seq(0:924)-1)
china_map_data<-join(interim,china_map_frame,type="full")

第二步就是将我们的数据与上面这个完整的中国地图数据合并起来，不过首先要知道我们的数据长什么样（点击这里下载，提取码：vtre）：

就是这个样子，我们的数据导入后命名为my_data，其中有个NAME（注意大写）变量，其中的省份就是和interim中的NAME一样的。因为完整的中国地图数据china_map_data中也有了NAME变量，所以这个时候就可以将完整的地图数据和我们的数据合并了。

my_data<-read.csv("data_dt.csv",header=T,as.is=T)
china_heatmap_data <- join(china_map_data, my_data, type="full")

这个时候的热图数据就已经准备好了china_heatmap_data。最终的她长这样：

16个变量，91040个观测，数据信息很多，但我们用到的只有4个：经纬度（long，lat），分组（group）和我们热力图的数据（ratio等自己需要的数据）。

第五步：生成热图

数据的清洗与准备永远是数据操作中最困难，最耗时的部分，接下来的操作，只要按照代码规规矩矩的走就可以了。

ggplot(china_heatmap_data,aes(x=long,y=lat,group=group,fill=ratio))+
  geom_polygon(colour="grey40")+
  scale_fill_gradient(low="white",high="red")+
  coord_map("polyconic")+
  theme(
    panel.grid=element_blank(),
    panel.background=element_blank(),
    axis.text=element_blank(),
    axis.ticks=element_blank(),
    axis.title=element_blank(),
    legend.position='right'
  )

其中，ggplot的映射参数fill就是我们需要反映在热图上的数据ratio；几何对象函数geom_polygon()中的colour参数对应地图上各个省边框的颜色；度量函数scale_fill_gradient()函数是填充渐变函数，其中的两个参数low和high对应渐变色的最高和最低色；主题参数theme()中的legend.poistion参数是调整图例的位置，不过大可不必一点点的调整，因为在图片导出之后，图上的所有元素都是可以再调整的，总之，做出来的图长这样：

如果需要添加图例什么的，可在theme()主题函数中调整。或者导出为esp图片格式之后，在Ai中自行添加。

到现在热图就已经做出来了，下面就是如何在这个热图的基础上添加扇形图。

二、复合气泡饼图

饼图能够反映数据多个变量的比例，在我们是示例图中，要想做出来在热图上绘制饼图，需要的信息是：饼图的位置、饼图半径的映射数据、几个变量所占的比例。知道了这些个数据，那么这个饼图就非常的容易画出来了。只要再添加一个scatterpie几何映射函数就可以了，因为有上面的基础，我尝试敲了代码……
在网上找到了各个省会的经纬度（点击这里下载，提取码：89an）（但是这里只有31个省会的经纬度数据，不过影响不大）。稍微查看了scatterpie函数的用法，可以了解到他需要的变量数据：经纬度、分类变量（每个饼图的多个扇形为一类，在这里就是省会名称）、饼图中每个扇形的数据（一个饼图中有几个扇形就有几个变量）、饼图大小对应的数据。这时候我就随便搞些数据了：

第一步：制作数据

long_lat_data <- read.csv("long_lat_data.csv",header=T,as.is=T)

n <- nrow(long_lat_data)

pie_data <- data.frame(
  NAME = long_lat_data$NAME,
  radius = abs(rnorm(n,sd=2)),
  A = abs(rnorm(n,sd=1)),
  B = abs(rnorm(n,sd=2)),
  C = abs(rnorm(n,sd=3)),
  D = abs(rnorm(n,sd=4)),
  stringsAsFactors = F
)

pie_data_complete <- join(pie_data,long_lat_data)

最终的pie_data_complete就是我们要用到的数据，长这样：

9个变量，31个观测，这里的ABCD就是我们饼图的四个扇形对应的变量，radius变量就是映射气泡饼图的半径（其中的name参数下面会说）。

第二步：画图

结合我们做热力地图时的变量和代码，添加上一个新的scatterpie几何映射函数，很容易就写出来下面的代码：

ggplot()+
  geom_polygon(data=china_heatmap_data,colour="grey40",
               aes(x=long,y=lat,group=group,fill=ratio))+
  scale_fill_gradient(low="white",high="red")+
  geom_scatterpie(data=pie_data_com,
                  aes(x=long,y=lat,group=NAME,r=radius/1.5),
                  cols=names(pie_data_com)[3:6],color=NA)+
  coord_map("polyconic")+
  theme(
    panel.grid=element_blank(),
    panel.background=element_blank(),
    axis.text=element_blank(),
    axis.ticks=element_blank(),
    axis.title=element_blank(),
    legend.position='right'
  )

奈何自己是个写bug小能手：

错误: Discrete value supplied to continuous scale

运行就报错了。经过反复的检查测试这几行，发现这个错误说的离散型变量用在了连续性变量上说的就应该是热力地图中的fill参数和气泡饼图中的cols参数冲突。如果把热力地图中的fill参数去掉，对应的fill参数的度量函数也去掉，然后运行之后他就长这样：

有点丑。。。并且还能发现，因为在做热力地图的时候，加入了地图投影方式的图层，这里面的气泡扇形图也变样了，在这个坐标系下，也变形了。
所以我有个想法，是不是这两种对象（热力地图和气泡饼图）本来就不是在一张图上画出来的？又或者是有其他更神奇的代码或者软件，不过这些短时间还找不到答案，所以我还是按照自己的想法，将这两个图分别画出来，然后将他们放在一起就行了。

第三步：拼图

到这里我们的新思路就出来的：示例图是由两种元素组合成的，一个是热力地图，另一个是气泡饼图，既然我们有气泡饼图的经纬度数据，那么在合成的时候微调一下位置就行了，或者说在气泡饼图上添加省份名称的图层，在合并的时候就不会出现调错位置的情况。分三步走：
第一步，画出单独的热力地图素材，上面已经画好了，跳过。
第二步，画出单独的气泡饼图素材：

ggplot()+
  geom_scatterpie(data=pie_data_com,aes(x=long,y=lat,group=NAME,r=radius/1.5),
                  cols=names(pie_data_com)[3:6],color=NA)+
  geom_text(data=pie_data_com,aes(x=long,y=lat,label=name),colour="black",size=2)+
  geom_scatterpie_legend(pie_data_com$radius/1.5, x=90, y=20)

其中有个图例映射函数geom_scatterpie_legend()，用于映射气泡的大小，geom_text()函数是文本映射函数，将文本放到指定的位置，但是不能用中文。。。会乱码，所以就用了name变量，都是拼音。画出来这样的图：

现在素材都准备好了，下面就开始合并了！
第三步，合并：这时候需要把我们做的图导出成ESP格式的图片，就是矢量图，并且图中的每个元素都是可编辑的，编辑的软件就是Ai（Adobe Illustrator）。因为是位置的调整，并不涉及数据映射的改变，所以不违反作图的原则。
然后他就长这样了：

忽略配色什么的，还有那几个大得可怕的气泡饼图。。。只是数据没找好的原因，收工。