注意
如果浏览器提示本文加载了不安全的脚本，请点允许。
今次主题比较简单。上个话题留了点冷饭，看起来还没馊，咱敲个鸡蛋炒个蛋炒饭。
炒个什么蛋炒饭呢？——动态图(dynamic charts)。这也是应留言要求额外发的番外。基本和分析关系不大，纯粹是可视化范畴。
动态图
通常我们看到的都是静态图，最常见的是.jpg、.png这类位图，逼格高一点的会用到.svg矢量图。但它们都是死图，所有图形元素都不会动。某些情况下，我们不仅要把统计结果映射到特定的视觉通道，还希望表现其历时性，或者允许用户自己进行挖掘。这就需要让图形部件动起来。
R有个名包，叫animation，可以用它压制.gif，用在社交媒体效果足够醒目。它的基本思路就是拿出一维来映射时间，基于时间点对数据切片、统计、制图，最后把静态图们合成一个动画。
- 还有一种方法是交互图(interactive charts)。
把统计数据绑定到JavaScript控件上，定义好交互方法，用户即可在网页上通过控件操作来调整视觉呈现（切片、缩放、改变类型等）。RStudio发过一个工具框架包htmlwidgets，可以很方便地把已有的JavaScript可视化库移植到R。我们今天就要用到其中的两个：ECharts2和leaflet。
如果再进一步，就是数据交互面板了。R有shiny及其系列衍生品，比如flexboard。想象一下作战室交互图仪表盘面板，几行命令就做出来了。简直酷炫。但是这需要部署在shiny服务器上。
动画
还是用前次的数据集，再画一下各朝的进士地理热图。但这次变个花样，我们先把北宋、明、清都按前、中、后分一下期。
- 用cut函数把连续变量切分成分段因子；
- 用sprintf格式化因子标签，否则会自动套用前开后闭科学计数法区间样式，表示年代比较怪；
- 各朝的分期不见得都对，我就随便分了下：
- 北宋(960-1127)：以1021、1085分期，大体对应真宗/仁宗、神宗/哲宗；
- 明朝(1668-1644)：以1434、1572分期，大体对应宣宗(宣德)/英宗、穆宗(隆庆)/神宗(万历)；
- 清朝(1644-1911)：以1735、1850分期，大体对应世宗(雍正)/高宗(乾隆)、宣宗(道光)/文宗(咸丰)。
1
2
3
4
5
6
7
8
9
10
11
12
13
|
library(animation)
nsong.js$period <- cut(
nsong.js$EntryYear, c(960, 1021, 1085, 1127), include.lowest=TRUE,
labels=paste(sprintf("%4d", cutyears[1:(length(cutyears)-1)]),
sprintf("%4d", cutyears[2:length(cutyears)]), sep="-"))
ming.js$period <- cut(
ming.js$EntryYear, c(1368, 1434, 1572, 1644), include.lowest=TRUE,
labels=paste(sprintf("%4d", cutyears[1:(length(cutyears)-1)]),
sprintf("%4d", cutyears[2:length(cutyears)]), sep="-"))
qing.js$period <- cut(
qing.js$EntryYear, c(1644, 1735, 1850, 1911), include.lowest=TRUE,
labels=paste(sprintf("%4d", cutyears[1:(length(cutyears)-1)]),
sprintf("%4d", cutyears[2:length(cutyears)]), sep="-"))
|
各自切片统计，再做成.gif动画。这里设置动画间隔为2秒，每一帧大小都是640×480像素。
1
2
3
4
5
6
7
8
9
|
oopt <- ani.options(interval=2, ani.width=640, ani.height=480)
saveGIF({
for (i in levels(nsong.js$period)){
print(p.song + stat_density_2d(aes(x_coord, y_coord, fill=..level..),
data=nsong.js[nsong.js$period==i,], geom="polygon", alpha=0.5) +
scale_fill_gradient(low="cyan", high="darkblue")+
ggtitle(paste0("北宋进士来源地 (", i, ")")))
}
}, "song.gif")
|
图 | 北宋进士来源地
很明显，前期京畿还很有优势，仁宗开始江南大盛，到神宗以后就是闽人大盛。热力分布持续东南移。
1
2
3
4
5
6
7
8
|
saveGIF({
for (i in levels(ming.js$period)){
print(p.ming + stat_density_2d(aes(x_coord, y_coord, fill=..level..),
data=ming.js[ming.js$period==i,], geom="polygon", alpha=0.5) +
scale_fill_gradient(low="cyan", high="darkblue")+
ggtitle(paste0("明朝进士来源地 (", i, ")")))
}
}, "ming.gif")
|
图 | 明代进士来源地
明和北宋热力分布是反向变动的，神宗以后进士来源地更弥散了。赣闽次第没落，而京畿有了起色。
1
2
3
4
5
6
7
8
9
|
saveGIF({
for (i in levels(qing.js$period)){
print(p.qing + stat_density_2d(aes(x_coord, y_coord, fill=..level..),
data=qing.js[qing.js$period==i,], geom="polygon", alpha=0.5) +
scale_fill_gradient(low="cyan", high="darkblue")+
ggtitle(paste0("清朝进士来源地 (", i, ")")))
}
}, "qing.gif")
ani.options(oopt)
|
图 | 清代进士来源地
清和明类似，中前期江南独大，鸦片战争以后分布逐渐弥散。
和静态图相比，引入历时维度后，信息量明显就更大了。
交互图
个体散点图
在哈佛人文地理平台上，进士数据集本身就以散点和热力图形式展示。这有一个问题：同一个来源地的坐标是完全一样的，所以这些个体在空间上都叠在一起，无法区分。在用leaflet再现散点分布之前，我们用随机数把他们打散开。
1
2
3
4
5
6
|
nsong.js$long <- nsong.js$x_coord + rnorm(nrow(nsong.js))/100
nsong.js$lat <- nsong.js$y_coord + rnorm(nrow(nsong.js))/100
ming.js$long <- ming.js$x_coord + rnorm(nrow(ming.js))/100
ming.js$lat <- ming.js$y_coord + rnorm(nrow(ming.js))/100
qing.js$long <- qing.js$x_coord + rnorm(nrow(qing.js))/100
qing.js$lat <- qing.js$y_coord + rnorm(nrow(qing.js))/100
|
这样一来，低级缩放时这些点还是叠在一起的，只要持续放大，就能散开。这当然很不精确，但是我感觉能忍。
构造一个工作函数make_leaflet，把参数传进去，就能画个leaflet交互图出来。leaflet的语法很简单，就是个leaflet js库的wrapper(咋翻译？包裹函数？)，参数列表几乎都是从API文档继承来的。基本的使用顺序是初始化->添加地图图层->添加主题图层->添加交互控件。我没多玩花活，加了个label和popup。
leaflet是基于svg矢量图的，源代码里就包含很多矢量代码。加上这里数据点多，文件尺寸就都很大。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
|
library(leaflet)
make_leaflet <- function(refMap, dyn, bgColor="red", dataset, cutyears){
dataset$period <- cut(
dataset$EntryYear, cutyears, include.lowest=TRUE,
labels=paste(sprintf("%4d", cutyears[1:(length(cutyears)-1)]),
sprintf("%4d", cutyears[2:length(cutyears)]), sep="-"))
dataset <- dataset[order(dataset$period),]
g <- leaflet() %>% addProviderTiles(providers$OpenStreetMap)
for (i in unique(refMap$id)){
g <- g %>%
addPolygons(~long, ~lat, data=refMap[refMap$id==i,],
label=dyn, labelOptions=labelOptions(textsize="20px"),
weight=1, fillColor=bgColor, color=bgColor)
}
pal <- colorFactor(
viridis::plasma(2*nlevels(dataset$period)),
levels=levels(dataset$period), ordered=TRUE)
g <- g %>% addCircleMarkers(
~long, ~lat, data=dataset, radius=2,
label=paste(dataset$NameChn, dataset$Name),
popup=paste(dataset$EntryYear, "<br>", dataset$Prov,
dataset$AddrChn, dataset$AddrName, sep=" "),
color=~pal(period), group=~period)
g %>% addLayersControl(
overlayGroups=rev(levels(dataset$period)),
options=layersControlOptions(
collapse=FALSE)
)
}
|
1
|
make_leaflet(nsong.bou, "北宋", "red", nsong.js, c(960, 1021, 1085, 1127))
|
点开查看源文件
1
|
make_leaflet(ming.bou, "明朝", "red", ming.js, c(1368, 1434, 1572, 1644))
|
点开查看源文件
1
|
make_leaflet(qing.bou, "清朝", "black", qing.js, c(1644, 1735, 1850, 1911))
|
点开查看源文件
古今地名一致性
前次分析，直接把所有县级地名拿出来跑频数，第一名是莆田。这当然很不精准——地名一直会变嘛。像南京这种历史上频繁改名的地方，就会吃很大亏。最好统一起来。找地名的历史沿革数据库不是很容易，干脆统一到今名如何？（请都说“好”）
统一到今名有三个好处：
- 能按照今天的区划统计哪个省市出了更多进士（其实谁稀罕这个）；
- 便于我操作（-口_口-：仄寺坠重要的）；
- 能引出一个知识点：sp::over()。
既然我们有每个省/市级空间多边形数据，又有每个进士的位置坐标，恰好又都是同一套大地坐标系测量值，历史沿革问题不就转化成了一个纯粹的GIS转换问题吗？只要用算法评价每个点落在哪个多边形内，它自然就属于这个行政区划。
这个算法有现成的，就是sp包的over函数。
构造一个which_polygon函数，把一个点轮流放进大陆和台湾的多边形数据集里，看它到底落在哪个多边形，输出这个多边形的元数据。注意：数据都得是sp包的工作类型，比如SpatialPoints、SpatialPolygons之类。
1
2
3
4
5
6
7
8
9
10
11
|
which_polygon <- function(point, cn.cities, tw.cities){
library(sp)
out <- over(SpatialPoints(
matrix(point, nrow=1),
proj4string = cn.cities@proj4string), cn.cities)
if (is.na(out$ID_0))
out <- over(SpatialPoints(
matrix(point, nrow=1),
proj4string = tw.cities@proj4string), tw.cities)
return(out)
}
|
每个进士个体都丢进去算太耗时，我们先统计下频数。
1
2
3
4
|
library(data.table)
nsong.js.stat <- dcast(nsong.js, x_coord+y_coord~., length)
ming.js.stat <- dcast(ming.js, x_coord+y_coord~., length)
qing.js.stat <- dcast(qing.js, x_coord+y_coord~., length)
|
再次动用并行计算parallel。因为要按行遍历矩阵，所以要用apply的并行版本parApply。获得结果是个list，需要再do.call调用bind_rows，合并成一个数据框。当然还有种办法是直接改写which_polygon，把进士位置整理成SpatialPointsDataFrame，一次性丢进去，效率上可能更高。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
|
library(parallel)
cl <- makeCluster(getOption("cl.cores", 2))
nsong.js.belong <- parApply(cl, nsong.js.stat[,c("x_coord", "y_coord")],
1, which_polygon, cn.cities, tw.cities)
ming.js.belong <- parApply(cl, ming.js.stat[,c("x_coord", "y_coord")],
1, which_polygon, cn.cities, tw.cities)
qing.js.belong <- parApply(cl, qing.js.stat[,c("x_coord", "y_coord")],
1, which_polygon, cn.cities, tw.cities)
stopCluster(cl)
library(dplyr)
nsong.js.belong <- do.call('bind_rows', nsong.js.belong)
ming.js.belong <- do.call('bind_rows', ming.js.belong)
qing.js.belong <- do.call('bind_rows', qing.js.belong)
nsong.js.belong$num <- nsong.js.stat$`.`
ming.js.belong$num <- ming.js.stat$`.`
qing.js.belong$num <- qing.js.stat$`.`
|
进士数排名
把三个朝代的数据集做一下深加工，合并到一起（DYNASTY列标识朝代属性）。
1
2
3
4
5
6
7
8
9
|
library(data.table)
nsong.js.belong <- dcast(nsong.js.belong, NAME_1+NAME_2+NL_NAME_2~., sum, value.var="num")
nsong.js.belong$DYNASTY <- "北宋"
ming.js.belong <- dcast(ming.js.belong, NAME_1+NAME_2+NL_NAME_2~., sum, value.var="num")
ming.js.belong$DYNASTY <- "明朝"
qing.js.belong <- dcast(qing.js.belong, NAME_1+NAME_2+NL_NAME_2~., sum, value.var="num")
qing.js.belong$DYNASTY <- "清朝"
library(dplyr)
js.belong <- do.call("bind_rows", list(nsong.js.belong, ming.js.belong, qing.js.belong))
|
再清洗一下数据，把英文地名里的无关符号删掉（“内江”出现了两种写法），再把0都转为NA（可视化时不要涂色）。
1
2
|
js.belong$NAME_2 <- stringr::str_replace(js.belong$NAME_2, "^([^]]+)\\]+$", "\\1")
js.belong$.[js.belong$.==0] <- NA
|
排名来了。
先是市级排名
1
2
3
|
js.order <- dcast(js.belong, NAME_1+NAME_2+NL_NAME_2~DYNASTY, sum,
value.var=".", margins="DYNASTY")
knitr::kable(js.order[order(js.order$`(all)`, decreasing=TRUE),])
|
地名归集后，第一名是苏州，然后是福州，第三才是莆田，第四是杭州。前四名实际难分伯仲。
| 序号 |
NAME_1 |
NAME_2 |
NL_NAME_2 |
北宋 |
明朝 |
清朝 |
(all) |
| 1 |
Jiangsu |
Suzhou |
苏州市 |
157 |
60 |
179 |
396 |
| 2 |
Fujian |
Fuzhou |
福州市 |
273 |
36 |
83 |
392 |
| 3 |
Fujian |
Putian |
莆田市 |
264 |
109 |
8 |
381 |
| 4 |
Zhejiang |
Hangzhou |
杭州市 |
115 |
38 |
221 |
374 |
| 5 |
Fujian |
Nanping |
南平市 |
224 |
15 |
7 |
246 |
| 6 |
Jiangxi |
Ji’an |
吉安市 |
128 |
97 |
11 |
236 |
| 7 |
Zhejiang |
Ningbo |
宁波市 |
94 |
97 |
39 |
230 |
| 8 |
Jiangsu |
Changzhou |
常州市 |
97 |
20 |
97 |
214 |
| 9 |
Zhejiang |
Jiaxing |
嘉兴市 |
19 |
68 |
117 |
204 |
| 10 |
Jiangxi |
Fuzhou |
抚州市 |
120 |
27 |
52 |
199 |
| 11 |
Zhejiang |
Huzhou |
湖州市 |
86 |
17 |
78 |
181 |
| 12 |
Zhejiang |
Shaoxing |
绍兴市 |
75 |
48 |
51 |
174 |
| 13 |
Fujian |
Quanzhou |
泉州市 |
110 |
28 |
31 |
169 |
| 14 |
Jiangsu |
Wuxi |
无锡市 |
73 |
34 |
55 |
162 |
| 15 |
Jiangxi |
Shangrao |
上饶市 |
93 |
26 |
28 |
147 |
给张图，同样是基于leaflet的。注意：文件非常很大，极费流量。
点开查看源文件
然后是省级排名
1
2
|
js.order.prov <- dcast(js.belong, NAME_1~DYNASTY, sum, value.var=".", margins="DYNASTY")
knitr::kable(js.order.prov[order(js.order.prov$`(all)`, decreasing=TRUE),])
|
| 序号 |
NAME_1 |
北宋 |
明朝 |
清朝 |
(all) |
| 1 |
Zhejiang |
636 |
337 |
530 |
1503 |
| 2 |
Fujian |
950 |
202 |
147 |
1299 |
| 3 |
Jiangsu |
476 |
179 |
536 |
1191 |
| 4 |
Jiangxi |
559 |
233 |
232 |
1024 |
| 5 |
Henan |
276 |
97 |
185 |
558 |
| 6 |
Shandong |
99 |
88 |
311 |
498 |
| 7 |
Anhui |
142 |
68 |
197 |
407 |
| 8 |
Sichuan |
269 |
30 |
65 |
364 |
| 9 |
Hebei |
27 |
86 |
200 |
313 |
| 10 |
Hunan |
79 |
27 |
152 |
258 |
| 11 |
Hubei |
31 |
52 |
133 |
216 |
| 12 |
Shanxi |
24 |
63 |
125 |
212 |
| 13 |
Guangdong |
40 |
54 |
105 |
199 |
| 14 |
Shaanxi |
18 |
49 |
89 |
156 |
| 15 |
Shanghai |
26 |
47 |
68 |
141 |
| 16 |
Beijing |
2 |
12 |
115 |
129 |
| 17 |
Yunnan |
0 |
1 |
69 |
70 |
| 18 |
Guangxi |
2 |
10 |
56 |
68 |
| 19 |
Tianjin |
0 |
2 |
61 |
63 |
| 20 |
Guizhou |
0 |
4 |
51 |
55 |
| 21 |
Chongqing |
8 |
7 |
24 |
39 |
| 22 |
Gansu |
6 |
6 |
19 |
31 |
| 23 |
Hainan |
0 |
9 |
3 |
12 |
| 24 |
Liaoning |
0 |
0 |
12 |
12 |
| 25 |
Jilin |
0 |
0 |
8 |
8 |
| 26 |
Ningxia Hui |
0 |
0 |
7 |
7 |
| 27 |
Taiwan |
0 |
0 |
2 |
2 |
同样给张图。这次是基于ECharts2的。由于这个包是我写的，所以我要额外安利一下http://madlogos.github.io/recharts。
1
2
3
4
5
6
7
8
9
|
library(recharts)
dict <- geoNameMap[geoNameMap$FKEY==31, c("EN", "CN")]
js.order.prov <- merge(js.order.prov, dict, by.x="NAME_1", by.y="EN", all.x=TRUE)
js.order.ec <- melt(js.order.prov[,c(2:4,6)], id="CN")
js.order.ec$value[js.order.ec$value==0] <- NA
echartR(js.order.ec, CN, value, t=variable, type="map_china", subtype="average") %>%
setDataRange(splitNumber=0, color=c('darkblue','cyan')) %>%
setTitle("各省进士数", pos=11) %>%
setTimeline(autoPlay=TRUE) %>% setLegend(FALSE)
|
点开查看源文件
图动起来了。感觉棒呆。
自己试着去和它们交互吧。
[完]
扫码关注我的公众号
