用R画地图数据
首先,从这里下载中国地图的GIS数据,这是一个压缩包,完全解压后包含三个文件(bou2_4p.dbf、bou2_4p.shp和bou2_4p.shx),将这三个文件解压到同一个目录下。
用R绘制地图比较简单。比如画一下全国范围的区域,可以用如下代码:
library(maptools) mydat = readShapePoly("china-province-border-data.tar/china/bou2_4p.shp") #地图包位置,根据自己的角标位置设置 plot(mydat)
生成的图如下:
但是,可以看出这样绘制的地图的形状有些扁平。这是因为,在绘图的过程中,默认把经度和纬度作为普通数据,均匀平等对待,绘制在笛卡尔坐标系上造成的。其实,地球的球面图形如何映射到平面图上,在地理学上是有一系列不同的专业算法的。地图不应该画在普通的笛卡尔坐标系上,而是要画在地理学专业的坐标系上。
也可以安装maps
和mapdata
这两个包,然后输入下面的命令:
install.packages(mapdata) install.packages("mapdata") library(maps) library(mapdata) map("china")
生成的图如下:
其中map()
函数还可以加上很多参数,大致如下:
map(database = "world", regions = ".", exact = FALSE, boundary = TRUE, interior = TRUE, projection = "", parameters = NULL, orientation = NULL, fill = FALSE, col = 1, plot = TRUE, add = FALSE, namesonly = FALSE, xlim = NULL, ylim = NULL, wrap = FALSE, resolution = if (plot) 1 else 0, type = "l", bg = par("bg"), mar = c(4.1, 4.1, par("mar")[3], 0.1), myborder = 0.01, ...)
可以设置数据库、地区、精确度、边界等,在这里就不一一详述,具体的用法可以?map。
给地图上色
在实际使用的过程中,我们往往会根据自己的需要对地图中的某些省份着以特定的颜色,这时就可以通过调节plot
命令中的fg
参数来予以实现。
首先看看R绘制地图的原理:
在绘制地图时,每一个省市自治区或者岛屿都是用一个多边形来表示的。之前的GIS数据,其实就是提供了每一个行政区其多边形逐点的坐标,然后R软件通过顺次连接这些坐标,就绘制出了一个多边形区域。在上面的数据中,一共包含了925个多边形的信息,之所以有这么多是因为一些省份有很多小的附属岛屿。在这925个多边形中,每一个都对应一个唯一的ID,编号分别从1到925。
plot
命令中的fg
参数在本例中应该是一个长度为925的向量,其第i个分量的取值就代表了地图中第i个多边形的颜色。一个简单的尝试是运行下面这个命令看看效果:
plot(x,col=gray(924:0/924));
生成的图如下所示:
于是自然就产生了一个问题:如何获取某一个特定地区的ID,进而设置我们想要的颜色?事实上,在变量x中,就已经存储了我们想要的信息。在R中输入“x[[2]]
”或“x$att.data
”,会得到一个925行7列的数据框,这其实是bou2_4p.dbf这个文件中存储的信息,之前的read.shape()
函数虽然读取的是bou2_4p.shp文件,但在默认情况下会把dbf文件的信息也放到变量之中。对于这个数据框,其行名就是每一个区域的ID编号,第一列和第二列分别是面积和周长,最后一列是该区域所属的行政区名,其它的列应该也是一些编号性质的变量。于是,通过查找相应的行政区对应的行名,就可以对fg
参数进行赋值了。下面是我编的一个函数,用来生成所需的col向量:
getColor=function(mapdata,provname,provcol,othercol) { f=function(x,y) ifelse(x %in% y,which(y==x),0); colIndex=sapply(mapdata@data$NAME,f,provname); col=c(othercol,provcol)[colIndex+1]; return(col); }
其中mapdata
是存放地图数据的变量,在上面的例子中就是x,provname
是需要改变颜色的地区的名称,provcol
是对应于provname
的代表颜色的向量(名称和数字均可),othercol
是其它地区的颜色。举例如下:
provname=c("北京市","天津市","上海市","重庆市"); provcol=c("red","green","yellow","purple"); plot(x,col=getColor(x,provname,provcol,"white"));
生成的图数据如下:
生成人口数据分布图
利用类似的方法就可以根据自己的需要对不同的区域进行着色,下面再举一例。从国家统计局获取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,col=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, col = getColor(x, midchina, rep("green", 6), "white"), border = "white", xlab = "", ylab = "")
生成的图:
用R画中国地图并标注城市位置
画出的图上仍然可以用points()
函数和text()
函数加上点和文字,而maptools
包中还提供了一个pointLabel()
函数,用来解决文本标签的重叠问题。
par(mar=rep(0,4)) dat = read.csv(text = "城市,jd,wd 北 京,116.4666667,39.9 上 海,121.4833333,31.23333333 天 津,117.1833333,39.15 重 庆,106.5333333,29.53333333 哈尔滨,126.6833333,45.75 长 春,125.3166667,43.86666667 沈 阳,123.4,41.83333333 呼和浩特,111.8,40.81666667 石家庄,114.4666667,38.03333333 太 原,112.5666667,37.86666667 济 南,117,36.63333333 郑 州,113.7,34.8 西 安,108.9,34.26666667 兰 州,103.8166667,36.05 银 川,106.2666667,38.33333333 西 宁,101.75,36.63333333 乌鲁木齐,87.6,43.8 合 肥,117.3,31.85 南 京,118.8333333,32.03333333 杭 州,120.15,30.23333333 长 沙,113,28.18333333 南 昌,115.8666667,28.68333333 武 汉,114.35,30.61666667 成 都,104.0833333,30.65 贵 阳,106.7,26.58333333 福 州,119.3,26.08333333 台 北,121.5166667,25.05 广 州,113.25,23.13333333 海 口,110.3333333,20.03333333 南 宁,108.3333333,22.8 昆 明,102.6833333,25 拉 萨,91.16666667,29.66666667 香 港,114.1666667,22.3 澳门,113.5,22.2") library(maps) library(mapdata) map("china", col = "darkgray", ylim = c(18, 54), panel.first = grid()) points(dat$jd, dat$wd, pch = 19, col = rgb(0, 0, 0, 0.5)) text(dat$jd, dat$wd, dat[, 1], cex = 0.9, col = rgb(0, 0, 0, 0.7), pos = c(2, 4, 4, 4, 3, 4, 2, 3, 4, 2, 4, 2, 2, 4, 3, 2, 1, 3, 1, 1, 2, 3, 2, 2, 1, 2, 4, 3, 1, 2, 2, 4, 4, 2)) axis(1, lwd = 0); axis(2, lwd = 0); axis(3, lwd = 0); axis(4, lwd = 0)
生成的图:
ggplot2绘制地图
以中国地图为例,下载最新的ArcGIS矢量地图数据,这种地图数据包含了很多信息,这是画地图的基础数据。下载地址:
http://download.csdn.net/detail/lgstarzkhl/9427677
用以下代码进行地图绘制:
library(maptools) library(ggplot2) library(plyr) #读取地图文件 china_map<-readShapePoly("C:/Users/feng/Desktop/chinaprovinceborderdata_tar_gz/china-province-border-data.tar/Lambert/省级行政区.shp") #提取用于绘图的地图数据 x<-china_map@data xs<-data.frame(x,id.1=seq(0:33)-1) #将地图数据转换为数据框 china_map1<-fortify(china_map) #添加一个id.1字段,用于和上面的xs(各省市数据)糅合,合并 china_map1$id.1<-china_map1$id #去掉china_map1中的id字段,避免在糅合数据的时候,出现两个相同字段id和id.1,保证只用id.1来糅合 china_map2<-china_map1[,-7] #糅合地图数据 china_mapdata<-join(china_map2, xs, type = "full") #绘制地图 ggplot(china_mapdata, aes(x = long, y = lat, group = group,fill=NAME))+ geom_polygon( )+ geom_path(colour = "grey40")+ scale_fill_manual(values=colours(),guide=FALSE)
生成的图:
总结
使用R的地图扩展包可以画出各种样式的地图,具体的展现形态及结合方式还有待进一步挖掘。