使用 R 语言绘图展示中国土地覆盖类型分布和各省市区县耕地占比分布

在之前的数据分享中我们分享了「1985～2021 年中国各省市区县土地覆盖类型占比面板数据」，数据介绍中我绘制了一些图表展示中国土地覆盖类型分布和各省市区县占比分布，为了方便大家在论文中使用这些图表，这里我简要介绍下这些图表的绘制方法。

首先是中国土地覆盖类型分布。我使用的是 R 语言基础绘图系统绘制的：

注意这里需要先在 Profile 里面设置基础绘图系统主题，可以参考系列课程「R 语言数据科学」的第一次课进行配置。

加载所需 R 包：

library(tidyverse)
library(terra)

绘图统一使用栅格数据的 crs：

mycrs <- "+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"

读取所需数据：

rast("CLCD_v01_2020_albert.tif") -> rst
vect("九段线.geojson") %>%
  terra::project(mycrs) -> jdx
vect("海岸线/海岸线.shp") %>%
  terra::project(mycrs) -> hax
readxl::read_xlsx("CLCD_classificationsystem.xlsx") %>%
  select(1:2) %>%
  set_names(c("kind", "class")) -> classdf
rst[rst == 0] <- NA

如果想使用 png 格式保存绘图结果，可以：

png("2020年中国土地覆盖类型分布.png",
    width = 4000, height = 4000, res = 600)
terra::plot(jdx, axes = F)
title(main = "2020 年中国土地覆盖类型分布", cex.main = 1.5)
terra::plot(rst, col = c("#FAE39C", "#446F33", "#33A02C", "#ABD37B", "#1E69B4", "#A6CEE3", "#CFBDA3", "#E24290", "#289BE8"),
            type = "classes", legend = "bottomleft", axes = F,
            add = T, plg = list(legend = c(classdf$class)))
terra::plot(jdx, add = T, lwd = 1.5, axes = F)
terra::plot(hax, add = T, lwd = 1, col = "#709ae1", axes = F)
dev.off()

这里 axes = F 表示取消横轴和纵轴的文本标签。

如果想使用 pdf 格式保存绘图结果，可以：

cairo_pdf("2020年中国土地覆盖类型分布.pdf",
          width = 7, height = 7)
terra::plot(jdx, axes = F)
title(main = "2020 年中国土地覆盖类型分布", cex.main = 1.5)
terra::plot(rst, col = c("#FAE39C", "#446F33", "#33A02C", "#ABD37B", "#1E69B4", "#A6CEE3", "#CFBDA3", "#E24290", "#289BE8"),
            type = "classes", legend = "bottomleft", axes = F,
            add = T, plg = list(legend = c(classdf$class)))
terra::plot(jdx, add = T, lwd = 1.5, axes = F)
terra::plot(hax, add = T, lwd = 1, col = "#709ae1", axes = F)
dev.off()

然后我们再绘制各省市区县耕地占比分布，实际上就是填充地图的绘制。

加载所需 R 包和读取数据：

# 各省份耕地占比
library(ggplot2)
library(ggspatial)
library(hrbrthemes)
library(sf)
# 字体设置
# 这里可能需要把 song.otf 文件替换成大家自己电脑上的字体：
library(showtext)
showtext_auto(enable = TRUE)
font_add("Songti", regular = "song.otf")
cnfont <- "Songti"
read_sf("九段线.geojson") -> cn
read_sf("海岸线/海岸线.shp") -> hax
read_sf("2020行政区划/省.geojson") %>%
  st_make_valid() -> provmap

省份的。读取“2020年中国各省份土地覆被类型占比（百分比）.xlsx”：

readxl::read_xlsx("2020年中国各省份土地覆被类型占比（百分比）.xlsx") -> df

把 df 和地图数据连接起来：

provmap %>%
  left_join(df) -> prov1

把 Cropland 变量切分成分组变量进行绘图：

prov1 %>%
  mutate(group = cut(Cropland,
                       breaks = seq(0, 70, by = 10),
                       labels = c("0——10", "10——20",
                                  "20——30", "30——40",
                                  "40——50", "50——60",
                                  "60——70"))) %>%
  ggplot() +
  geom_sf(aes(fill = group), linewidth = 0.1,
          color = "black") +
  geom_sf(data = cn, linewidth = 0.5,
          color = "black", fill = NA) +
  geom_sf(data = hax, color = "#0055AA", linewidth = 0.5) +
  geom_sf_text(aes(label = 省), family = cnfont, size = 3,
               color = "gray50") +
  coord_sf(crs = mycrs,
           xlim = c(-3500000, 3090000)) +
  scico::scale_fill_scico_d(
    palette = "bamako", direction = -1,
    name = "2020 年中国各省份耕地占比分布（%）",
    guide = guide_legend(ncol = 2, byrow = T)
  ) +
  scale_x_continuous(expand = c(0.02, 0.02)) +
  scale_y_continuous(expand = c(0.02, 0.02)) +
  annotation_scale(
    width_hint = 0.2,
    text_family = "Songti",
    pad_x = unit(0.5, "cm")
  ) +
  annotation_north_arrow(
    location = "tr", which_north = "false",
    width = unit(1.6, "cm"),
    height = unit(2, "cm"),
    style = north_arrow_fancy_orienteering(
      text_family = "Songti"
    )
  ) +
  theme_ipsum(base_family = cnfont, grid = F) +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.position = c(0.25, 0.2),
        plot.background = element_rect(fill = "white",
                                       color = "white")) -> p1
ggsave(filename = "2020 年中国各省份耕地占比分布.png", device = png, width = 8, height = 8)

城市的：

# 城市
read_sf("2020行政区划/市.geojson") %>%
  st_make_valid() -> citymap

readxl::read_xlsx("2020年中国各城市土地覆被类型占比（百分比）.xlsx") %>%
  mutate(省代码 = as.numeric(省代码)) -> df2

citymap %>%
  left_join(df2) -> city1

city1 %>%
  mutate(group = cut(Cropland,
                     breaks = seq(0, 100, by = 10),
                     labels = c("0——10", "10——20",
                                "20——30", "30——40",
                                "40——50", "50——60",
                                "60——70", "70——80",
                                "80——90", "90——100"))) %>%
  dplyr::filter(!is.na(group)) %>%
  ggplot() +
  geom_sf(aes(fill = group), linewidth = 0.1,
          color = "black") +
  geom_sf(data = cn, linewidth = 0.5,
          color = "black", fill = NA) +
  geom_sf(data = hax, color = "#0055AA", linewidth = 0.5) +
  coord_sf(crs = mycrs ,
           xlim = c(-3500000, 3090000)) +
  scico::scale_fill_scico_d(
    palette = "bamako", direction = -1,
    name = "2020 年中国各城市耕地占比分布（%）",
    guide = guide_legend(ncol = 2, byrow = T)
  ) +
  scale_x_continuous(expand = c(0.02, 0.02)) +
  scale_y_continuous(expand = c(0.02, 0.02)) +
  annotation_scale(
    width_hint = 0.2,
    text_family = "Songti",
    pad_x = unit(0.5, "cm")
  ) +
  annotation_north_arrow(
    location = "tr", which_north = "false",
    width = unit(1.6, "cm"),
    height = unit(2, "cm"),
    style = north_arrow_fancy_orienteering(
      text_family = "Songti"
    )
  ) +
  theme_ipsum(base_family = cnfont, grid = F) +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.position = c(0.25, 0.2),
        plot.background = element_rect(fill = "white",
                                       color = "white")) -> p2
ggsave(filename = "2020 年中国各城市耕地占比分布.png", device = png, width = 8, height = 8)

区县的：

# 区县
read_sf("2020行政区划/县.geojson") %>%
  st_make_valid() -> countymap

readxl::read_xlsx("2020年中国各区县土地覆被类型占比（百分比）.xlsx") %>%
  mutate(省代码 = as.numeric(省代码),
            市代码 = as.numeric(市代码)) -> df3

countymap %>%
  left_join(df3) -> county1

county1 %>%
  mutate(group = cut(Cropland,
                     include.lowest = T,
                     breaks = seq(0, 100, by = 10),
                     labels = c("0——10", "10——20",
                                "20——30", "30——40",
                                "40——50", "50——60",
                                "60——70", "70——80",
                                "80——90", "90——100"))) %>%
  dplyr::filter(!is.na(group)) %>%
  ggplot() +
  geom_sf(aes(fill = group), linewidth = 0.1,
          color = "black") +
  geom_sf(data = cn, linewidth = 0.5,
          color = "black", fill = NA) +
  geom_sf(data = hax, color = "#0055AA", linewidth = 0.5) +
  coord_sf(crs = mycrs,
           xlim = c(-3500000, 3090000)) +
  scico::scale_fill_scico_d(
    palette = "bamako", direction = -1,
    name = "2020 年中国各区县耕地占比分布（%）",
    guide = guide_legend(ncol = 2, byrow = T)
  ) +
  scale_x_continuous(expand = c(0.02, 0.02)) +
  scale_y_continuous(expand = c(0.02, 0.02)) +
  annotation_scale(
    width_hint = 0.2,
    text_family = "Songti",
    pad_x = unit(0.5, "cm")
  ) +
  annotation_north_arrow(
    location = "tr", which_north = "false",
    width = unit(1.6, "cm"),
    height = unit(2, "cm"),
    style = north_arrow_fancy_orienteering(
      text_family = "Songti"
    )
  ) +
  theme_ipsum(base_family = cnfont, grid = F) +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.position = c(0.25, 0.2),
        plot.background = element_rect(fill = "white",
                                       color = "white")) -> p2
ggsave(filename = "2020 年中国各区县耕地占比分布.png", device = png, width = 8, height = 8)

点击这里跳转到 RStata 短书平台获取附件：使用 R 语言绘图展示中国土地覆盖类型分布和各省市区县耕地占比分布