使用 ggplot2 + sf 绘制中国地图

使用 ggplot2 + sf 绘制中国地图

这是使用 R 语言绘制地图的第一篇教程,学习之前大家需要对 R 语言中的数据框操作有所了解。建议阅读 《R 数据科学》 一书。

使用 ggplot2 + sf 绘制地图可以使用 shp 数据,我提供的附件里有两个文件夹,一个是 chinamap,这个文件夹里面有很多数据,包括:

视频讲解

  1. 中国县界;
  2. 中国市界;
  3. 中国省界;
  4. 中国湖泊;
  5. 主要公路;
  6. 主要河流;
  7. 主要铁路;
  8. 县城驻地;
  9. 国界线;
  10. 地级城市驻地;
  11. 省会城市;
  12. 线状县界;
  13. 线状省界;
  14. 经纬网。

这个地图数据比较老,仅作为教程示例演示使用。

视频讲解(上篇)

链接:https://pan.baidu.com/s/1o9yUdBr61AUbDvzX01timA 密码:elny

准备工作

首先是导入所需的 R 包:

R
1
2
3
4
5
6
7
8
9
10
11
library(ggplot2)
library(sf)
library(tidyverse)
#> ── Attaching packages ──────────────────────────── tidyverse 1.3.0 ──
#> ✓ tibble 3.0.0 ✓ purrr 0.3.3
#> ✓ tidyr 1.0.2 ✓ stringr 1.4.0
#> ✓ readr 1.3.1 ✓ forcats 0.5.0
#> ── Conflicts ─────────────────────────────── tidyverse_conflicts() ──
#> x stats::filter() masks dplyr::filter()
#> x stats::lag() masks dplyr::lag()
library(hrbrthemes)

这四个包中,ggplot2 是用于绘图的,里面提供了 geom_sf() 函数用于绘制地图;sf 包是用于读取处理地理数据的,tidyverse 包可以理解为包管理器,调用 tidyverse 包的时候会同时调用一些常用的数据处理包并且检查包与包之间函数名称冲突的问题。

绘制省级填充地图

首先读入中国省界的 shp 数据:

R
1
mapdata <- read_sf("chinamap/中国省界.shp")

查看 mapdata 的类别:

R
1
2
class(mapdata)
#> [1] "sf" "tbl_df" "tbl" "data.frame"

可以看出我们得到的是个 sf 类数据,同时它也是 data.frame 类的文件,这非常重要,因为我们知道,R 里面最好处理的数据叫数据框。

如果你想深入学习 ggplot2 的绘图方法,可以学习:ggplot2: Elegant Graphics for Data Analysis

R
1
2
ggplot(mapdata) + 
geom_sf()

可以看到,没有九段线,然后再读入国界数据:

R
1
2
3
4
mapborder <- read_sf("chinamap/国界线.shp")
ggplot() +
geom_sf(data = mapdata, aes(geometry = geometry)) +
geom_sf(data = mapborder, aes(geometry = geometry))

注意到这里我使用了 aes(geometry = geometry)aes() 里面是设定映射的,这里设定的映射是把 mapdata/mapborder 数据集里的 geometry 变量映射给 geometry 参数。ggplot2 包里有很多图层,不同的图层有不同的功能,例如这里的 geom_sf 图层是用于绘制地图图层的,每种图层都有自己需要的一些参数,geometry 参数就是 geom_sf 图层的必备参数,大家也会注意到,我刚才第一次绘制地图的时候没有指定 aes(geometry = geometry),这是因为 geometry = geometry 就是默认设置,所以省略也可。

我现在手头上也没有各个省的数据,我们可以先随机生成一些数据:

R
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
44
45
46
mapdata
#> Simple feature collection with 34 features and 1 field
#> geometry type: MULTIPOLYGON
#> dimension: XY
#> bbox: xmin: -2578822 ymin: 2367106 xmax: 2092054 ymax: 6385320
#> epsg (SRID): NA
#> proj4string: +proj=lcc +lat_1=30 +lat_2=62 +lat_0=0 +lon_0=105 +x_0=0 +y_0=0 +ellps=krass +units=m +no_defs
#> # A tibble: 34 x 2
#> NAME geometry
#> <chr> <MULTIPOLYGON [m]>
#> 1 黑龙江 (((1329152 5619034, 1323328 5628245, 1317709 5635382, 1…
#> 2 新疆 (((-2189253 4611401, -2202922 4598047, -2217335 4598429…
#> 3 山西 (((761692.1 4443125, 760999.9 4441019, 758470.6 4437833…
#> 4 宁夏 (((-34477.05 4516814, -41105.13 4519021, -49016.8 45177…
#> 5 西藏 (((-2189253 4611401, -2187862 4611655, -2184182 4612719…
#> 6 山东 (((915805.7 4438425, 917551.3 4441336, 916454.1 4444794…
#> 7 河南 (((915805.7 4438425, 913870 4434178, 914650.2 4427691, …
#> 8 江苏 (((1261146 4381810, 1264549 4379112, 1262796 4375998, 1…
#> 9 安徽 (((1016688 4289115, 1018854 4289081, 1023048 4288111, 1…
#> 10 湖北 (((547929.6 4087822, 551532.2 4087997, 553822.7 4081664…
#> # … with 24 more rows

# 生成随机数据

mapdata$pop <- round(runif(34, 100, 1000))
mapdata
#> Simple feature collection with 34 features and 2 fields
#> geometry type: MULTIPOLYGON
#> dimension: XY
#> bbox: xmin: -2578822 ymin: 2367106 xmax: 2092054 ymax: 6385320
#> epsg (SRID): NA
#> proj4string: +proj=lcc +lat_1=30 +lat_2=62 +lat_0=0 +lon_0=105 +x_0=0 +y_0=0 +ellps=krass +units=m +no_defs
#> # A tibble: 34 x 3
#> NAME geometry pop
#> <chr> <MULTIPOLYGON [m]> <dbl>
#> 1 黑龙江 (((1329152 5619034, 1323328 5628245, 1317709 5635… 729
#> 2 新疆 (((-2189253 4611401, -2202922 4598047, -2217335 4… 900
#> 3 山西 (((761692.1 4443125, 760999.9 4441019, 758470.6 4… 980
#> 4 宁夏 (((-34477.05 4516814, -41105.13 4519021, -49016.8… 188
#> 5 西藏 (((-2189253 4611401, -2187862 4611655, -2184182 4… 459
#> 6 山东 (((915805.7 4438425, 917551.3 4441336, 916454.1 4… 786
#> 7 河南 (((915805.7 4438425, 913870 4434178, 914650.2 442… 938
#> 8 江苏 (((1261146 4381810, 1264549 4379112, 1262796 4375… 781
#> 9 安徽 (((1016688 4289115, 1018854 4289081, 1023048 4288… 483
#> 10 湖北 (((547929.6 4087822, 551532.2 4087997, 553822.7 4… 278
#> # … with 24 more rows

前面我们也说了,mapdata 是个数据框,数据框里面的变量可以使用 $ 符号进行选择和生成,从上面的代码中我们可以看到,本来 mapdata 数据集只有两个变量:NAME 和 geometry,运行 mapdata$pop <- round(runif(34, 100, 1000)) 之后该数据集里面多了个 pop 变量,而 pop 变量是什么呢,round(runif(34, 100, 1000))。runif() 函数生成了 34 个服从 100~1000 上的均匀分布的值(它们组合成了一个向量),round() 函数默认保留 0 位小数。

下面我们就可以使用这个 pop 变量进行地图填充了:

R
1
2
3
4
5
6
7
8
ggplot() + 
geom_sf(data = mapdata, aes(geometry = geometry,
fill = pop),
color = "white",
size = 0.05) +
geom_sf(data = mapborder, aes(geometry = geometry),
size = 0.05) +
scale_fill_viridis_c()

这里的代码和上面的代码相比,有三个地方发生了改动:1. fill = pop,注意这个是放在 aes() 里面的,表明我们是把 pop 变量映射给 fill。那么 pop 变量是个数值向量,如何映射给 fill(使用颜色填充)呢?这就是 scale_fill_viridis_c() 的作用了,这个函数用于指明我们的映射规则,大家看到的这个图就是使用这种映射规则绘图的结果:pop 的最小值是深紫色,最大值是明黄色。。。2. color = "white",注意这个是放在 aes() 外面的,表明我们并不是要建立映射,而是设定 color = "white",大家看到刚刚 mapdata 打印出来的结果中变量 geometry 的下面有个 <MULTIPOLYGON [m]>,这是变量的类别,表明这里面的 geometry 变量是个 MULTIPOLYGON 型的(类似的,pop 是 double 型的,NAME 是 character 型的)。使用 MULTIPOLYGON 型的数据绘制出的是面,而对于面来说,里面的填充色是 fill 指定的,边缘的颜色是 color 指定的,边缘的宽度是 size 指定的。3. size = 0.05,同样要注意到这个是在 aes() 外面的,表明我们是想设定边缘(这里就是省界)的宽度是 0.05。

还可以再添加一些文本:

R
1
2
3
4
5
6
7
8
9
10
11
12
ggplot() + 
geom_sf(data = mapdata, aes(geometry = geometry,
fill = pop),
color = "white",
size = 0.05) +
geom_sf(data = mapborder, aes(geometry = geometry),
size = 0.05) +
scale_fill_viridis_c() +
labs(title = "China Map",
subtitle = "TidyFriday",
caption = "Source: Random data",
fill = "Population")

再调整下让这个图更好看些:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
ggplot() + 
geom_sf(data = mapdata, aes(geometry = geometry,
fill = pop),
color = "white",
size = 0.05) +
geom_sf(data = mapborder, aes(geometry = geometry),
size = 0.05) +
scale_fill_viridis_c() +
labs(title = "China Map",
subtitle = "TidyFriday",
caption = "Source: Random data",
fill = "Population") +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
legend.position = c(0.2, 0.3),
legend.direction = "horizontal") +
guides(fill = guide_colorbar(title.position = "top", label.position = "bottom"))

这里我添加了很多设置,theme() 函数里面有 6 个,前两个的功能是去掉经纬网,中间两个的功能是去掉经纬度标签,最后两个的功能分别是把图例放置在离左侧 20%,下侧 30% 的位置和让图例水平放置。guides() 函数可以用于图例的更精细设置,注意这幅图的图例是使用 fill 映射生成的,所以是 fill = guide_colorbar,然后我在里面做了两个设置,标题放在图例的上方,图例上的文本标签放置在图例的下方。

组合地图

经常需要组合多幅地图,可以使用 patchwork 包:

R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
library(patchwork)
p <- ggplot() +
geom_sf(data = mapdata, aes(geometry = geometry,
fill = pop),
color = "white",
size = 0.05) +
geom_sf(data = mapborder, aes(geometry = geometry),
size = 0.05) +
scale_fill_viridis_c() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank())
{p + theme(legend.position = "none",
plot.margin = grid::unit(c(1, 0, 1, 1), "cm"))} +
{p + theme(plot.margin = grid::unit(c(1, 1, 1, 0), "cm"))}

这里为了让图更紧凑,我分别调整了两个图的边距,边距的调整是对 plot.margin 进行设置,grid::unit(c(1, 0, 1, 1), "cm")) 是一个四个元素的向量,四个数分别对应图表的“上右下左”边距(也就是顺时针方向)。legend.position = "none" 的设置可以去除图例。

质心的计算和标签的添加

还可以为每个省添加标签,标签图层需要经纬度,所以我们先根据 geometry 计算每个省的质心,然后再叠加标签图层:

R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
centroid <- st_centroid(mapdata)
centroid
#> Simple feature collection with 34 features and 2 fields
#> geometry type: POINT
#> dimension: XY
#> bbox: xmin: -1609429 ymin: 2489613 xmax: 1671853 ymax: 5866667
#> epsg (SRID): NA
#> proj4string: +proj=lcc +lat_1=30 +lat_2=62 +lat_0=0 +lon_0=105 +x_0=0 +y_0=0 +ellps=krass +units=m +no_defs
#> # A tibble: 34 x 3
#> NAME geometry pop
#> <chr> <POINT [m]> <dbl>
#> 1 黑龙江 (1623938 5866667) 162
#> 2 新疆 (-1609429 5113455) 261
#> 3 山西 (621816 4559114) 284
#> 4 宁夏 (99683.5 4501151) 617
#> 5 西藏 (-1532343 4038557) 967
#> 6 山东 (1142137 4495014) 125
#> 7 河南 (780158.7 4174676) 796
#> 8 江苏 (1325927 4156136) 830
#> 9 安徽 (1142714 3994883) 159
#> 10 湖北 (688817.5 3846593) 445
#> # … with 24 more rows

centroid 里面也有一个 geometry 变量,这也是一个 sf 对象,不过它的 geometry 变量是个 POINT 型的,[m] 是单位,米制。

R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
ggplot() + 
geom_sf(data = mapdata, aes(fill = pop),
color = "white", size = 0.05) +
geom_sf(data = mapborder, size = 0.1) +
geom_sf(data = centroid, size = 0.05) +
geom_sf_label(data = centroid,
aes(label = NAME),
family = cnfont) +
scale_fill_viridis_c() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank())

以后会单独讲解如何设置字体。Windows 系统的用户可以先抛开字体设置,因为电脑默认支持中文绘图。但是上图中的标签遮盖问题很严重,可以使用 ggrepel 可以解决这个问题,不过 ggrepel 包的 geom_label_repel 图层只能接收 XY 坐标,所以我们需要把 geometry 变量转成 XY 坐标:

R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
centroid_coords <- cbind(
centroid,
st_coordinates(centroid)
)
centroid_coords
#> Simple feature collection with 34 features and 4 fields
#> geometry type: POINT
#> dimension: XY
#> bbox: xmin: -1609429 ymin: 2489613 xmax: 1671853 ymax: 5866667
#> epsg (SRID): NA
#> proj4string: +proj=lcc +lat_1=30 +lat_2=62 +lat_0=0 +lon_0=105 +x_0=0 +y_0=0 +ellps=krass +units=m +no_defs
#> First 10 features:
#> NAME pop X Y geometry
#> 1 黑龙江 162 1623937.8 5866667 POINT (1623938 5866667)
#> 2 新疆 261 -1609428.6 5113455 POINT (-1609429 5113455)
#> 3 山西 284 621816.0 4559114 POINT (621816 4559114)
#> 4 宁夏 617 99683.5 4501151 POINT (99683.5 4501151)
#> 5 西藏 967 -1532343.0 4038557 POINT (-1532343 4038557)
#> 6 山东 125 1142136.6 4495014 POINT (1142137 4495014)
#> 7 河南 796 780158.7 4174676 POINT (780158.7 4174676)
#> 8 江苏 830 1325927.1 4156136 POINT (1325927 4156136)
#> 9 安徽 159 1142714.0 3994883 POINT (1142714 3994883)
#> 10 湖北 445 688817.5 3846593 POINT (688817.5 3846593)

rbind() 函数可以把两个等行数的数据框合并在一起。

R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
library(ggrepel)
ggplot() +
geom_sf(data = mapdata, aes(fill = pop),
color = "white", size = 0.05) +
geom_sf(data = mapborder, size = 0.1) +
geom_sf(data = centroid, size = 0.05) +
geom_label_repel(data = centroid_coords,
aes(x = X, y = Y, label = NAME),
family = cnfont) +
scale_fill_viridis_c() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank())

添加比例尺和指北针可以使用 ggspatial 包:

R
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
library(ggspatial)
ggplot() +
geom_sf(data = mapdata, aes(fill = pop),
color = "white", size = 0.05) +
geom_sf(data = mapborder, size = 0.1) +
geom_sf(data = centroid, size = 0.05) +
geom_label_repel(data = centroid_coords,
aes(x = X, y = Y, label = NAME),
family = cnfont) +
scale_fill_viridis_c() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank()) +
annotation_scale(location = "bl", width_hint = 0.3,
bar_cols = c("#0D0887FF", "#F0F921FF")) +
annotation_north_arrow(location = "bl", which_north = "true",
pad_x = unit(0.75, "cm"),
pad_y = unit(0.5, "cm"),
style = north_arrow_fancy_orienteering(
line_col = "#0D0887FF",
text_col = "#F0F921FF",
fill = c("#0D0887FF", "#F0F921FF")
))

县级地图

县级地图的绘制方法类似。

首先读入县级数据:

R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
county <- read_sf("chinamap/中国县界.shp")
# 生成随机数据
county$pop <- round(runif(2391, 100, 1000))
ggplot() +
geom_sf(data = county, aes(fill = pop),
color = "white", size = 0.05) +
geom_sf(data = mapborder, size = 0.1) +
scale_fill_viridis_c() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank())

中国主要河流

如果想展示中国的主要河流,只需要再叠加一层河流图层就好了,需要注意河流是的 geometry 是 <LINESTRING [m]> 型的。

R
1
2
3
4
5
6
7
8
9
10
11
12
river <- read_sf("chinamap/主要河流.shp")
river
ggplot() +
geom_sf(data = mapdata, aes(fill = pop),
color = "white", size = 0.05) +
geom_sf(data = mapborder, size = 0.1) +
geom_sf(data = river, size = 0.1, color = "lightblue") +
scale_fill_viridis_c() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank())

因为河流是线元素,所以使用 color 进行映射。

地图描点

上面我们演示了面元素和线元素的使用,最后就是点元素了:

R
1
2
3
4
5
6
7
8
9
10
11
12
13
city <- read_sf("chinamap/县城驻地.shp")
city$pop <- round(runif(2423, 100, 1000))
ggplot() +
geom_sf(data = mapdata, aes(geometry = geometry),
color = "black", fill = NA, size = 0.1) +
geom_sf(data = mapborder, aes(geometry = geometry),
size = 0.1) +
geom_sf(data = city, aes(geometry = geometry),
size = 0.1) +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank())

#

评论

Your browser is out-of-date!

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

×