画地图的Notebook

情况提要

2019年1月7日中午接到TeamLeader钟的电话，果然有大事发生 我们的建模论文有点点小问题，地图要重新画 本来以为可以ps截图解决，但TeamLeader要求比较高，只好重画 GEODA我又没用过，于是，捡起几年前的R画地图的东西，7点开工 好在TeamLeader非常能干找到了底图，于事加班到11点半，完成

需要的基础数据

要画地图，毫无疑问，首先需要的是底图， 因为我们需要着色，必须是多边形地图，一般文件名为bou2_4p.shp 全名的话，就是boundary_polygon_polygon.shp，当然还有几个配套的文件

另外需要的就是自已的数据文件，根据数据中的变量来叠加颜色、标注等内容

主要过程

1. 把shp底图中的数据分解成数据框
2. 找到底图数据框与自已数据的对应关系
3. 将两者合并
4. 画图，调整样式，保存

正式开始啦

0. 载入3个包，maptools处理地图，tidyverse整理数据，stringr处理文字

library(maptools)
library(tidyverse)
library(stringi)

1.1 读入底图数据

citysmap <- readShapePoly("e:/gitR/map/zhong17/state/STATE_poly.shp")

1.2 读入要用到的数据

cdata <- read_csv('cdata.csv', locale = locale(encoding='GB18030'))

1.3 初始化用的基础数据

stdnames <- c("高水平", "较高水平", "中等水平", "较低水平", "低水平", '无数据')
std <- c(1.5,0.8,0.5,0.3,0)
cb_palette <- c('1'="#D32F2F", '2'="#CC79A7", '3'="#43A047", '4'="#76D7C4", '5'="#0072B2", '9'="#CFD8DC")
years=c(2009, 2012, 2014, 2016)

2.1 赋发展程度的函数

ranklevel <- function(value){
  return(min(which((c(1.5,0.8,0.5,0.3,0)-value)<=0)))
}

2.2 数据简单处理， 提取需要变量，赋发展程度，并检验

csdata <- select(cdata, CITYCN, CITYEN, datayear, cityrank)
csdata$level <- as.factor( sapply(csdata$cityrank, ranklevel))
table(csdata$level)

   1    2    3    4    5 
  36  178  453 1049  244 

3.1 把底图数据用fortify函数提取数据框， 提取后数据结构如下

mapnp <- fortify(citysmap)

Regions defined for each Polygons

head(mapnp)

提取数据框后的底图，只有经纬度数据和分组数据， 标明地域的只有id，但id只是没有实际意义的序列

3.2 所以必须取出底图的地域信息

mdata <- citysmap@data
head(mdata)

可以发现，这个基本信息比较全，有地区的中文名称，行政中心位置面积等 另外可以发现，底图数据框的id就是这个信息数据框的行号 这样，经过简单处理，两部分数据就可以拼合了 不过，鉴于底图数据框行数很大，有523987行， 最好一次只拼合一年的数据，以免内存溢出

mdata <- mutate(mdata, id=rownames(mdata), CITYCN=as.character(NAME99))

4.1 拼合数据并画图， 用中文名称将自有数据与底图信息数据拼合， 用id将底图信息与底图数据框拼合， 这样就能用ggplot2画图了 调整好配色，数据标签位置，标图完整后就可以保存了

plotyear <- function(year){
  mdata.y <- left_join(mdata, filter(csdata, datayear==year), by='CITYCN')
  mapnp.y <- left_join(mapnp, mdata.y, by='id')
  
  pm <- ggplot(data = mapnp.y) 
  pm <- pm +   geom_polygon(aes(x = long, y = lat, fill=level, group = group), color = "black")
  pm <- pm + labs( x = "经度", y = "纬度",title =paste(year,'年'))
  pm <- pm + theme(legend.position=c(0.1,0.15))
  pm <- pm + scale_fill_manual(values=cb_palette, name = "高质量发展水平", labels = stdnames)
  print(pm)
  ggsave(paste(year, '.png'), pm, width = 10, height = 8, dpi = 150)
}

5.1 循环执行需要作图的年份，即可得到多年的图了

for(year in years) plotyear(year)

后记

TeamLeader的标准果然不同凡响，所以她是TeamLeader， 8号一早，我又画了2个小时将图P成可以用的样子

地图杂谈

其实我们画中国的图要烦死了，因为R并没有自带中国内部的底图， 如果要画美国的地图简单死了，因为已有R包内置了美国的底图， 比如画个美国各州的地图，只要读入相应数据画下就OK 此部分参考自Eric C. Anderson 的 Making Maps with Rhttp://eriqande.github.io/rep-res-web/lectures/making-maps-with-R.html

states <- map_data("state")

Attaching package: 愼㸱愼㹥maps愼㸱愼㹦

The following object is masked from 愼㸱愼㹥package:purrr愼㸱愼㹦:

    map

ggplot(data = states) + 
  geom_polygon(aes(x = long, y = lat, fill = region, group = group), color = "white")+
  guides(fill=FALSE)

而且，更加变态的是，美国的地图数据可以分县（ps，美国没有地级行政区划）， 演示一个，把加州各县的图画一下，就是这么easy

counties <- map_data("county")
ca_county <- subset(counties, region == "california")
ggplot(data = ca_county) + 
  geom_polygon(aes(x = long, y = lat, fill = region, group = group), color = "white")+
  guides(fill=FALSE)

但是，很抱歉，如果想读中国的是没有的，如下：

china <- map_data("china")

Error: 'chinaMapEnv' is not an exported object from 'namespace:maps'

不过用这个数据还可以画世界地图倒是没问题

world <- map_data("world")
ggplot(data = world) + 
  geom_polygon(aes(x = long, y = lat, fill = region, group = group), color = "white")+
  guides(fill=FALSE)

现在是Wed Jan 09 13:10:10 2019，该就寝了，bye my305