标签：
一个合格的数据分析师必定能够从大量数据中洞察需求，并把数据和需求展示给团队人员，获得资源支持。数据可视化技术对于产品经理来说也是相当重要的一项技能，通过数据可视化技术我们可以把数据分析结果以“人性化”、“友好”的方式反馈给各个合作方，使接收者能够较为愉快地了解整个行业、产品或其他方面的相关发展状况，从而有利于项目的顺利进行。在调研过程中我们经常会遇到与中国地区区域有关的数据显示问题，今天笔者在此讨论一下如何通过R软件来绘制中国地图及与之相关的类似于热力图/密度图等相关内容。
1.绘制中国地图
要想绘制中国地图，首先想到的是在环境中加载maps和mapdata这两个包，R软件安装时是默认不安装这两个包的。代码：
>> path.package()
可以用来查看当前加载了哪些包，以及这些包的存放路径。接下来我们下载安装maps和mapdata这两个包，相关的代码如下，注意路径选择要根据自己R软件的安装情况，服务器要选择China-Beijing的服务器：
>> install.packages("maps","C:/ProgramFiles/R/R-3.2.3/library")
>> install.packages("mapdata","C:/ProgramFiles/R/R-3.2.3/library")
安装完之后我们需要将其加载到环境中，相关代码是这样的：
>> library(maps)
>> library(mapdata)
完成加载包之后我们就可以绘制中国地图了，代码和绘制出来的中国地图如下：
>> map("china")
我们不仅可以绘制中国地图我们还可以绘制美国地图、日本地图、世界地图等等，如下：
>> map("world")
如果我们仔细观察一下上面绘制的中国地图的话，可以发现四川省和重庆市是一体的，可见该数据比较老，为此我们重新下载中国地图的GIS数据（http://pan.baidu.com/s/1pLOSmM3），该数据包含三个文件：bou2_4p.dbf、bou2_4p.shp和bou2_4p.shx。数据来源于国家基础地理信息中心早年公开的数据。把这三个文件放置到当前的工作目录中，以下代码用来获取当前的工作目录：
>> getwd()
我们还需要用到readShapePoly()函数,因此我们应用刚才的方法加载maptools包：
>> install.packages("maptools","C:/ProgramFiles/R/R-3.2.3/library")
>> library(maptools)
加载完成之后我们就可以直接使用该包中所包含的一些与地图有关的函数。下面我们重新绘制中国地图：
>> map_data -> readShapePoly(‘bou2_4p.shp‘)
>> plot(map_data)
可以看出此地图上四川省和重庆市已然是两个单位。在此我们有必要讨论一下刚才绘制地图的原理，事实上，地图上的每一个区域（包括省、直辖市、特别行政区、岛屿等等）都是用多边形来表示的，所谓的GIS数据，实际上包含了构成每个区域的多边形的顶点坐标，通过plot函数顺次连接这些坐标就构成了地图。map_data中包含了925组数据，之所以有这么多是因为我国有许多零零散散的小岛。这925个多边形每个都有唯一的一个ID来对应。
-------------------------------------------------------------cut-----------------------------------------------------------------
探索map_data:
>> names(map_data)
[1] "AREA" "PERIMETER" "BOU2_4M_" "BOU2_4M_ID" "ADCODE93"
[6] "ADCODE99" "NAME"
可以看出map_data中有7列，对应的字段名如上面显示。
>> map_data$AREA #925个区域单元的面积
>> map_data$PERIMETER #925个区域单元的周长
>> map_data$BOU2_4M_ #没有重复的数字，2~926，可作为区域单元ID
>> map_data$BOU2_4M_ID #有重复数字，特定情况下可作为区域单元ID
>> map_data$ADCODE93 #93版ADCODE地理编码
>> map_data$ADCODE93 #99版ADCODE地理编码
>> map_data$NAME #各区域单元所隶属的省级行政单元的名称
>> uinique(map_data$NAME) #查看各区域的名称是什么文本
[1] 黑龙江省 内蒙古自治区 新疆维吾尔自治区
吉林省
[5] 辽宁省 甘肃省 河北省 北京市
[9] 山西省 天津市 陕西省 宁夏回族自治区
[13] 青海省 山东省 西藏自治区 河南省
[17] 江苏省 安徽省 四川省 湖北省
[21] 重庆市 上海市 浙江省 湖南省
[25] 江西省 云南省 贵州省 福建省
[29] 广西壮族自治区 台湾省 广东省 香港特别行政区
[33] 海南省 <NA>
33 Levels: 安徽省北京市
福建省甘肃省 广东省广西壮族自治区 ...重庆市
-------------------------------------------------------------cut-----------------------------------------------------------------
2.为地图着色
>> plot(map_data,col="red") #整个地图都变成了红色，当然也可以变成yellow,green等颜色
我们也可以为整个地图涂上深浅不一的颜色，用gray()函数。
>> plot(map_data,col=gray(1:925/925))
刚才是为整个地图着色，然而现实中我们经常遇到为特定的省份着色，或者根据用户数量密度不同在不同地区涂上深浅不同的颜色。下面我们为不同的省份涂色。首先创建一个getID函数，该函数可通过省份名称获得包含于该省份的区域的所有ID。代码如下：
getID <- function(mapData,prov_vec){
ID_vec <- NULL
for (i in (1:length(prov_vec))) {
NAME_match <- which(mapData$NAME == prov_vec[i])
ID_vec <- c(ID_vec,mapData$BOU2_4M_[NAME_match])
}
return(ID_vec)
} #由省份名称得到包含于该省份的所有ID
mapData即为地图数据，在这里是map_data；prov_vec为省份名称向量，如c("河北省","北京市")。假如我们想获得包含于河北省、北京市的所有区域ID，我们只需要：
>> getID(map_data,c("河北省","北京市"))
[1] 8 30 70 73 92 94 96 107 121 9
可以看出河北省和北京市总共包含了10个区域，笔者猜想应该是河北省包含了些许小岛或海岸区域。具体河北省包含了几个区域、北京市包含了几个区域通过该代码是无法获取的。
紧接着我们获取为地图着色所需要用到的颜色向量，构建getcols函数，该函数可以获得你想要涂色的颜色向量：
getcols <- function(prov_vec,col_vec,other_col){
col_list <- NULL
IDD <- getID(map_data,prov_vec)
for (i in (1:length(map_data$BOU2_4M_))){
if (map_data$BOU2_4M_[i] %in% IDD)
col_list[i] <- col_vec[which(prov_vec==map_data$NAME[i])]
else col_list[i] <- other_col
}
return(col_list)
} #得到颜色向量
该函数的作用是：我们输入省份向量（长度为自定义，根据需要填写省份）以及对应的颜色向量（长度跟省份向量长度一致并且对应），则可以获取涂色所需要的整个颜色向量（即为plot函数中所用到的向量，长度为925）。其中prov_vec和函数getID中的prov_vec意义一样；col_vec为颜色向量，如c("yellow","red")；other_col为除了涂色目标省份之外其他区域要涂的颜色，如"white"。假如现在我们分别想为河北省和北京市涂上黄色和红色，则代码只需按如下写：
>> plot(map_data,col=getcols(c("河北省","北京市"),c("yellow","red"),"white"))
3.人口密度图的例子
笔者获取了中国2010年第六次人口普查的数据（虽然由于国情问题数据不包含港澳台地区，但是香港、澳门、台湾自古以来都是中国不可分割的一部分，下面地图中港澳台颜色被设置为white），如下：
不妨用灰色来表示密度，利用gray函数来画人口密度图。我们首先建立getAREA函数，该函数的功能是通过省份名称得到省份的面积，代码如下：
getAREA <- function(prov_vec){
AREA_vec <- NULL
for (i in (1:length(prov_vec))){
Each_AREA <- map_data$AREA[which(map_data$BOU2_4M_ %in%getID (map_data,prov_vec[i]))]
AREA_vec <- c(AREA_vec,sum(Each_AREA))
}
return(AREA_vec)
} #由省份名称得到对应的省份面积
其中的prov_vec是省份名称向量。假如现在我们想知道河北省和北京市的面积，我们只需要如下代码：
>> getAREA(c("河北省","北京市"))
[1] 19.636 1.733
根据地理常识可知数据的单位应该是万平方千米。我们也可以给出34个省级行政单位的面积：
>> getAREA(unique(map_data$NAME))
[1] 54.447 129.113 175.591 21.315 15.650 41.508 19.636 1.733 15.961
[10] 1.214 20.394 5.272 71.363 15.400114.331 16.135 9.741 13.365
[19] 45.534 17.563 7.716 0.592 9.439 19.390 15.277 34.276 15.975
[28] 10.962 20.941 3.193 15.598 0.089 2.903 0.000
注意：最后的0是数据缺陷。接下来我们输入作图所需的省份向量和人口向量（注意一定要对应）：
>> prov_name <- c("黑龙江省","内蒙古自治区","新疆维吾尔自治区","吉林省","辽宁省","甘肃省","河北省","北京市","山西省","天津市","陕西省","宁夏回族自治区","青海省","山东省","西藏自治区","河南省","江苏省","安徽省","四川省","湖北省","重庆市","上海市","浙江省","湖南省","江西省","云南省","贵州省","福建省","广西壮族自治区","广东省","海南省")
>> popul <-c(38313991,24706291,21815815,27452815,43746323,25575263,71854210,19612368,35712101,
12938693,37327379,6301350,5626723,95792719,3002165,94029939,78660941,59500468,80417528,
57237727,28846170,23019196,54426891,65700762,44567797,45966766,34748556,36894217,46023761,
104320459,8671485)
求出密度并且构造出col向量：
>> dens <- popul/getAREA(prov_name)
>> plot(map_data,col=getcols(prov_name,gray(1-dens/max(dens)),"white"))
从上图可以看出，颜色越深的地方表示人口密度越大。也可以验证一个常识，即：东部沿海地区的人口密度要比西北部的人口密度大。
在其他的许多分析中，我们没必要用到dens，绝对的数量就可以直接反映到地图上并且意义重大：
>> plot(map_data,col=getcols(prov_name,gray(1-popul/max(popul)),"white"))
我们也可以把颜色做的更加绚丽一点，利用rgb()函数：
>> plot(map_data,col=getcols(prov_name,rgb(red=1-popul/max(popul),green=1-popul/max(popul),blue=0),"white"))
>> plot(map_data,col=getcols(prov_name,rgb(red=1-popul/max(popul),green=0,blue=1-popul/max(popul)),"white"))
读者还可以为地图添加标注信息（名称、数值等），也可以举一反三研究更小区域的问题，希望读者可以做进一步的探索。随着中国乃至世界范围内的气候、环境、资源等问题的日益突出，对中国地图与可视化结合问题的探索早已变得很重要。
参考文章：Yixuan，2009，http://cos.name/2009/07/drawing-china-map-using-r/
标签：
原文地址：http://blog.csdn.net/zhenglit/article/details/51104332
