R for map(笔记1)

1、地图数据准备

首先，下载 中国省级行政划分的地图GIS数据，解压后文件夹内含三个地图文件：bou2_4p.dbf、bou2_4p.shp、bou2_4p.shx。也可以直接从http://www.gadm.org/country下载中国的地图数据包CHN_adm_shp，内含国家、省级、县区级的详细地图信息。地图数据基本可以分为点、线、面三种数据，在maptools包内分别有对应的函数来读取（readShapePoints、readShapeLines和readShapePoly函数）。

#加载maptools包，通过readShapePoly函数读入省级地图据，并画出最基础的中国地图：
library(maptools)
x <- readShapePoly(file.choose()) #读取bou2_4p.shp
plot(x)

#加载ggplot2包，用ggplot绘制，并使用polyconic投影方式，使地图显示正常：
library(ggplot2)
ggplot(x,aes(x=long,y=lat,group=group)) +
  geom_polygon(fill="white",colour="grey") +
  coord_map("polyconic")

2、GIS 地图的数据结构

#检测数据结构：
class(x)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
length(x) 
[1] 925
names(x)
[1] "AREA"       "PERIMETER"  "BOU2_4M_"   "BOU2_4M_ID" "ADCODE93"   "ADCODE99"  
[7] "NAME"      

此时，x中保存的是各个省/直辖市的多边形面图，数据类型是SpatialPolygonsDataFrame。x中有925条记录,每条记录中含有面积(AREA)、周长(PERIMETER)、各种编号、中文名(NAME)等字段。

其中中文名（NAME）字段是以GBK编码的,而我的工作环境是 utf-8,所以直接提取各个省/直辖市的名称得到的会是乱码，需要iconv函数转化一下才能正常显示。

#提取每个省/直辖市的多边形数目
table(iconv(x$NAME, from = "GBK"))
          安徽省           北京市           重庆市           福建省           甘肃省 
               1                1                1              168                1 
          广东省   广西壮族自治区           贵州省           海南省           河北省 
             154                6                2               79                9 
          河南省         黑龙江省           湖北省           湖南省           吉林省 
               1                1                1                1                1 
          江苏省           江西省           辽宁省     内蒙古自治区   宁夏回族自治区 
               5                1               94                1                1 
          青海省           山东省           山西省           陕西省           上海市 
               1               86                1                1               12 
          四川省           台湾省           天津市       西藏自治区   香港特别行政区 
               1               57                1                1               53 
新疆维吾尔自治区           云南省           浙江省 
               1                1              179 

3、作图

#提取出北京市的地图数据并plot 出北京市地图：
tmp <- iconv(x$NAME, from = "GBK")
grep("北京", tmp, value = TRUE)
[1] "北京市"
grep("北京", tmp)
[1] 8
x$ADCODE99[grep("北京", tmp)]
[1] 110000
Beijing <- x[x$ADCODE99 == 110000,]
plot(Beijing)

#进一步加工，使其显示正常：
ggplot(Beijing,aes(x=long,y=lat,group=group)) +
  geom_polygon(fill="white",colour="grey") +
  coord_map("polyconic")

按县区划分北京市:

#从CHN_adm3_shp读入中国县区级的地图信息：
z <- readShapePoly(file.choose()) 
library(ggplot2)
z1 <- fortify(z) #将SpatialPolygonsDataFrame转化为DataFrame
z2 <- z@data #读取行政信息
zs <- data.frame(z2,id=seq(0:2407)-1)
#合并形状数据和行政信息：
library(plyr)
china_map_data1 <- join(z1, zs, type = "full") #提示：Joining by: id

#取出北京的子集：
beijing1 <-subset(china_map_data1,NAME_1=="Beijing")
#分地区着色：
ggplot(beijing1, aes(x = long, y = lat, group = group,fill=NL_NAME_3))+
  geom_polygon( )+
  geom_path(colour = "grey40")+
  scale_fill_manual(values=colours(),guide=FALSE)+
  coord_map("polyconic")

#标注县区名称
region <- readShapePoly(file.choose()) #CHN_adm3_shp 读取经纬度失败
#取形状内的平均坐标:
midpos <- function(x) mean(range(x,na.rm=TRUE)) 
centres <- ddply(beijing1,.(NAME_3),colwise(midpos,.(long,lat)))
#按 longitude 渐变作图，并标注县区名称
ggplot(beijing1,aes(long,lat))+       
  geom_polygon(aes(group=group,fill=long),colour="black")+
  scale_fill_gradient(low="white",high="steelblue") +
  coord_map("polyconic") +
  geom_text(aes(label=NAME_3),data=centres)+
  #去除不需要的元素
  theme(
    panel.grid = element_blank(),
    panel.background = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    axis.title = element_blank()
  )