Create Complete China Maps Using GGPLOT2 and SF

Create Complete China Maps Using GGPLOT2 and SF

I always feel disappointment about how to plot a beautiful China map. Yesterday, I finally drew a beautiful China map, at least in my opinion.

I make a random data file for this map: dataset.xlsx, and here is map data: chinamap.zip.

The final result is:

Map theme

theme.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
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
theme_map <- function(...) {
theme_minimal() +
theme(
text = element_text(family = "TimesNewRomanPSMT",
color = "black"),
# remove all axes
axis.line = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
# add a subtle grid
panel.grid.major = element_line(color = "#dbdbd9", size = 0.2),
panel.grid.minor = element_blank(),
# background colors
plot.background = element_rect(fill = "white",
color = NA),
panel.background = element_rect(fill = "white",
color = NA),
legend.background = element_rect(fill = "white",
color = NA),
# borders and margins
plot.margin = unit(c(.5, .5, .2, .5), "cm"),
panel.border = element_blank(),
panel.spacing = unit(c(-.1, 0.2, .2, 0.2), "cm"),
# titles
legend.title = element_text(size = 11),
legend.text = element_text(size = 9, hjust = 0.1, vjust = 1,
color = "black"),
plot.title = element_text(size = 15, hjust = 0.1,vjust = 1,
color = "black"),
plot.subtitle = element_text(size = 10, hjust = 0.1,
color = "black",
margin = margin(b = -0.1,
t = -0.1,
l = 2,
unit = "cm"),
debug = F),
# captions
plot.caption = element_text(size = 7,
hjust = .5,
margin = margin(t = 0.2,
b = 0,
unit = "cm"),
color = "#939184"),
...
)
}

grid_arrange_shared_legend <- function(..., ncol = length(list(...)), nrow = 1, position = c("bottom", "right")) {

plots <- list(...)
position <- match.arg(position)
g <- ggplotGrob(plots[[2]] + theme(legend.position = position, legend.title = element_blank()))$grobs
legend <- g[[which(sapply(g, function(x) x$name) == "guide-box")]]
lheight <- sum(legend$height)
lwidth <- sum(legend$width)
gl <- lapply(plots, function(x) x + theme(legend.position="none"))
gl <- c(gl, ncol = ncol, nrow = nrow)

combined <- switch(position,
"bottom" = arrangeGrob(do.call(arrangeGrob, gl),
legend,
ncol = 1,
heights = unit.c(unit(1, "npc") - lheight, lheight)),
"right" = arrangeGrob(do.call(arrangeGrob, gl),
legend,
ncol = 2,
widths = unit.c(unit(1, "npc") - lwidth, lwidth)))

grid.newpage()
grid.draw(combined)

# return gtable invisibly
invisible(combined)
}

grid_arrange_shared_legend() function is used to compined two gg object which share a common legend.

Load Revelant Packages and Source theme.R

1
2
3
4
5
6
7
8
9
10
11
12
13
14
library(ggplot2)
library(maptools)
library(sf)
library(dplyr)
library(readxl)
library(ggspatial)
library(fishualize)
library(viridis)
library(gridExtra)
library(grid)
library(cowplot)
library(magrittr)
library(purrr)
source("theme.R")

Load map data:

1
2
mapdata <- st_read("chinamap/中国省界.shp")
mapborder <- st_read("chinamap/国界线.shp")

Figure 2

Figure 2 display data in 2017, values in which are bigger than that in 1990 (Figure 1). Since I want to use one common legend for figure 1 and figure 2, so I plot figure 2 first.

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
df2 <- readxl::read_xlsx("dataset.xlsx", 
sheet = 2) %>%
select(1:2)

colnames(df2) <- c("NAME", "a2017")
df2$NAME[df2$NAME == "宁复"] <- "宁夏"

df2 <- df2 %>%
add_row(NAME = "台湾", a2017 = NA) %>%
left_join(mapdata) %>%
as_tibble() %>%
dplyr::filter(NAME != "全国") %>%
dplyr::filter(NAME != "跨区")
df2$a2017[df2$a2017 < 1] <- NA

# define number of classes
no_classes <- 6

# extract quantiles
quantiles2 <- df2 %>%
pull(a2017) %>%
quantile(probs = seq(0, 1, length.out = no_classes + 1), na.rm = TRUE) %>%
as.vector() # to remove names of quantiles, so idx below is numeric
quantiles2[1] <- 0
# here we create custom labels
labels2 <- imap_chr(quantiles2, function(., idx){
return(paste0(round(quantiles2[idx] / 1000000, 1),
"m",
" – ",
round(quantiles2[idx + 1] / 1000000, 1),
"m"))
})

# we need to remove the last label
# because that would be something like "478k - NA"
labels2 <- labels2[1:length(labels2) - 1]

# here we actually create a new
# variable on the dataset with the quantiles
df2 %<>%
mutate(a2017 = cut(a2017,
breaks = quantiles2,
labels = labels2,
include.lowest = T))

p2 <- ggplot(mapborder) +
geom_sf(size = 0.1) +
geom_sf(data = df2,
aes(fill = a2017, geometry = geometry), size = 0.1, alpha = 0.8) +
scale_fill_manual(
values = RColorBrewer::brewer.pal(6, "RdBu"),
guide = guide_legend(
keyheight = unit(5, units = "mm"),
title.position = "top",
reverse = T
)) +
theme_map() +
labs(title = "(b)")

Figure 1

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
df1 <- readxl::read_xlsx("dataset.xlsx", 
sheet = 1) %>%
select(1:2)

colnames(df1) <- c("NAME", "a1990")
df1$NAME[df1$NAME == "宁复"] <- "宁夏"

df1 <- df1 %>% left_join(mapdata) %>%
as_tibble() %>%
dplyr::filter(NAME != "全国") %>%
dplyr::filter(NAME != "跨区")
df1$a1990[df1$a1990 < 1] <- NA

# here we actually create a new
# variable on the dataset with the quantiles
df1 %<>%
mutate(a1990 = cut(a1990,
breaks = quantiles2,
labels = labels2,
include.lowest = T))

p1 <- ggplot(mapborder) +
geom_sf(size = 0.1) +
geom_sf(data = df1,
aes(fill = a1990, geometry = geometry), size = 0.1, alpha = 0.8) +
scale_fill_manual(
values = RColorBrewer::brewer.pal(6, "RdBu"),
guide = guide_legend(
keyheight = unit(5, units = "mm"),
title.position = "top",
reverse = T
)) +
theme_map() +
labs(title = "(a)")

Figure 4

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
df4 <- readxl::read_xlsx("dataset.xlsx", 
sheet = 4) %>%
select(1:2)

colnames(df4) <- c("NAME", "a2017")
df4$NAME[df4$NAME == "宁复"] <- "宁夏"

df4 <- df4 %>% left_join(mapdata) %>%
as_tibble() %>%
dplyr::filter(NAME != "全国") %>%
dplyr::filter(NAME != "跨区")
df4$a2017[df4$a2017 < 1] <- NA

# extract quantiles

no_classes <- 6
quantiles4 <- df4 %>%
pull(a2017) %>%
quantile(probs = seq(0, 1, length.out = no_classes + 1), na.rm = TRUE) %>%
as.vector() # to remove names of quantiles, so idx below is numeric
quantiles4[1] <- 0

# here we create custom labels
library(purrr)
labels4 <- imap_chr(quantiles4, function(., idx){
return(paste0(round(quantiles4[idx] / 1000, 2),
"k",
" – ",
round(quantiles4[idx + 1] / 1000, 2),
"k"))
})

# we need to remove the last label
# because that would be something like "478k - NA"
labels4 <- labels4[1:length(labels4) - 1]

# here we actually create a new
# variable on the dataset with the quantiles
df4 %<>%
mutate(a2017 = cut(a2017,
breaks = quantiles4,
labels = labels4,
include.lowest = T))

p4 <- ggplot(mapborder) +
geom_sf(size = 0.1) +
geom_sf(data = df4,
aes(fill = a2017, geometry = geometry), size = 0.1, alpha = 0.8) +
scale_fill_manual(
values = RColorBrewer::brewer.pal(6, "RdBu"),
guide = guide_legend(
keyheight = unit(5, units = "mm"),
title.position = "top",
reverse = T
)) +
theme_map() +
labs(title = "(d)")

Figure 3

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
df3 <- readxl::read_xlsx("dataset.xlsx", 
sheet = 3) %>%
select(1:2)

colnames(df3) <- c("NAME", "a1990")
df3$NAME[df3$NAME == "宁复"] <- "宁夏"

df3 <- df3 %>% left_join(mapdata) %>%
as_tibble() %>%
dplyr::filter(NAME != "全国") %>%
dplyr::filter(NAME != "跨区")
df3$a1990[df3$a1990 < 1] <- NA


# here we actually create a new
# variable on the dataset with the quantiles
df3 %<>%
mutate(a1990 = cut(a1990,
breaks = quantiles4,
labels = labels4,
include.lowest = T))

p3 <- ggplot(mapborder) +
geom_sf(size = 0.1) +
geom_sf(data = df3,
aes(fill = a1990, geometry = geometry), size = 0.1, alpha = 0.8) +
scale_fill_manual(
values = RColorBrewer::brewer.pal(6, "RdBu"),
guide = guide_legend(
keyheight = unit(5, units = "mm"),
title.position = "top",
reverse = T
)) +
theme_map() +
labs(title = "(c)")

Figure 5

Figure 5 is a combined plot with two y axes. Actually ggplot2 doesn’t provide this function, this is accomplished by cowplot:

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
library(tidyr)
(df5 <- read_excel("dataset.xlsx", sheet = 5, skip = 53) %>%
slice(1:8) %>%
select(2:30) %>%
rename(region = year) %>%
gather(`1990`, `1991`, `1992`, `1993`, `1994`,
`1995`, `1996`, `1997`, `1998`, `1999`,
`2000`, `2001`, `2002`, `2003`, `2004`,
`2005`, `2006`, `2007`, `2008`, `2009`,
`2010`, `2011`, `2012`, `2013`, `2014`,
`2015`, `2016`, `2017`, key = "year", value = "emission") %>%
type.convert())

(q1 <- ggplot(df5) +
geom_col(aes(x = factor(year), y = emission,
color = region, fill = region),
position = position_fill()) +
scale_fill_manual(
values = RColorBrewer::brewer.pal(8, "RdBu"),
guide = guide_legend(
keyheight = unit(5, units = "mm"),
title.position = "top",
reverse = T
)) +
scale_color_manual(
values = RColorBrewer::brewer.pal(8, "RdBu"),
guide = guide_legend(
keyheight = unit(5, units = "mm"),
title.position = "top",
reverse = T
)) +
scale_y_percent() +
theme_cowplot(font_family = "TimesNewRomanPSMT") +
theme(axis.text.x = element_text(angle = 90, size = 10),
axis.text.y = element_text(size = 10),
axis.title.y = element_text(size = 10),
legend.position = "bottom") +
labs(x = element_blank(), y = "Cross-region Power Grid",
color = element_blank(), fill = element_blank()) +
theme(panel.grid.major = element_blank()) +
theme(panel.grid.minor = element_blank()))

(df6 <- read_excel("dataset.xlsx", sheet = 5, skip = 53) %>%
slice(10) %>%
select(2:30) %>%
gather(`1990`, `1991`, `1992`, `1993`, `1994`,
`1995`, `1996`, `1997`, `1998`, `1999`,
`2000`, `2001`, `2002`, `2003`, `2004`,
`2005`, `2006`, `2007`, `2008`, `2009`,
`2010`, `2011`, `2012`, `2013`, `2014`,
`2015`, `2016`, `2017`, key = "year", value = "gini") %>%
type.convert())

(q2 <- df6 %>%
ggplot() +
geom_line(aes(x = year, y = gini)) +
labs(y = "GINI coefficient per capita") +
theme_cowplot(font_family = "TimesNewRomanPSMT") +
theme(axis.title.y = element_text(size = 10),
axis.text.x = element_blank(),
axis.title.x = element_blank(),
axis.text.y = element_text(size = 10)) +
scale_x_discrete(limits = c(1990, 2017)) +
scale_y_continuous(position = "right") +
theme(panel.grid.major = element_blank()) +
theme(panel.grid.minor = element_blank()) +
theme(axis.ticks.x = element_blank()))


aligned_plots <- align_plots(q1, q2, align="hv", axis="tblr")

ggdraw(aligned_plots[[1]]) + draw_plot(aligned_plots[[2]])

Combine all!

1
2
3
4
5
6
7
8
9
10
plot2by2 <- plot_grid(
grid_arrange_shared_legend(p1, p2, position = "right"),
grid_arrange_shared_legend(p3, p4, position = "right"),
ggdraw(aligned_plots[[1]]) + draw_plot(aligned_plots[[2]]),
align = "v", ncol = 1
)

save_plot("plot2by2.png", plot2by2,
ncol = 1, nrow = 2,
base_aspect_ratio = 2)

unsplash-logoJames Coleman

# R

Comments

Your browser is out-of-date!

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

×