用R软件绘制中国分省市地图
【注】新版本的maptools包对很多函数进行了修改，对于修改的内容，文章中用红色的文字进行了说明。
鉴于最近有不少人在讨论用R软件绘制地图的问题，我也就跟着凑了凑热闹，对相应的方法学习了一番。下面的这篇文章是一个初步的介绍，还有很多内容仍在学习和探索中，如果大家有什么意见或建议，我将根据自己学习的情况对文章进行进一步的补充。
在R中绘制地图其实是十分方便的，最直接的办法大概就是安装maps和mapdata这两个包，然后输入下面的命令：
library(maps)
library(mapdata)
map("china")
其中map()函数还可以加上很多参数，在这里就不一一详述，具体的用法只需问号之。然而仔细看一看这张地图你会发现重庆市和四川省仍然是浑然一体，可见该地图的数据应该是有些年头了。
幸运的是，通过谢益辉的这篇博文我们已经可以大体知道该如何操作了，下面就为大家介绍一下具体的步骤。
首先，从这里下载中国地图的GIS数据，这是一个压缩包，完全解压后包含三个文件（bou2_4p.dbf、bou2_4p.shp和bou2_4p.shx），将这三个文件解压到同一个目录下，并在R中设好相应的工作空间，然后安装maptools包，运行如下程序：
library(maptools);
x=read.shape('bou2_4p.shp');#下文中会继续用到x这个变量，
#如果你用的是其它的名称，
#请在下文的程序中也进行相应的改动。
plot(x);
【修改】新版本的maptools包不再提供read.shape()函数，请用readShapePoly()代替。
这时一张完整的中国地图就已经画好了。但是在实际使用的过程中，我们往往会根据自己的需要对地图中的某些省份着以特定的颜色，这时就可以通过调节plot命令中的fg参数来予以实现。然而为了清楚地说明这部分的内容，我需要插播一段R绘制地图的原理。
======================传说中的分割线=====================
在绘制地图时，每一个省市自治区或者岛屿都是用一个多边形来表示的。之前的GIS数据，其实就是提供了每一个行政区其多边形逐点的坐标，然后R软件通过顺次连接这些坐标，就绘制出了一个多边形区域。在上面的数据中，一共包含了925个多边形的信息，之所以有这么多是因为一些省份有很多小的附属岛屿。在这925个多边形中，每一个都对应一个唯一的ID，编号分别从1到925。
======================传说中的分割线=====================
回到刚才的话题，plot命令中的fg参数在本例中应该是一个长度为925的向量，其第i个分量的取值就代表了地图中第i个多边形的颜色。一个简单的尝试是运行下面这个命令看看效果：
plot(x,fg=gray(924:0/924));
【修改】新版本的maptools包的绘图参数也有所改变，请将fg换成col。
于是自然就产生了一个问题：如何获取某一个特定地区的ID，进而设置我们想要的颜色？事实上，在变量x中，就已经存储了我们想要的信息。在R中输入“x[[2]]”或“x$att.data”，会得到一个925行7列的数据框，这其实是bou2_4p.dbf这个文件中存储的信息，之前的read.shape()函数虽然读取的是bou2_4p.shp文件，但在默认情况下会把dbf文件的信息也放到变量之中。对于这个数据框，其行名就是每一个区域的ID编号，第一列和第二列分别是面积和周长，最后一列是该区域所属的行政区名，其它的列应该也是一些编号性质的变量。于是，通过查找相应的行政区对应的行名，就可以对fg参数进行赋值了。下面是我编的一个函数，用来生成所需的fg向量：
getColor=function(mapdata,provname,provcol,othercol)
{
f=function(x,y) ifelse(x %in% y,which(y==x),0);
colIndex=sapply(mapdata$att.data$NAME,f,provname);
fg=c(othercol,provcol)[colIndex+1];
return(fg);
}
【修改】地图数据的组织形式有所变化，上面函数中的mapdata$att.data$NAME需要替换为mapdata@data$NAME。
其中mapdata是存放地图数据的变量，在上面的例子中就是x，provname是需要改变颜色的地区的名称，provcol是对应于provname的代表颜色的向量（名称和数字均可），othercol是其它地区的颜色。举例如下：
provname=c("北京市","天津市","上海市","重庆市");
provcol=c("red","green","yellow","purple");
plot(x,fg=getColor(x,provname,provcol,"white"));
注意provname一定要写地区的全称，写法可以参照下面这条命令生成的向量：
as.character(na.omit(unique(x$att.data$NAME)));
由此生成的向量有33个元素，少了澳门特别行政区，这是这个数据中的一块瑕疵。在x$att.data的第899行有一个NA，不知道它代表的是否就是澳门。
利用类似的方法就可以根据自己的需要对不同的区域进行着色，下面再举一例。从国家统计局获取2007年我国各地区的人口数据，然后根据人口的多少对各省份进行着色。程序如下：
provname=c("北京市","天津市","河北省","山西省","内蒙古自治区",
"辽宁省","吉林省","黑龙江省","上海市","江苏省",
"浙江省","安徽省","福建省","江西省","山东省",
"河南省","湖北省","湖南省","广东省",
"广西壮族自治区","海南省","重庆市","四川省","贵州省",
"云南省","西藏自治区","陕西省","甘肃省","青海省",
"宁夏回族自治区","新疆维吾尔自治区","台湾省",
"香港特别行政区");
pop=c(1633,1115,6943,3393,2405,4298,2730,3824,1858,7625,
5060,6118,3581,4368,9367,9360,5699,6355,9449,
4768,845,2816,8127,3762,4514,284,3748,2617,
552,610,2095,2296,693);
provcol=rgb(red=1-pop/max(pop)/2,green=1-pop/max(pop)/2,blue=0);
plot(x,fg=getColor(x,provname,provcol,"white"),xlab="",ylab="");
其中颜色越深的地方代表人口数越多，反之为人口数越少。
此外，在绘制地图的过程中，还有一个比较有用的参数是recs，它是一个由多边形ID组成的向量，表示在地图中只画出这些ID所代表的区域。利用这个参数，就可以画出某一部分的地图，例如下面的例子是我国中部六省的地图：
getID=function(mapdata,provname)
{
index=mapdata$att.data$NAME %in% provname;
ids=rownames(mapdata$att.data[index,]);
return(as.numeric(ids));
}
midchina=c("河南省","山西省","湖北省","安徽省","湖南省","江西省");
plot(x,recs=getID(x,midchina),fg="green",ol="white",xlab="",
ylab="");
上面的getID()是我编写的一个功能与getColor()类似的函数，用来返回指定省份的ID。
【修改】新版本的maptools包的绘图函数已经取消了recs这个参数，现在要实现这个功能，可以在颜色上把不需要的省份变成白色，其中填充色用col参数，边界颜色用border参数。例如上面的例子可以用下面的函数来实现：
plot(x, col = getColor(x, midchina, rep("green", 6),
"white"), border = "white", xlab = "", ylab = "")
最后要说的是，在画出的图上仍然可以用points()函数和text()函数加上点和文字，而maptools包中还提供了一个pointLabel()函数，用来解决文本标签的重叠问题。这部分内容请参阅博文：用R画中国地图并标注城市位置，以及避免文本标签重叠：maptools中的pointLabel()。
从以上的内容来看，本文所述的都是一些最基本的绘图方法，还没有对地理信息数据进行更进一步的分析。如果有机会的话，这一主题的下一篇文章将为大家介绍地图数据的组成结构，并说明如何将不同格式的地理数据整合起来，例如如何在上面的地图上绘制出我国的铁路、水系分布等内容。
很好很强大，以前用maptools，发现里面的数据很支离破碎，无法画分省图。
这下可好了。
顶了，这篇文对空间统计分析很有用哈~
如果可以对地图上的距离长度信息进行提取的话那么就更完美了！
假定地球为半径为R的球体，只要知道地图上两点的坐标（经纬度），就可以算出两点之间的大圆弧长，也就是最短距离。
太好了！我太喜欢了！！！！
以前要做个报告想把中国地图放进去，然后把不同的省份运用不同的颜色表示出来简直就是一件非常让人头疼得事情。所以就扫描然后自己ps处理。可是遇到需要外文标注的时候就异常的痛苦！！所以这次总算有了解决的方案！
多谢多谢！！
对不起，刚才一激动，感叹号居然超过了3个，犯了谢老大的版规。自我反省一下。
“再激动感叹号也不要超过三个！！”
呵呵，这个不绝对
建议补充上澳门的数据，否则不好在学术或正式场合引用。
很好,很和谐,很强大!!!
不错不错，以前我做的时候用了很麻烦的笨方法，现在好了，语句简化但用途更强大了：）
导入GIS数据后，发现省名在R里显示是乱码，请问应该如何解决？ 谢谢！
我也遇到了与bo一样的问题。应该是语言设置的问题，当我读入x的时候发现：
Shapefile type: Polygon, (5), # of Shapes: 925
Warning message:
use readShapeSpatial:
objects other than Spatial objects defined in the sp package are deprecated
之后也就没有办法对每个省的颜色进行改变！
正在想法解决
R是中文的还是英文的？
没有加载到GIS数据所在的工作空间
我试着运行程序，提示
> x<-read.shape('bou2_4p.shp');
错误: 没有"read.shape"这个函数
求解答，先谢了
新版本的maptools包把很多原来的函数都改掉了，如果你要用的话可能得重新读一下它的帮助。
谢谢，修改后就没问题了
用
readShapePoly()函数。如x = readShapePoly('bou2_4p.shp')邱怡轩得更新正文了。我刚看了一下，新的函数读进来的数据的子元素直接就是
NAME那一级的，原来的x$att.data$NAME就是现在的x$NAME，比如plot(x, col = rainbow(length(levels(x$NAME)))[as.integer(x$NAME)])哈哈。太感谢了。终于成功了。我用的是2.10.1
plot(x,fg=gray(924:0/924))
错误于plot.window(xlim = xlim, ylim = ylim, asp = asp, …) :
图形参数”fg”的长度不对
参见你头上我的回复中的例子，这个包以及
sp包都更新了，现在应该用col参数。这个是怎么回事？
> getID=function(mapdata,provname)
+ {
+ index=mapdata$att.data$NAME %in% provname;
+ ids=rownames(mapdata$att.data[index,]);
+ return(as.numeric(ids));
+ }
> midchina=c(“河南省”,”山西省”,”湖北省”,”安徽省”,”湖南省”,”江西省”);
> plot(x,recs=getID(x,midchina),col=”green”,ol=”white”,xlab=”",
+ ylab=”");
警告多于50个(用warnings()来显示第一个到第五十个)
>
新版本的
maptools包改动很大，我已经把文章更新了一下，参见文中的红色部分。呵呵 getColor 好像找不到了
什么意思？
getColor()不是写在前面吗？只是需要把里面的mapdata$att.data$NAME改成mapdata@data$NAME。恩。找到了 。谢谢
有两个问题期待大家帮助：
1、原文“首先，从这里下载中国地图的GIS数据，这是一个压缩包，完全解压后包含三个文件（bou2_4p.dbf、bou2_4p.shp和bou2_4p.shx），将这三个文件解压到同一个目录下，并在R中设好相应的工作空间，然后安装maptools包，运行如下程序：”
请问“工作空间”指的是什么？如何设置？
2、这个地图能精细到地级市级别吗？是不是在“mapdata”包里可以看到都有什么样的数据？可我找不到包的路径。
1、用R软件时可能经常要读取一些外部文件，引用这些文件有两种方式，一种是绝对路径，比如“C:/AAA.txt”，另一种就是把工作空间设在C:/，然后只要用相对路径“AAA.txt”就行了。说白了就是少写了一大串路径名。R中的函数是
setwd()。2、地图的精确度与数据有关，
mapdata里自带了一些，如果要更精确的，可以到网上去搜GIS数据。要读取mapdata里的数据可以这样做：cnmapdata=map("china")，然后cnmapdata$x和cnmapdata$y就分别是经度和维度了。谢谢，弄明白了。
那个工作空间，原来R for beginners里就有提到。