第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)
R绘制点阵地图.

图 11.1: R绘制点阵地图.

看得出来,这些都是我们前面学过的基本作图方法,没什么新鲜的。

Example 11.1 在美国地图上标出三个你感兴趣的城市的位置。

Example 11.2 在美国地图上画出纬度线,间隔为5度。并找出北纬35度线穿过美国哪些州。

Example 11.3 重画美国地图,用不同颜色来区分时区。美国跨过几个时区?

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_ID

readOGR()函数用来读取地图文件。这个函数的自变量格式比较特别,dsn是数据源名称data source name的缩写,需要指定为地图文件所在文件夹路径,并且末尾不要斜线符号;layer需要指定为地图文件名, 但不要扩展名。

读入数据后,就可以直接画地图了:

plot(cm)
思考 11.1 矢量地图跟点阵地图相比,有什么优势?

下面,我们按照空气污染程度,给各省图上不同的颜色。先看看有哪些省。

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) <- mycol
plot(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")  # 添加图例。
R绘制矢量地图.

图 11.2: R绘制矢量地图.

Example 11.4 示例的地图用的是灰度,有的省份的大气污染级别从图例上不太容易区分,请将上图的颜色改为容易区分的彩色,重新绘制。

Example 11.5 从网上下载世界地图并画出来,给不同的国家地区涂上颜色。

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数值。

怎么样,是不是很炫?

Example 11.6 从网上下载全国各市的空气质量数据,并在交互地图上用不同颜色标示。

11.4 课外活动:绘制动画地图

在本章,我们绘制了一张中国大气质量地图,输入的数据是一组AQI数值。如果我们有多组数据,每个小时的数据做一张地图,有一年的数据,那么该怎么做?如果能做出来的话,连续播放,不就是动态显示我国空气质量的时空分布了吗?试试看。其中涉及的函数,我们都已经学过了。

提示:

  1. 将各省的AQI数据保存到一个数据框里,用循环语句每次读取其中的一组来作图。

  2. 图片保存为文件时,用字符串的处理方法为不同的文件取不同的名称。

  3. 可以用自定义函数的方法,将绘图过程写成一个函数,方便每次调用。

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.


  1. 矢量地图下载:http://www.diva-gis.org/gdata

  2. 中国空气质量在线监测分析平台:https://www.aqistudy.cn/