1 准备工作

1.1 数据加载

演示的数据放置在data文件夹下。

set.seed(123)
library(rgdal)
library(dplyr)
options(stringsAsFactors = F)
# bjcir <- read.csv("data/BJcir.csv")
bjcir<-matrix(data = runif(n = 7*31,min=1,max=20), nrow = 7, ncol = 31)
rownames(bjcir)<-c('114','101','106','109','228','112','229')
bjcir<-data.frame(bjcir)
bjcir$id<-rownames(bjcir)
data.lyr<-ogrListLayers(data.dsn<-'data/beijPM25.shp')[1]
beijingmap<-readOGR(dsn = data.dsn,layer = data.lyr)
## OGR data source with driver: ESRI Shapefile 
## Source: "data/beijPM25.shp", layer: "beijPM25"
## with 18 features
## It has 12 fields

行政区边界数据的空间多边形如下图所示(为了能显示在web地图上，需要转换到wgs84坐标系上)。鼠标点击任一行政区多边形会弹出该行政区划名称的提示。

library(leaflet)
popup <- paste0("<strong>Name: </strong>", 
                        beijingmap$NAME)
leaflet()%>%
    addPolygons(data=beijingmap%>%
                    spTransform('+init=epsg:4326'), popup=popup)%>%
    addTiles(urlTemplate = 'http://t0.tianditu.cn/DataServer?T=vec_w&X={x}&Y={y}&L={z}')%>%
    addMiniMap()

行政区边界的属性表如下所示(前两列)，

DT::datatable(beijingmap@data[,c(1,2)],rownames = F, class = 'cell-border compact stripe')

2 基本制图

2.1 地图主题

ggplot2提供了及其丰富的自定义出图方式，这里我们先定义一个用于出图的主题，方便后面使用。

theme_map<-function(base_size = 12, base_family = "Helvetica") {
    theme_bw(base_size = base_size, base_family = base_family) %+replace%
    theme(
        axis.line=element_blank(),
        # axis.text.x=element_blank(),axis.text.y=element_blank(),
        axis.ticks=element_blank(),
        # axis.title.x=element_blank(),axis.title.y=element_blank(),
        legend.position = "right",
        legend.text=element_text(size=base_size),
        legend.background =element_blank(),
        panel.background=element_blank(),
        panel.grid.major=element_blank(),panel.grid.minor=element_blank(),
        plot.background=element_blank())
}

2.2 数据处理

由于ggplot2在使用多边形会自动调用fortify清理结构，但fortify结果中只会包含分区信息

library(ggplot2)
bj.df<-fortify(beijingmap,region = 'object') # 指定分区
DT::datatable(bj.df,rownames = F, class = 'cell-border compact stripe')

这时我们需要在其结果上join属性表，而不是先join到多边形再出图

bjmap.df <- left_join(bj.df, bjcir,by='id')
DT::datatable(bj.df,rownames = F, class = 'cell-border compact stripe')

2.3 出图

准备工作完成，下面就是最后的出图步骤。因为我们使用fill的离散符号化，所以颜色的指定也需要用配套的函数指定色带。

ggplot()+
  geom_polygon(data = bjmap.df,
               aes(x = long, y = lat,group = group, fill = X1),
               color = "white", size = 0.25) +theme_map()+
    scale_fill_gradient(low = "yellow", high = "red")

3 动态地图

3.1 安装模块

首先安装gganimate包。由于这个包的代码发布在github上。 有些网络环境下，需要设置代理服务器才能顺畅的下载github上的代码。

library(httr)
set_config(use_proxy(url="127.0.0.1",port=8087))
set_config( config( ssl_verifypeer = 0L ) )

我们用devtools的功能直接从github上安装,完成后加载。

devtools::install_github("dgrtwo/gganimate")
library(gganimate)

3.2 数据变形

前面的日尺度数据格式是宽格式，

head(bjcir)
##            X1        X2        X3       X4        X5        X6       X7
## 114  6.463973 17.955962  2.955569 14.16326  6.494035 10.078123 8.860762
## 101 15.977798 11.477265 18.096674 13.16963  3.795159 15.410731 8.008064
## 106  8.770562  9.675680  5.675667 19.89113 19.297460  5.111751 3.896450
## 109 17.777331 19.179834  1.799131 13.45841 18.143682  7.045439 3.637315
## 228 18.868878  9.613349  7.230494 14.46208 14.123400  5.400890 5.427648
## 112  1.865573 13.873842 19.135569 11.33725 16.113881  3.713200 9.853287
##            X8        X9       X10       X11       X12       X13       X14
## 114 17.298727  3.423101  6.213289 15.335028 12.642649  2.954428 13.408937
## 101  1.870792 15.312849 16.478161 12.955201  7.684160  9.262962  7.526813
## 106  9.401801 18.005862  9.521810 14.493466  3.111573 19.714183 13.478404
## 109 16.179572  8.114793 16.391223  1.011871  5.628770 17.967971  7.087092
## 228  3.316086 13.637189 16.435401 10.031015 13.693056 17.842912  4.566131
## 112 11.658012  2.801973 16.092504  5.182259  8.935289  4.326000 15.863592
##           X15       X16       X17       X18       X19      X20       X21
## 114  9.868802 17.916654  2.153691  8.685695  3.929844 15.00448  6.922342
## 101 10.718604 18.374326 19.006812 13.309976  2.729836 10.90158  8.780024
## 106 12.399790 12.565965 14.691329  7.076592  3.696231 13.53693  1.198875
## 109  7.323647  8.803106  3.703592  6.846680 14.110135 16.61430  4.493141
## 228 10.283648  3.794799 11.436408  5.175585 12.765873 15.93935 17.011857
## 112 19.135003 18.770696 19.127734  8.020288 17.936488 19.61662  5.392074
##           X22       X23       X24       X25       X26       X27       X28
## 114  2.457132  3.110833  7.724187  8.842176 12.742340 14.457516 18.135613
## 101  5.668750  8.409894 13.349718  6.048047  8.072523  6.035338  6.209166
## 106 14.910569 11.866771  8.119565 12.969488 11.066878 12.292521  7.108172
## 109 17.101610  5.120962  7.753462  4.492741 17.618965 10.144506 19.727177
## 228 10.453018  9.450592 11.140071 17.409238 12.053252  6.035622 12.779873
## 112  8.370272  5.141823 15.066353 15.184792 16.955588 11.727218 18.808968
##           X29       X30       X31  id
## 114  8.729819 10.785565 10.158810 114
## 101 13.525376  8.648894  5.806334 101
## 106  3.894586 17.724684  5.108841 106
## 109 11.884474  7.917745 13.813151 109
## 228  5.535795  6.476546  1.905609 228
## 112 19.284820  4.242259 14.316209 112

而ggnanimate读取的是长格式数据。用reshape2转换格式神器(友情提示下要和reshape的函数区分开)看一下长格式，就明白宽和长的区别。

library(reshape2)
bjcir.long<-melt(bjcir,id.vars = 'id',variable.name = "day")
head(bjcir.long)
##    id day     value
## 1 114  X1  6.463973
## 2 101  X1 15.977798
## 3 106  X1  8.770562
## 4 109  X1 17.777331
## 5 228  X1 18.868878
## 6 112  X1  1.865573

显而易见，长格式就是把每一天的数值都放在一列，然后多出一个时间列以区分。在宽格式里面，不同时间的属性则是分别放在不同列里面。

3.3 渲染动态地图

bjmap.df.time <- left_join(bj.df, bjcir.long,by='id')
library(stringr)
bjmap.df.time$day<-  bjmap.df.time$day%>%
                            str_replace('X','')%>% as.numeric()

用gganimate制作动态图的最大优点，就是能最大限度的保持原有ggplot2语法结构。只需要在原有的aes中指明gganimate所需的frame，调用gganimate函数会自动识别为时空动态地图。

p<-ggplot()+
  geom_polygon(data = bjmap.df.time,
               aes(x = long, y = lat,group = group, fill = value,frame = day),
               color = "white", size = 0.25) +theme_map()+
    scale_fill_gradient(low = "yellow", high = "red") 
gganimate(p)

unnamed-chunk-14