案例如何用R实现空间数据可视化.docx
《案例如何用R实现空间数据可视化.docx》由会员分享，可在线阅读，更多相关《案例如何用R实现空间数据可视化.docx（9页珍藏版）》请在冰豆网上搜索。
案例如何用R实现空间数据可视化
【案例】如何用R实现空间数据可视化.
流行病学的数据讲究“三间分布”，即人群分布、时间分布和空间分布。
其中的“空间分布”最好是在地图上展示，才比较清楚。
R软件集统计分析与高级绘图于大成，是最适合做这项工作了。
关于地图的绘制过程，谢益辉、邱怡轩和陈丽云等人都早有文章讲述，开R地图中文教程之先河。
由于目前指导毕业论文用到，因此研究了一下。
本来因为网上教程很多，曾打消了写些文字的计划，但怡轩版主鼓励说“教程者众，整合者鲜”，所以才战胜拖延症，提起拙笔综述整合一下，并对DIY统计GIS地图提出了一点自己的想法。
1地图GIS数据的来源与R绘制软件包中国地图GIS数据的官方数据可以在国家基础地理信息中心的网站（）里面可以免费下载。
官方公开的数据包括：
地图数据，及居住地、交通、河流等辅助数据。
今年6月开始，官方正组织开始制作新版数据。
老数据暂时无法下载，读者要自行XX搜索，本文以旧版数据为例。
旧版地图数据中部分地名和地市区划已经过时，使用时需注意。
地图数据有4个压缩文件：
bou1_4m.zip、bou2_4m.zip、bou3_4m.zip和bou4_4m.zip。
bou代表边界的意思，数字1~4代表国家、省、市、县的4级行政划分；4m代表比例是400万分之一，这个比例的图形是公开的。
每个文件解压缩后含有两类文件：
以字母p结尾的表示多边形数据，用来绘制区域；以字母l结尾的文件是线形数据，用来绘制边界。
但是老版数据中，市级数据中缺少绘制区域的多边形数据，让市级分布图的绘制稍麻烦一些，新版中也许会有改进。
用R绘制地图比较简单。
比如画一下全国范围的区域，可以用如下代码：
library（maptools）
mydat=readShapePoly（"maps/bou1/bou1_4p.shp"）
plot（mydat）但是，可以看出这样绘制的地图的形状有些扁平。
这是因为，在绘图的过程中，默认把经度和纬度作为普通数据，均匀平等对待，绘制在笛卡尔坐标系上造成的。
其实，地球的球面图形如何映射到平面图上，在地理学上是有一系列不同的专业算法的。
地图不应该画在普通的笛卡尔坐标系上，而是要画在地理学专业的坐标系上。
在这一点上，R的ggplot2包提供了专门的coord_map（）函数。
所以推荐R的ggplot2包来绘制地图。
library（ggplot2）
mymap=ggplot（data=fortify（mydat））+
geom_polygon（aes（x=long,y=lat,group=id）,colour="black",
fill=NA）+
theme_grey（）
print（mymap+coord_map（））这次中国地图的形状与XX地图一样了。
ggplot2包的coord_map函数默认的映射类型是mercator。
如果有其他需要，可以使用其他的映射类型来绘制地图，如：
mymap+coord_map（projection="azequidistant"）coord_map函数的映射类型及其含义可以通过下列代码查询帮助，一般我们用默认的就可以。
library（mapproj）
?
mapproject
2GIS地图的数据结构及省市地图的绘制GIS地图有很多种存储格式，其中shapefile格式（.shp）可以通过R的maptools包打开。
其他格式可以去R官网查询相应的软件包。
地图数据基本可以分为点、线、面三种数据，在maptools包内分别有对应的函数来读取（readShapePoints、readShapeLines和readShapePoly函数）。
首先以面（poly）型数据介绍。
先看代码，通过readShapePoly函数读入省级地图：
library（maptools）
mydat=readShapePoly（"maps/bou2/bou2_4p.shp"）此时，mydat中保存的是各个省/直辖市的多边形面图，数据类型是SpatialPolygonsDataFrame。
我们可以：
length（mydat）
##[1]925
names（mydat）
##[1]"AREA""PERIMETER""BOU2_4M_""BOU2_4M_ID""ADCODE93"
##[6]"ADCODE99""NAME"可以发现mydat中有925条记录，每条记录中含有面积（AREA）、周长（PERIMETER）、各种编号、中文名（NAME）等字段。
其中中文名（NAME）字段是以GBK编码的。
这个SpatialPolygonsDataFrame类型并不是真正的data.frame类型，而是一个sp包定义的类，只不过重载了[]和$运算符，使得一些行为上与data.frame相类似。
可以进一步统计一下，每个省/直辖市的多边形数目。
table（iconv（mydat$NAME,from="GBK"））
##
##上海市云南省内蒙古自治区北京市
##12111
##台湾省吉林省四川省天津市
##57111
##宁夏回族自治区安徽省山东省山西省
##11861
##广东省广西壮族自治区新疆维吾尔自治区江苏省
##154615
##江西省河北省河南省浙江省
##191179
##海南省湖北省湖南省甘肃省
##79111
##福建省西藏自治区贵州省辽宁省
##1681294
##重庆市陕西省青海省香港特别行政区
##11153
##黑龙江省
##1我的环境是UTF-8，所以需要iconv函数转化一下才能正常显示。
结果显示多数省的地图都是由一个多边形构成，少数临海省/直辖市由于有很多附属岛屿，多边形数目比较多。
利用与data.frame相似的[]和$运算符操作，我们可以迅速提取出一个省市的数据，比如上海及附属崇明岛：
Shanghai=mydat[mydat$ADCODE99==310000,]
plot（Shanghai）其中ADCODE99是国家基础地理信息中心定义的区域代码，共有6位数字，由省、地市、县各两位代码组成。
为了进一步在ggplot2包中绘图，需要把SpatialPolygonsDataFrame数据类型转化为真正的data.frame类型才可以。
ggplot2包专门针对地理数据提供了特化版本的fortify函数来做这个工作：
head（fortify（Shanghai））
##longlatorderholepiecegroupid
##1121.331.851FALSE1208.1208
##2121.331.852FALSE1208.1208
##3121.331.853FALSE1208.1208
##4121.331.854FALSE1208.1208
##5121.331.845FALSE1208.1208
##6121.431.836FALSE1208.1208
3在地图上展示流行病学数据3.1一地名对应一区域，长沙为例首先把长沙所辖地区找到，这个可以根据ADCODE99编码的前4位定位长沙，去查表就可以了。
但是这个地名是99年的标准，新版正在制定过程中，随时会变。
我们权且以此为例。
如果找不到表，可以通过代码在命令行下手工查找：
mydat=readShapePoly（"maps/bou4/BOUNT_poly.shp"）
tmp=iconv（mydat$NAME99,from="GBK"）
grep（"长沙",tmp,value=TRUE）
##[1]"长沙县""长沙市市辖区"
grep（"长沙",tmp）
##[1]21222183
mydat$ADCODE99[grep（"长沙",tmp）]
##[1]430121430101
##2368Levels:
0110100110112110113110221110224110226110227...820000这样我们就知道了长沙ADCODE99编码的前4位是4301，其中43代表湖南省，01就是长沙市。
接着就可以筛选出长沙的地图数据：
Changsha=mydat[substr（as.character（mydat$ADCODE99）,1,4）=="4301",]
mysh=fortify（Changsha,region='NAME99'）
mysh=transform（mysh,id=iconv（id,from='GBK'）,group=iconv（group,from='GBK'））
head（mysh,n=2）
##longlatorderholepiecegroupid
##1113.128.181FALSE1长沙市市辖区.1长沙市市辖区
##2113.128.182FALSE1长沙市市辖区.1长沙市市辖区
names（mysh）[1:
2]=c（"x","y"）#这句是不得已而为之的黑魔法接着我们给一串随机数当成是流行病学数据，并用颜色填充到地图上。
myepidat=data.frame（id=unique（sort（mysh$id）））
myepidat$rand=runif（length（myepidat$id））
myepidat
##idrand
##1宁乡县0.98076
##2望城县0.32123
##3浏阳市0.66957
##4长沙县0.09655
##5长沙市市辖区0.19437
csmap=ggplot（myepidat）+
geom_map（aes（map_id=id,fill=rand）,color="white",map=mysh）+
scale_fill_gradient（high="darkgreen",low="lightgreen"）+
expand_limits（mysh）+coord_map（）
print（csmap）接下来的工作就是添加地名，sp包提供了coordinates函数，来计算地图的中心坐标：
tmp=coordinates（Changsha）
print（tmp）
##[,1][,2]
##2121113.228.32
##2134113.728.23
##2136112.828.29
##2149112.328.13
##2182113.028.17
tmp=as.data.frame（tmp）
tmp$names=iconv（Changsha$NAME99,from='GBK'）
print（tmp）
##V1V2names
##2121113.228.32长沙县
##2134113.728.23浏阳市
##2136112.828.29望城县
##2149112.328.13宁乡县
##2182113.028.17长沙市市辖区
csmap+geom_text（aes（x=V1,y=V2,label=names）,family="GB1",data=tmp）如果需要支持更多字体，可
升级会员 