使用R、ggplot2和sf绘制地图

使用R、ggplot2和sf绘制地图

该篇是三篇博客的学习笔记:

  1. Drawing beautiful maps programmatically with R, sf and ggplot2 — Part 1: Basics
  2. Drawing beautiful maps programmatically with R, sf and ggplot2 — Part 2: Layers
  3. Drawing beautiful maps programmatically with R, sf and ggplot2 — Part 3: Layouts

全文一共有高达500+行的代码。。。敲的累死个人!但是也学到了不少新东西!

入门

rworldmap包中提供了世界地图的数据,rworldxtra包提供了分辨率更高的地图。

1
2
3
4
5
6
library(ggplot2)
library(sf)
theme_set(theme_bw(base_size = 20, base_family = "STSong") + theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5)))

library(rworldmap)
library(rworldxtra)

1
2
3
4
5
6
7
8
9
10
11
12
13
world <- getMap(resolution = "high")
# 查看world数据集的类别
class(world)

# [1] "SpatialPolygonsDataFrame"
# attr(,"package")
# [1] "sp"

# 现在的world是个sp(SpatialPolygonsDataFrame)可以使用sf包中的st_as_sf()函数将其转换成简单的sf
world <- st_as_sf(world)
class(world)

# [1] "sf" "data.frame"

数据和基本绘图

首先我们绘制一个简单的世界地图:

1
2
ggplot(data = world) +
geom_sf()

这个图要绘制蛮长时间,需要耐心等候。

设置地图边界颜色和填充颜色:

1
2
3
4
ggplot(data = world) +
geom_sf(color = "black", fill = "lightgreen") +
labs(title = "图:世界地图",
subtitle = paste0("(", length(unique(world$NAME)), "个国家", ")"))

绘制各国人口的热力地图:

1
2
3
4
5
ggplot(data = world) +
geom_sf(aes(fill = POP_EST)) +
labs(title = "图:世界各国人口",
subtitle = paste0("(", length(unique(world$NAME)), "个国家", ")")) +
scale_fill_viridis_c("人口", option = "plasma", trans = "sqrt")

投影与范围

这可以使用任何有效的PROJ4字符串(这里是以欧洲为中心的ETRS89 Lambert Azimuthal Equal-Area投影)来实现:

1
2
3
ggplot(data = world) +
geom_sf() +
coord_sf(crs = "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs ")

暨南大学的经纬度为:(113, 23),以此为中心绘制地图:

1
2
3
ggplot(data = world) +
geom_sf() +
coord_sf(crs = "+proj=laea +lat_0=23 +lon_0=113 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs ")

可以选择感兴趣的区域实现缩放:

1
2
3
ggplot(data = world) +
geom_sf() +
coord_sf(xlim = c(140, 70), ylim = c(3, 60), expand = F)

值得注意的是,xlim和ylim参数是非常有趣的,如果是西经,xlim参数正数表示东经,负数表示东经,顺序是自西向东。对于纬度,正数表示北纬,负数表示南纬,顺序是从赤道向两级延伸。

比例尺和指北针

1
2
3
4
5
6
7
8
9
library(ggspatial)
ggplot(data = world) +
geom_sf() +
annotation_scale(location = "bl", width_hint = 0.3) +
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) +
coord_sf(xlim = c(140, 60), ylim = c(3, 60), expand = F)

注意不准确比例尺的警告:由于地图在等距圆柱投影(所有经线是平行的)上使用经度/纬度(WGS84)的未投影数据,地图上的(千米)长度直接取决于数学上的经纬度。小区域或投影数据的图通常允许更精确的比例尺。

国名和地理名

1
2
3
4
5
6
7
8
9
10
ggplot(data = world) +
geom_sf() +
geom_text(aes(x = LON, y = LAT, label = NAME), size = 4,
hjust = "left", color = "darkblue", fontface = "bold",
check_overlap = T) +
annotate(geom = "text", x = 115, y = 15, label = "中国南海",
fontface = "italic", color = "grey22",
size = 4, family = "STSong") +
coord_sf(xlim = c(140, 60), ylim = c(3, 60), expand = F) +
labs(x = "", y = "")

填充

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
ggplot(data = world) +
geom_sf(fill = "antiquewhite1") +
geom_text(aes(LON, LAT, label = NAME),
size = 4, hjust = "left", color = "darkblue",
fontface = "bold", check_overlap = T) +
annotate(geom = "text", x = 115, y = 15, label = "中国南海",
fontface = "italic", color = "grey22",
size = 4, family = "STSong") +
coord_sf(xlim = c(140, 60), ylim = c(3, 60), expand = F) +
annotation_scale(location = "bl", width_hint = 0.5) +
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) +
labs(x = "经度", y = "纬度", title = "中国地图") +
theme(panel.grid.major = element_line(color = gray(0.5),
linetype = "dashed",
size = 0.5),
panel.background = element_rect(fill = "aliceblue"))

完全中文化

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
47
48
49
50
51
52
53
54
55
country_name_df <- data.frame(
c("China", "中国"),
c("Kazakhsta", "哈萨克斯坦"),
c("Mongolia", "蒙古"),
c("Baykonur", "贝克诺"),
c("Uzbekistan", "乌兹别克斯坦"),
c("Kyrgyzstan", "吉尔吉斯斯坦"),
c("Tajikistan", "塔吉克斯坦"),
c("Afghanistan", "阿富汗"),
c("Pakistan", "巴基斯坦"),
c("N.Korea", "朝鲜"),
c("S.Korea", "韩国"),
c("Japan", "日本"),
c("Nopal", "尼泊尔"),
c("Bhutan", "不丹"),
c("Taiwan", "台湾"),
c("Arb Emirates", "阿联酋"),
c("India", "印度"),
c("Bangladesh", "孟加拉国"),
c("Myanmar", "缅甸"),
c("Laos", "老挝"),
c("Vietnam", "越南"),
c("Thailand", "台湾"),
c("Cambodia", "柬埔寨"),
c("Sri Lank", "斯里兰卡"),
c("Hong Kong", "香港"),
c("Philippines", "菲律宾"),
c("Brunei", "文莱"),
c("Palau", "帕劳")
)

country_name_df <- t(country_name_df)
colnames(country_name_df) <- c("NAME", "cname")
rownames(country_name_df) <- NULL

df <- merge(world, country_name_df, by = "NAME")
ggplot(data = df) +
geom_sf(fill = "antiquewhite1") +
geom_text(aes(LON, LAT, label = cname),
size = 4, hjust = "left", color = "darkblue",
fontface = "bold", check_overlap = T, family = "STSong") +
annotate(geom = "text", x = 115, y = 15, label = "中国南海",
fontface = "italic", color = "grey22",
size = 4, family = "STSong") +
coord_sf(xlim = c(140, 60), ylim = c(3, 60), expand = F) +
annotation_scale(location = "bl", width_hint = 0.5) +
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) +
labs(x = "经度", y = "纬度", title = "中国地图") +
theme(panel.grid.major = element_line(color = gray(0.5),
linetype = "dashed",
size = 0.5),
panel.background = element_rect(fill = "aliceblue"))

添加点

这里添加两个佛罗里达的观测点:

1
2
3
4
5
6
7
(sites <- data.frame(longitude = c(-80.144005, -80.109), latitude = c(26.479005, 26.83)))
ggplot(data = world) +
geom_sf() +
geom_point(data = sites, aes(x = longitude,
y = latitude),
size = 4, shape = 23, fill = "darkred") +
coord_sf(xlim = c(-88, -78), ylim = c(24.5, 33), expand = F)

添加多边形

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# 添加多边形
library(maps)
states <- st_as_sf(map("state", plot = F, fill = T))
# 计算每个州质心
states <- cbind(states, st_coordinates(st_centroid(states)))

# 把州名首字母大写
library(tools)
states$ID <- toTitleCase(states$ID)

ggplot(data = world) +
geom_sf() +
geom_sf(data = states, fill = NA) +
geom_text(data = states, aes(X, Y, label = ID), size = 5) +
coord_sf(xlim = c(-88, -78), ylim = c(24.5, 33), expand = F)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# 稍微移动州名,以便更好地显示南卡罗来纳州和佛罗里达州
# 为此,我们需要创建一个新变量nudge_y,对于所有州略微向南移动1,佛罗里达州向南移动0.5, 南卡罗来纳州向南移动1.5
states$nudge_y <- -1
states$nudge_y[states$ID == "Florida"] <- -0.5
states$nudge_y[states$ID == "South Carolina"] <- -1.5

# 为了提高可读性,我们在州的外边绘制一个矩形,使用geom_label代替geom_text,然后再次绘图
ggplot(data = world) +
geom_sf() +
geom_sf(data = states, fill = NA) +
geom_label(data = states, aes(X, Y, label = ID),
size = 5, fontface = "bold",
nudge_y = states$nudge_y) +
coord_sf(xlim = c(-88, -78), ylim = c(24.5, 33), expand = F)

添加县

1
2
3
4
5
6
7
8
9
10
11
# 添加县
counties <- st_as_sf(map("county", plot = F, fill = T))
counties <- subset(counties, grepl("florida", counties$ID))
# 计算县的面积
counties$area <- as.numeric(st_area(counties))

ggplot(data = world) +
geom_sf() +
geom_sf(data = counties, aes(fill = area)) +
scale_fill_viridis_c("面积", trans = "sqrt", alpha = 0.4) +
coord_sf(xlim = c(-88, -78), ylim = c(24.5, 33), expand = F)

添加佛罗里达州的五个大城市

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
flcities <- data.frame(state = rep("Florida", 5),
city = c("Miami", "Tampa",
"Orlando", "Jacksonville",
"Sarasota"),
lat = c(25.7616798, 27.950575,
28.5383355, 30.3321838, 27.3364347),
lng = c(-80.1917902, -82.4571776, -81.3792365,
-81.655651, -82.5306527))

# 把数据框转换成sf格式
(flcities <- st_as_sf(flcities, coords = c("lng", "lat"), remove = F,
crs = 4326, agr = "constant"))

# 绘图
ggplot(data = world) +
geom_sf() +
geom_sf(data = counties, fill = NA, color = gray(0.5)) +
geom_sf(data = flcities) +
geom_text(data = flcities, aes(x = lng, y = lat, label = city),
size = 3.9, col = "black", fontface = "bold") +
coord_sf(xlim = c(-88, -78), ylim = c(24.5, 33), expand = F)

显然这个标签是让人很不舒服的,下面使用ggrepel包添加标签:

1
2
3
4
5
6
7
8
9
library(ggrepel)
ggplot(data = world) +
geom_sf() +
geom_sf(data = counties, fill = NA, color = gray(0.5)) +
geom_sf(data = flcities) +
geom_text_repel(data = flcities, aes(x = lng, y = lat, label = city),
fontface = "bold", nudge_x = c(1, -1.5, 2, 2, -1),
nudge_y = c(0.25, -0.25, 0.5, 0.5, -0.5)) +
coord_sf(xlim = c(-88, -78), ylim = c(24.5, 33), expand = F)

最终的地图

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
library(ggspatial)
ggplot(data = world) +
geom_sf(fill = "antiquewhite1") +
geom_sf(data = counties, aes(fill = area)) +
geom_sf(data = states, fill = NA) +
geom_point(data = sites, aes(x = longitude, y = latitude), size = 4, shape = 23, fill = "darkred") +
geom_sf(data = flcities) +
geom_text_repel(data = flcities,
aes(x = lng, y = lat, label = city),
fontface = "bold", nudge_x = c(1, -1.5, 2, 2, -1),
nudge_y = c(0.25, -0.25, 0.5, 0.5, -0.5)) +
geom_label(data = states, aes(X, Y, label = ID),
size = 5, fontface = "bold",
nudge_y = states$nudge_y) +
scale_fill_viridis_c("面积", trans = "sqrt", alpha = 0.4) +
annotation_scale(location = "bl", width_hint = 0.4) +
annotation_north_arrow(location = "bl", which_north = "true",
pad_x = unit(0.75, "in"),
pad_y = unit(0.5, "in"),
style = north_arrow_fancy_orienteering) +
coord_sf(xlim = c(-88, -78), ylim = c(24.5, 33), expand = F) +
labs(x = "经度", y = "纬度",
title = "图:观测地点",
subtitle = "位于佛罗里达州棕榈海滩的两个观测点") +
theme(panel.grid.major = element_line(color = gray(0.5),
linetype = "dashed",
size = 0.5),
panel.background = element_rect(fill = "aliceblue"))

创建组合图形

有两种组合子图的方案:

  1. 使用Grobs(图形对象,仅允许绘图区域中的绘图,基于坐标),直接使用ggplot2;
  2. ggdraw从包中使用(允许任何地方的地块,包括外边距,基于相对位置),使用cowplot实现。

使用grob

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
g1 <- qplot(0:10, 0:10)
g1_void <- g1 + theme_void() + theme(panel.border = element_rect(colour = "black", fill = NA))
g1 + annotation_custom(
grob = ggplotGrob(g1_void),
xmin = 0,
xmax = 3,
ymin = 5,
ymax = 10
) +
annotation_custom(
grob = ggplotGrob(g1_void),
xmin = 5,
xmax = 10,
ymin = 0,
ymax = 3
)
1
2
3
4
5
6
7
8
9
10
11
12
13
# 使用ggdraw
library(cowplot)
ggdraw(g1) +
draw_plot(g1_void,
width = 0.25,
height = 0.5,
x = 0.02,
y = 0.48) +
draw_plot(g1_void,
width = 0.5,
height = 0.25,
x = 0.75,
y = 0.09)

使用ggdraw

1
2
3
4
5
6
7
8
9
10
11
12
library(cowplot)
ggdraw(g1) +
draw_plot(g1_void,
width = 0.25,
height = 0.5,
x = 0.02,
y = 0.48) +
draw_plot(g1_void,
width = 0.5,
height = 0.25,
x = 0.75,
y = 0.09)

组合地图

准备子图1:世界地图

1
2
3
4
5
6
7
8
9
10
11
12
13
# 简化REGION[7]为
levels(world$REGION)[7] <- "South America"
(
gworld <- ggplot(data = world) +
geom_sf(aes(fill = REGION)) +
geom_rect(xmin = -102.15, xmax = -74.12,
ymin = 7.65, ymax = 33.97,
fill = NA, colour = "black", size = 1.5) +
scale_fill_viridis_d("大洲", option = "plasma") +
theme_bw(base_size = 20, base_family = "STSong") +
theme(panel.background = element_rect(fill = "azure"),
panel.border = element_rect(fill = NA))
)

子图2: 墨西哥湾

1
2
3
4
5
6
7
8
9
10
11
12
(
ggulf <- ggplot(data = world) +
geom_sf(aes(fill = REGION)) +
annotate(geom = "text", x = -90, y = 26, label = "墨西哥湾", family = "STSong", fontface = "italic", color = "grey22", size = 6) +
coord_sf(xlim = c(-102.5, -74.12), ylim = c(7.65, 33.97), expand = F) +
scale_fill_viridis_d(option = "plasma") +
theme(legend.position = "none",
axis.title.x = element_blank(),
axis.title.y = element_blank(),
panel.background = element_rect(fill = "azure"),
panel.border = element_rect(fill = NA))
)

合并

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
ggplot() +
coord_equal(xlim = c(0, 3.3),
ylim = c(0, 1),
expand = F) +
annotation_custom(ggplotGrob(gworld),
xmin = 0,
xmax = 2.3,
ymin = 0,
ymax = 1) +
annotation_custom(ggplotGrob(ggulf),
xmin = 2.3,
xmax = 3.3,
ymin = 0,
ymax = 1) +
theme_void(base_size = 20, base_family = "STSong")

美国地图的组合

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
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
# 首先是美国主体
# EPSG是一套地理坐标系统,2163对应的区域可以参考这个网页:https://epsg.io/2163
usa <- subset(world, ADMIN == "United States of America")
(mainland <- ggplot(data = usa) +
geom_sf(fill = "cornsilk") +
coord_sf(crs = st_crs(2163),
xlim = c(-2500000, 2500000),
ylim = c(-2300000, 730000)))

# 阿拉斯加地图:3467
# datum = NA用于删除刻度和坐标
(alaska <- ggplot(data = usa) +
geom_sf(fill = "cornsilk") +
coord_sf(crs = st_crs(3467),
xlim = c(-2400000, 1600000),
ylim = c(200000, 2500000),
expand = F, datum = NA))

# 夏威夷地图:4135
(hawaii <- ggplot(data = usa) +
geom_sf(fill = "cornsilk") +
coord_sf(crs = st_crs(4135),
xlim = c(-161, -154),
ylim = c(18, 23),
expand = F, datum = NA))

# 使用cowplot包中的ggdraw进行组合:
# 首先计算插图长宽比率
(ratioAlaska <- (25 - 2)/(16 + 24))
(ratioHawaii <- (23 - 18)/(-154 + 161))

ggdraw(mainland) +
draw_plot(alaska, width = 0.26,
height = 0.26 * 10/6 *ratioAlaska,
x = 0.05, y = 0.05) +
draw_plot(hawaii, width = 0.15, height = 0.15 * 10/6 * ratioHawaii, x = 0.3, y = 0.05)

# 可以使用grobs创建相同类型的绘图ggplotGrob,(注意使用xdiff / ydiff和任意比率):
mainland +
annotation_custom(
grob = ggplotGrob(alaska),
xmin = -2750000,
xmax = -2750000 + (1600000 - (-2400000))/2.5,
ymin = -2450000,
ymax = -2450000 + (2500000 - 200000)/2.5
) +
annotation_custom(
grob = ggplotGrob(hawaii),
xmin = -1250000,
xmax = -1250000 + (-154 - (-161))*120000,
ymin = -2450000,
ymax = -2450000 + (23 - 18)*120000
)

# 还可以使用print()函数进行组合
library(grid)
vp1 <- viewport(width = 1,
height = 0.30,
x = 0.15,
y = 0.10,
just = c("bottom"))
vp2 <- viewport(width = 1,
height = 0.30,
x = 0.30,
y = 0.10,
just = c("bottom"))
print(mainland)
print(alaska, vp = vp1)
print(hawaii, vp = vp2)

几个地图用箭头连接

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
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
sites <- st_as_sf(
data.frame(
lon = c(-80.15, -80.1),
lat = c(26.5, 26.8)
),
coords = c("lon", "lat"), crs = 4326,
agr = "constant"
)

# 佛罗里达地图
theme_set(theme_bw(base_family = "STSong", base_size = 15))
(
florida <- ggplot(data = world) +
geom_sf(fill = "antiquewhite1") +
geom_sf(data = sites, size = 4, shape = 23,
fill = "darkred") +
annotate("text", x = -85.5, y = 27.5, label = "墨西哥湾", family = "STSong", color = "grey22", size = 4.5) +
coord_sf(xlim = c(-87.35, -79.5), ylim = c(24.1, 30.8)) +
labs(x = "经度", y = "纬度") +
theme(panel.grid.major = element_line(colour = gray(0.5), linetype = "dashed", size = 0.5),
panel.background = element_rect(fill = "aliceblue"),
panel.border = element_rect(fill = NA))
)

# 站点A
(
siteA <- ggplot(data = world) +
geom_sf(fill = "antiquewhite1") +
geom_sf(data = sites, size = 4, shape = 23, fill = "darkred") +
coord_sf(xlim = c(-80.25, -79.95),
ylim = c(26.65, 26.95),
expand = F) +
annotate("text", x = -80.18, y = 26.92,
label = "站点A", family = "STSong",
size = 6) +
theme_void(base_family = "STSong", base_size = 15) +
theme(panel.grid.major = element_line(colour = gray(0.5), linetype = "dashed", size = 0.5),
panel.background = element_rect(fill = "aliceblue"),
panel.border = element_rect(fill = NA))
)

# 站点B
(
siteB <- ggplot(data = world) +
geom_sf(fill = "antiquewhite1") +
geom_sf(data = sites, size = 4, shape = 23, fill = "darkred") +
coord_sf(xlim = c(-80.3, -80),
ylim = c(26.35, 26.65),
expand = F) +
annotate("text", x = -80.23, y = 26.62,
label = "站点B", family = "STSong",
size = 6) +
theme_void(base_family = "STSong", base_size = 15) +
theme(panel.grid.major = element_line(colour = gray(0.5), linetype = "dashed", size = 0.5),
panel.background = element_rect(fill = "aliceblue"),
panel.border = element_rect(fill = NA))
)

# 两个箭头
arrowA <- data.frame(x1 = 18.5, x2 = 23, y1 = 9.5, y2 = 14.5)
arrowB <- data.frame(x1 = 18.5, x2 = 23, y1 = 8.5, y2 = 6.5)

# 组合
ggplot() +
coord_equal(xlim = c(0, 28), ylim = c(0, 20), expand = F) +
annotation_custom(ggplotGrob(florida),
xmin = 0, xmax = 20,
ymin = 0, ymax = 20) +
annotation_custom(ggplotGrob(siteA),
xmin = 20, xmax = 28,
ymin = 11.25, ymax = 19) +
annotation_custom(ggplotGrob(siteB),
xmin = 20, xmax = 28,
ymin = 2.5, ymax = 10.25) +
geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2), data = arrowA, arrow = arrow(), lineend = "round") +
geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2), data = arrowB, arrow = arrow(), lineend = "round") +
theme_void()

## ggdraw也可以产生类似的效果
ggdraw(xlim = c(0, 28), ylim = c(0, 20)) +
draw_plot(florida, x = 0, y = 0, width = 20, height = 20) +
draw_plot(siteA, x = 20, y = 11.25, width = 8, height = 8) +
draw_plot(siteB, x = 20, y = 2.5, width = 8, height = 8) +
geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2), data = arrowA, arrow = arrow(), lineend = "round") +
geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2), data = arrowB, arrow = arrow(), lineend = "round")

# R

评论

程振兴

程振兴 @czxa.top
截止今天,我已经在本博客上写了329.8k个字了!

庆祝
计算机二级MySQL考试通过!
成功
网站成功开启https!
Your browser is out-of-date!

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

×

keyboard_arrow_up