如何绘制一幅漂亮的路网地图? & 绘制阜阳市确诊病例的分布

如何绘制一幅漂亮的路网地图? & 绘制阜阳市确诊病例的分布

之前发现过一个很有意思的网站:Draw all roads in a city at once,可以绘制一个城市的路网地图,而且还绘制的很漂亮,就想用 R 画一下。另外我还从一则微信推文里面爬取了阜阳市新冠肺炎确诊病例的分布情况并将他们绘制在路网地图上。

首先导入我们需要的 R 包:

1
2
3
4
5
library(osmdata)
library(dodgr)
library(tidyverse)
library(sf)
library(lwgeom)

osmdata 包可以用来把 OpenStreetMap 数据导入 R 中。dodgr 包中的 dodgr_streetnet() 函数可以用来从 sf 表单中提取路网数据。lwgeom 包中的 st_make_valid() 可以用于矫正 sf 对象中 invalid 的 geometry。

首先我们读取一份阜阳的地图数据,这份地图可以从这里下载:GEOJSON 数据下载

1
2
3
# 阜阳市的地图
fy <- read_sf('阜阳市.json') %>%
st_make_valid()

然后获取阜阳市的路网数据:

1
2
3
4
5
6
7
8
9
# 获取安徽省阜阳市的经纬度范围
estimate_box <- osmdata::getbb("Fuyang, Anhui, China")
estimate_box
#> min max
#> x 118.99574 119.31574
#> y 36.54716 36.86716

# 提取这个经纬度范围内的街道路线
streets_raw <- dodgr_streetnet(estimate_box, expand = 1)

但是需要注意,这里的 streets_raw 里面包含的路网实际上是 estimate_box 这个经纬度网络包含的路网,我们可以从其中再提取和 fy 相交的部分,也就是阜阳市的路网了:

1
2
3
4
# 提取阜阳地区
streets <- streets_raw %>%
st_make_valid() %>%
st_intersection(fy)

我们再使用 ggplot2 把这个路网绘制出来:

1
2
3
4
5
6
7
8
9
10
11
12
ggplot() + 
geom_sf(data = streets, alpha = 0.4,
show.legend = F, size = 0.1,
color = "black") +
labs(caption = "中国 · 安徽 · 阜阳",
y = "") +
theme_ipsum(base_family = cnfont) +
worldtilegrid::theme_enhance_wtg() +
theme(plot.background = element_rect(fill = "#F7F2F2",
color = "#F7F2F2"),
panel.background = element_rect(fill = "#F7F2F2",
color = "#F7F2F2"))

可以把阜阳地图放在下面作为底图:

1
2
3
4
5
6
7
8
9
10
11
12
13
ggplot() + 
geom_sf(data = fy) +
geom_sf(data = streets, alpha = 0.4,
show.legend = F, size = 0.1,
color = "black") +
labs(caption = "中国 · 安徽 · 阜阳",
y = "") +
theme_ipsum(base_family = cnfont) +
worldtilegrid::theme_enhance_wtg() +
theme(plot.background = element_rect(fill = "#F7F2F2",
color = "#F7F2F2"),
panel.background = element_rect(fill = "#F7F2F2",
color = "#F7F2F2"))

当然可以对不同类型的路使用不同的颜色:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# devtools::install_github('gadenbuie/ggpomological')
library(ggpomological)
ggplot() +
geom_sf(data = fy) +
geom_sf(data = streets, alpha = 0.4, size = 0.1,
aes(color = highway), show.legend = FALSE) +
scale_color_manual(values = rep(ggpomological:::pomological_palette, 3)) +
labs(caption = "中国 · 安徽 · 阜阳",
y = "") +
theme_ipsum(base_family = cnfont) +
worldtilegrid::theme_enhance_wtg() +
theme(plot.background = element_rect(fill = "#F7F2F2",
color = "#F7F2F2"),
panel.background = element_rect(fill = "#F7F2F2",
color = "#F7F2F2"))

我们再回顾一下之前关于地区分组聚合的知识,计算阜阳市各个区县的道路总里程(我就不分道路的类型了)。

crs = 3055 对应 UTM Zone 27(EPSG:3055)坐标系,单位为“米”,更适合距离的计算,所以首先把 streets 的坐标系转换成 UTM Zone 27(EPSG:3055)坐标系,然后再计算道路长度:

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
30
31
32
33
34
streets3055 <- streets %>% 
st_transform(crs = 3055) %>%
mutate(length = st_length(streets))
streets_len = aggregate(x = dplyr::select(streets3055,
length),
by = st_transform(fy, 3055),
FUN = sum)
streets_len$name = fy$name
streets_len
library(units)
# 计算每个市的质心
streets_len <- cbind(streets_len, st_centroid(streets_len)) %>%
mutate(length = as.numeric(length) / 1000,
length = set_units(length, km))
# 绘图展示
ggplot() +
geom_sf(data = st_transform(streets_len, 4326),
aes(fill = as.numeric(length)),
size = 0.2, color = "black") +
geom_sf_text(data = st_transform(streets_len, 4326,
geometry = geometry.1),
aes(geometry = geometry,
label = name),
family = cnfont) +
scale_fill_gradient(low = "#f5c04a", high = "#c03728") +
labs(caption = "中国 · 安徽 · 阜阳",
title = "阜阳市各个区县的道路里程",
fill = "道路里程(km)") +
theme_ipsum(base_family = cnfont) +
worldtilegrid::theme_enhance_wtg() +
theme(plot.background = element_rect(fill = "#F7F2F2",
color = "#F7F2F2"),
panel.background = element_rect(fill = "#F7F2F2",
color = "#F7F2F2"))

新冠肺炎确诊病例在阜阳市的分布

今天早上看到阜阳发布公众号里面的推文枚举了阜阳所有病例的分布:https://mp.weixin.qq.com/s/K-wbnK6sXHW-E0vmNmyECQ ,下面爬取这个页面的数据并整理:

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
library(rvest)
# 从 https://mp.weixin.qq.com/s/K-wbnK6sXHW-E0vmNmyECQ 爬取整理所有阜阳的确诊病例分布
url <- "https://mp.weixin.qq.com/s/K-wbnK6sXHW-E0vmNmyECQ"
html <- read_html(url)
html %>%
html_nodes(xpath = '//*[@id="js_content"]/section') %>%
html_text() %>%
tibble::enframe(name = NULL) %>%
dplyr::filter(str_detect(value, "[\u4e00-\u9fa5]")) %>%
dplyr::filter(str_detect(value, "([1-9]*例")) %>%
slice(1:105) -> df
df

#> # A tibble: 105 x 2
#> value row
#> <chr> <dbl>
#> 1 一、颍州区(18例) 1
#> 2 1.文峰街道办事处奎星社区博物馆巷(1例,已出院) 2
#> 3 2.文峰街道办事处万达社区万达华府(1例,已出院) 3
#> 4 3.文峰街道办事处三里岗社区碧桂园(1例) 4
#> 5 4.鼓楼街道办事处双龙桥社区汇鑫小区(1例,已出院) 5
#> 6 5.鼓楼街道办事处白衣桥社区临泉路老师范学院院内(1例) 6
#> 7 6.颍西街道办事处西湖社区西城御景小区(2例,已出院1例) 7
#> 8 7.颍西街道办事处城郊社区徽州家园小区(1例) 8
#> 9 8.颍西街道办事处单桥社区(1例) 9
#> 10 9.清河街道办事处九里社区锦华第一郡小区(1例) 10
#> # … with 95 more rows

还需要再整理整理:

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
30
31
32
33
34
35
36
37
38
39
40
41
42
43
df$row <- as.numeric(rownames(df))

df %>%
mutate(
district = case_when(
between(row, 2, 15) ~ "颍州区",
between(row, 20, 26) ~ "颍东区",
between(row, 28, 33) ~ "颍泉区",
between(row, 35, 61) ~ "临泉县",
between(row, 63, 68) ~ "太和县",
between(row, 70, 90) ~ "阜南县",
between(row, 92, 100) ~ "颍上县",
between(row, 102, 103) ~ "界首市",
T ~ ""
)
) %>%
dplyr::filter(district != "") %>%
mutate(
count = str_match(value, "(([1-9]*)例")[,2],
value = str_remove(value, "[0-9]*\\."),
value = str_remove(value, "((.*))")
) %>%
mutate(
address = paste0("阜阳市", district, value)
) %>%
select(address, count, district) %>%
type_convert() -> fy_confirm
fy_confirm

#> # A tibble: 92 x 3
#> address count district
#> <chr> <dbl> <chr>
#> 1 阜阳市颍州区文峰街道办事处奎星社区博物馆巷 1 颍州区
#> 2 阜阳市颍州区文峰街道办事处万达社区万达华府 1 颍州区
#> 3 阜阳市颍州区文峰街道办事处三里岗社区碧桂园 1 颍州区
#> 4 阜阳市颍州区鼓楼街道办事处双龙桥社区汇鑫小区 1 颍州区
#> 5 阜阳市颍州区鼓楼街道办事处白衣桥社区临泉路老师范学院院内… 1 颍州区
#> 6 阜阳市颍州区颍西街道办事处西湖社区西城御景小区 2 颍州区
#> 7 阜阳市颍州区颍西街道办事处城郊社区徽州家园小区 1 颍州区
#> 8 阜阳市颍州区颍西街道办事处单桥社区 1 颍州区
#> 9 阜阳市颍州区清河街道办事处九里社区锦华第一郡小区… 1 颍州区
#> 10 阜阳市颍州区清河街道办事处九里社区九鼎华苑小区 1 颍州区
#> # … with 82 more rows

下面再使用百度地图接口解析每个患者地址的经纬度:

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
30
31
32
33
34
35
36
37
38
39
40
41
# devtools::install_github('badbye/baidumap')
# 使用百度地图接口解析每个患者地址的经纬度
# 设置百度密钥
options(baidumap.key = Sys.getenv("baidumap"))
library(progress)
library(baidumap)
fy_confirm$lon = NA
fy_confirm$lat = NA
pb <- progress_bar$new(total = nrow(fy_confirm))
for(i in 1:nrow(fy_confirm)){
pb$tick()
tempv <- getCoordinate(fy_confirm$address[i], formatted = T)
fy_confirm$lon[i] = tempv['longtitude']
fy_confirm$lat[i] = tempv['latitude']
}

# 把 fy_confirm 转成 sf 对象
fy_confirm %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326) -> fy_confirm_sf
fy_confirm_sf

#> Simple feature collection with 92 features and 3 fields
#> geometry type: POINT
#> dimension: XY
#> bbox: xmin: 114.9952 ymin: 32.44401 xmax: 116.3776 ymax: 33.45647
#> epsg (SRID): 4326
#> proj4string: +proj=longlat +datum=WGS84 +no_defs
#> # A tibble: 92 x 4
#> address count district geometry
#> <chr> <dbl> <chr> <POINT [°]>
#> 1 阜阳市颍州区文峰街道办事处奎星社区博物馆巷… 1 颍州区 (115.8262 32.89543)
#> 2 阜阳市颍州区文峰街道办事处万达社区万达华府… 1 颍州区 (115.8244 32.88011)
#> 3 阜阳市颍州区文峰街道办事处三里岗社区碧桂园… 1 颍州区 (115.8286 32.884)
#> 4 阜阳市颍州区鼓楼街道办事处双龙桥社区汇鑫小区… 1 颍州区 (115.8107 32.91112)
#> 5 阜阳市颍州区鼓楼街道办事处白衣桥社区临泉路老… 1 颍州区 (115.8112 32.8996)
#> 6 阜阳市颍州区颍西街道办事处西湖社区西城御景小… 2 颍州区 (115.8014 32.90791)
#> 7 阜阳市颍州区颍西街道办事处城郊社区徽州家园小… 1 颍州区 (115.7952 32.91192)
#> 8 阜阳市颍州区颍西街道办事处单桥社区… 1 颍州区 (115.7686 32.90729)
#> 9 阜阳市颍州区清河街道办事处九里社区锦华第一郡… 1 颍州区 (115.8172 32.87158)
#> 10 阜阳市颍州区清河街道办事处九里社区九鼎华苑小… 1 颍州区 (115.8218 32.86637)
#> # … with 82 more rows

再绘图展示:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
ggplot() + 
geom_sf(data = fy) +
geom_sf_text(data = fy_confirm_sf, aes(size = count),
label = "", family = "FontAwesome") +
geom_sf(data = streets, alpha = 0.4, size = 0.1,
aes(color = highway), show.legend = FALSE) +
scale_color_manual(values = rep(ggpomological:::pomological_palette, 3)) +
labs(caption = "阜阳市确诊病例分布",
y = "", x = "", size = "确诊\n数量") +
theme_ipsum(base_family = cnfont) +
worldtilegrid::theme_enhance_wtg() +
theme(plot.background = element_rect(fill = "#F7F2F2",
color = "#F7F2F2"),
panel.background = element_rect(fill = "#F7F2F2",
color = "#F7F2F2"))

蛮严重的。

潍坊市路网地图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
wf <- read_sf('潍坊市.json') %>% 
st_make_valid()
# 获取潍坊市的经纬度范围
estimate_box <- osmdata::getbb("Weifang, Shandong, China")
# 提取这个经纬度范围内的街道路线
streets_raw <- dodgr_streetnet(estimate_box, expand = 1)

# 提取潍坊地区
streets <- streets_raw %>%
st_make_valid() %>%
st_intersection(wf)

ggplot() +
geom_sf(data = streets, alpha = 0.4,
show.legend = F, size = 0.1,
color = "black") +
labs(caption = "中国 · 山东 · 潍坊",
y = "") +
theme_ipsum(base_family = cnfont) +
worldtilegrid::theme_enhance_wtg() +
theme(plot.background = element_rect(fill = "#F7F2F2",
color = "#F7F2F2"),
panel.background = element_rect(fill = "#F7F2F2",
color = "#F7F2F2"))

知识星球附件链接:https://t.zsxq.com/37M33FM

# R

评论

Your browser is out-of-date!

Update your browser to view this website correctly. Update my browser now

×