第11章 地图
学习R花的精力就是代价，就像买个商业软件所支付的美元或欧元一样。— Felix Grant, November 2004
地图数据，经常用到的是地理信息系统GIS数据，很多人都用昂贵的 ArcGIS 软件来处理。其实，免费的R 配上强大的扩展包，也能够处理很多GIS问题，有时甚至更灵活。这里,我们举几个最简单的例子，来看看R是如何处理地图数据的。
11.1 绘制点阵地图
最简单的地图数据，莫过于用我们熟悉的数据框形式存储的经纬度坐标了。下面，我们用R来处理这样的点阵地图数据。
示例数据文件可以来我们的服务器下载。我们先用download.file()函数，将文件下载到本地电脑硬盘里。
download.file(url = "http://dapengde.com/r4rookies/us.csv",
destfile = "c:/r4r/us.csv")好了，已经下载到硬盘上了。请用Excel或记事本打开本地的这个文件，看看内容是什么。
这个数据框有3列，依次是经度、纬度、州名。州名包含了美国除夏威夷和阿拉斯加之外各州。
首先，我们读入数据，并看看总结报告：
us <- read.csv("c:/r4r/us.csv")
summary(us)第三列是州名，我们可以用因子水平函数来看看有几个州：
nlevels(us$state)## [1] 49这样的数据，我们只要像前面接触的普通数据一样，绘制以经度为x，纬度为y的散点图，那么就会画出一张美国点阵地图。
plot(us$lon, us$lat)用第3章学到的方法，我们可以给某个州指定个颜色，比如德克萨斯州：
plot(us$lon, us$lat)
points(us$lon[us$state == "Texas"],
us$lat[us$state == "Texas"], col="blue")类似地，可以给每个州涂上不同的彩虹色。
cols <- rainbow(nlevels(us$state))
plot(us$lon, us$lat, col = cols[us$state], pch = 20)可想而知，如果我们有各州的某个统计数据，例如人口密度、气温、选票支持率、PM2.5浓度等，以不同的颜色表示其大小，那么就都可以用类似的方法绘制在地图上。
利用某个城市的经纬度，就可以把城市在图上标出来。例如，我们用添加数据点的方法，把纽约(经纬度40.73, -74.02)用三角形标出来。我们还可以用添加文字的方法，给每个州标上州名，标注的位置是各州图形的中点：
lon.median <- tapply(us$lon, us$state, median)
lat.median <- tapply(us$lat, us$state, median)
text(labels=levels(us$state), x=lon.median, y=lat.median,
cex=0.5, col = 'White')
points(-74.02, 40.73, pch = 17, cex = 2)
图 11.1: R绘制点阵地图.
看得出来，这些都是我们前面学过的基本作图方法，没什么新鲜的。
11.2 绘制矢量地图
矢量地图是平时更常见的地图文件。shape格式的地图文件，是美国环境系统研究所公司（ESRI）开发的一种空间数据开放格式，是地理信息软件界的一个开放标准。
下面，我们用矢量地图文件，结合网上找到的各省空气质量数据，来绘制一张中国空气质量地图。
中国地图文件可以在网上下载34。为了方便，我们已经事先把找到的文件放到了我们自己的服务器上。所以，下面的代码，先从我们服务器下载一个名为cm.zip的压缩文件，然后用unzip()函数，将这个压缩包文件解压缩。
download.file(url = "http://dapengde.com/r4rookies/cm.zip",
destfile = "c:/r4r/cm.zip")
unzip(zipfile = "c:/r4r/cm.zip", exdir = "c:/r4r")得到的文件夹下应该有三个文件，即 bou2_4p.dbf，bou2_4p.shp，bou2_4p.shx。这就是我们的地图文件。读取这套文件，需要先安装和调用rgdal扩展包(Bivand, Keitt, and Rowlingson 2016)。
install.packages('rgdal') require(rgdal)
cm <- readOGR(dsn = "c:/r4r/cm", layer = "bou2_4p") ## OGR data source with driver: ESRI Shapefile
## Source: "c:/r4r/cm", layer: "bou2_4p"
## with 925 features
## It has 7 fields
## Integer64 fields read as strings: BOU2_4M_ BOU2_4M_IDreadOGR()函数用来读取地图文件。这个函数的自变量格式比较特别，dsn是数据源名称data source name的缩写，需要指定为地图文件所在文件夹路径，并且末尾不要斜线符号；layer需要指定为地图文件名， 但不要扩展名。
读入数据后，就可以直接画地图了：
plot(cm)下面，我们按照空气污染程度，给各省图上不同的颜色。先看看有哪些省。
summary(cm)
is.factor(cm$NAME)
levels(cm$NAME) # 省市名称。cm$NAME里存放的是各省名称，是个因子。我们只需把各省的空气质量以不同颜色或灰度表示，然后跟省名称的因子联系起来就可以了。
各省的空气质量，我们可以从中国空气质量在线监测分析平台35得到，这里以2016年12月30日14时的空气质量指数(AQI)为例，当时的AQI数据可以来我们的服务器下载：
## download.file(url =
## "http://dapengde.com/r4rookies/aqi.csv",
## destfile = "c:/r4r/aqi.csv")
aqi <- read.csv("c:/r4r/aqi.csv")
aqi## aqi province
## 1 77 福建省
## 2 90 甘肃省
## 3 190 北京市
## 4 157 安徽省
## 5 79 吉林省
## 6 88 江苏省
## 7 90 江西省
## 8 61 广东省
## 9 81 广西壮族自治区
## 10 72 贵州省
## 11 98 辽宁省
## 12 129 宁夏回族自治区
## 13 71 内蒙古自治区
## 14 54 青海省
## 15 163 新疆维吾尔自治区
## 16 112 山东省
## 17 172 山西省
## 18 112 四川省
## 19 160 陕西省
## 20 39 上海市
## 21 -1 台湾省
## 22 45 西藏自治区
## 23 -1 香港特别行政区
## 24 227 天津市
## 25 46 海南省
## 26 93 浙江省
## 27 45 云南省
## 28 244 河北省
## 29 226 河南省
## 30 80 重庆市
## 31 158 湖北省
## 32 119 湖南省
## 33 47 黑龙江省-1表示无数据。我们按照空气质量的分级标准，用cut()函数来计算各省空气属于哪一级（0 – 50优，50 –100 良，100 – 150轻污，150 – 200 中污，200 – 300 重污，300 – 500严重）：
aqstandard <- c(0, 50, 100, 150, 200, 300, 500, Inf)
aqilevel <- cut(aqi$aqi[match(levels(cm$NAME),aqi$province)],
aqstandard)我们用match()函数调整了一下顺序，以便跟每个省区一一对应。
根据级别的多少，来确定分几个灰度显示，然后把灰度设为各省名称的因子，就可以作图了：
mycol <- grey(seq(1,0,
length.out = nlevels(aqilevel)))[aqilevel]
col <- cm$NAME
levels(col) <- mycolplot(cm, col = as.character(factor(col)), axes = TRUE)最后，我们用plotrix扩展包(Lemon 2006)的color.legend()函数，给我们的地图添加个漂亮的图例：
install.packages('plotrix')library(plotrix)
legendn <- character((length(aqstandard) - 2) * 2 + 1)
legendn[seq(2, length(legendn), by = 2)] <-
aqstandard[2:(length(aqstandard)-1) ]
color.legend(
xl = 135, yb = 10, xr = 137, yt = 30,legend = legendn,
rect.col = grey(seq(1, 0, length.out = nlevels(aqilevel))),
align = "lb", gradient = "y") # 添加图例。
color.legend(
xl = 135, yb = 10, xr = 137, yt = 30,
legend = c('Excellent','Good','Lightly', 'Moderately',
'Heavily', 'Severely', 'OMG'),
rect.col = grey(seq(1, 0, length.out = nlevels(aqilevel))),
align = "rb", gradient = "y") # 添加图例。
图 11.2: R绘制矢量地图.
11.3 绘制交互地图
交互地图，就是像谷歌、百度、高德那种鼠标点击可以放大缩小的地图，看起来很炫，料想大约会更难画出来吧？其实不然，因为我们有现成的扩展包——leafletCN (Lang 2017)。如果事先安装好这个扩展包，并且按要求准备好数据，那么做出交互式地图只需一条指令，简单得让人不画张地图都不好意思。
我们先来安装和调用leafletCN扩展包，然后将上一节使用的数据框aqi稍微转化一下格式：
install.packages('leafletCN')
require(leafletCN)
aqi$aqi[aqi$aqi == -1] <- NA
pvs <- regionNames("china")
loc <- match(pvs, aqi$province)
aqi2 <- data.frame(name = pvs, value = aqi$aqi[loc])aqi2数据框就是作图要用的数据，含有两列：名为name的一列是全国省区的名称，取值来自leafletCN扩展包的regionNames()函数；名为value的一列是各省区的AQI数值，来自先前用用过的aqi数据框，我们用match()函数调整了一下顺序，以便跟name的每个省区一一对应。regionNames()的具体用法，请咨询F1小助理。
好了，万事俱备，只差一句话：
geojsonMap(dat = aqi2, mapName = "china",
popup = paste(aqi2$name, aqi2$value),
legendTitle = "AQI")
图 11.3: 交互地图示例.
geojsonMap()函数里，我们设置了popup参数，这样，鼠标放在地图上点击省区，就会弹出AQI数值。
怎么样，是不是很炫？
11.4 课外活动：绘制动画地图
在本章，我们绘制了一张中国大气质量地图，输入的数据是一组AQI数值。如果我们有多组数据，每个小时的数据做一张地图，有一年的数据，那么该怎么做？如果能做出来的话，连续播放，不就是动态显示我国空气质量的时空分布了吗？试试看。其中涉及的函数，我们都已经学过了。
提示：
将各省的AQI数据保存到一个数据框里，用循环语句每次读取其中的一组来作图。
图片保存为文件时，用字符串的处理方法为不同的文件取不同的名称。
可以用自定义函数的方法，将绘图过程写成一个函数，方便每次调用。
References
Bivand, Roger, Tim Keitt, and Barry Rowlingson. 2016. Rgdal: Bindings for the Geospatial Data Abstraction Library. https://CRAN.R-project.org/package=rgdal.
Lemon, J. 2006. “Plotrix: A Package in the Red Light District of R.” R-News 6 (4): 8–12.
Lang, Dawei. 2017. LeafletCN: An R Gallery for China and Other Geojson Choropleth Map in Leaflet. https://CRAN.R-project.org/package=leafletCN.
中国空气质量在线监测分析平台：https://www.aqistudy.cn/↩