Rstudio学习|用已有的底图填充数据绘制地图_底图在绘制地图

1.矢量栅格数据常用两个包

2.设置坐标系

3.导入设置好的底图

library(tidyverse) 
library(sf)
 
# 中国地图通常使用这样的坐标系
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" 
 
# 1. 带小地图版本
 
# 读取小地图版本的中国省级地图
read_sf("E:/R_codedata/SO2与机器人安装密度/SO2与机器人安装密度/chinaprov2019mini/chinaprov2019mini.shp") %>% 
  filter(!is.na(省代码)) -> provmap
# 线条
read_sf("E:/R_codedata/SO2与机器人安装密度/SO2与机器人安装密度/chinaprov2019mini/chinaprov2019mini_line.shp") %>% 
  filter(class %in% c("九段线", "海岸线", "小地图框格")) %>% 
  select(class) -> provlinemap

4.读入数据

readxl::read_xlsx("画图.xlsx") %>% 
  gather(2:15, key = "key", value = "value") %>% #横坐标轴数据（年份）
  mutate(key = str_remove_all(key, "jq")) %>% #去掉关键词中的jq
  rename(year = key) %>% #将关键词命名为年份
  mutate(prov = str_sub(city, 1, 2)) %>% #Create, modify, and delete columns列
  type_convert() -> provdf1 

（1）gather（） #Gather columns into key-value pairs，将xlsx中的列读取进来

gather(
  data,
  key = "key",
  value = "value",
  ...,
  na.rm = FALSE,
  convert = FALSE,
  factor_key = FALSE
)

5.填充数据

# 机器人安装密度
readxl::read_xlsx("画图.xlsx") %>% 
  gather(2:15, key = "key", value = "value") %>% #横坐标轴数据（年份）
  mutate(key = str_remove_all(key, "jq")) %>% #去掉关键词中的jq
  rename(year = key) %>% #将关键词命名为年份
  mutate(prov = str_sub(city, 1, 2)) %>% #Create, modify, and delete columns列
  type_convert() -> provdf1 
 
provmap %>% 
  mutate(prov = str_sub(省, 1, 2)) %>% 
  left_join(provdf1) %>% 
  mutate(value = if_else(is.na(value), -1, value)) -> provmap2 
   #对应
# quantile(provmap2$value, probs = 1:8/8, digits = 2) -> cutlist#分位数，通用函数分位数产生对应于给定概率的样本分位数。最小的观察对应于概率0，最大的对应于概率1
# cutlist <- c(-1, 0, cutlist)# c,Combine Values into a Vector or List
# cutlist

 left_join()     |   keeps all observations in x.

The three outer joins keep observations that appear in at least one of the data frames:

• A left_join() keeps all observations in x.

• A right_join() keeps all observations in y.

• A full_join() keeps all observations in x and y.

inner_join(
  x,
  y,
  by = NULL,
  copy = FALSE,
  suffix = c(".x", ".y"),
  ...,
  keep = NULL
)
 
## S3 method for class 'data.frame'
inner_join(
  x,
  y,
  by = NULL,
  copy = FALSE,
  suffix = c(".x", ".y"),
  ...,
  keep = NULL,
  na_matches = c("na", "never"),
  multiple = "all",
  unmatched = "drop",
  relationship = NULL
)
 
left_join(
  x,
  y,
  by = NULL,
  copy = FALSE,
  suffix = c(".x", ".y"),
  ...,
  keep = NULL
)
 
## S3 method for class 'data.frame'
left_join(
  x,
  y,
  by = NULL,
  copy = FALSE,
  suffix = c(".x", ".y"),
  ...,
  keep = NULL,
  na_matches = c("na", "never"),
  multiple = "all",
  unmatched = "drop",
  relationship = NULL
)
 
right_join(
  x,
  y,
  by = NULL,
  copy = FALSE,
  suffix = c(".x", ".y"),
  ...,
  keep = NULL
)
 
## S3 method for class 'data.frame'
right_join(
  x,
  y,
  by = NULL,
  copy = FALSE,
  suffix = c(".x", ".y"),
  ...,
  keep = NULL,
  na_matches = c("na", "never"),
  multiple = "all",
  unmatched = "drop",
  relationship = NULL
)
 
full_join(
  x,
  y,
  by = NULL,
  copy = FALSE,
  suffix = c(".x", ".y"),
  ...,
  keep = NULL
)
 
## S3 method for class 'data.frame'
full_join(
  x,
  y,
  by = NULL,
  copy = FALSE,
  suffix = c(".x", ".y"),
  ...,
  keep = NULL,
  na_matches = c("na", "never"),
  multiple = "all",
  relationship = NULL
)

6.分组分段（图例）

# 分组
provmap2 %>% 
  mutate(group = cut(value, breaks = cutlist, 
                     include.lowest = T,
                     labels = c("No data", "<= 186", "186~402", "402~765",
                                "765~1588", "1588~3009", "3009~6230",
                                "6230~13740", "> 13740"))) -> provmap2
#去除多余名称
provmap %>% 
  mutate(label = str_remove_all(省, "回族|省|维吾尔自治区|壮族自治区|自治区|特别行政区|市|中朝共有")) -> provmap

 7.绘图|比例尺、图例、指北针、经纬网、边框、字体等

 
# 绘图
library(ggspatial) 
library(ggnewscale) 
 
provmap2 %>% 
  filter(year == 2006 | is.na(year)) %>% 
  ggplot() + 
  geom_sf(aes(fill = group), color = "gray", linewidth = 0.01) + 
  geom_sf(data = provlinemap, 
          aes(color = class, linewidth = class),
          show.legend = F) + 
  scale_fill_manual(values = c("gray", "#e0f2f1", "#b2dfdb", "#80cbc4", 
                               "#4db6ac", "#26a69a", "#009688", 
                               "#00897b", "#00796b")) + 
  stat_sf_coordinates(data = provmap,
                      geom = "text", color = "#90a4ae",
                      aes(label = label), family = cnfont,
                      fun.geometry = st_point_on_surface,
                      size = 3) + 
  scale_color_manual(
    values = c("九段线" = "#A29AC4",
               "海岸线" = "#0055AA",
               "小地图框格" = "black")
  ) + 
  scale_linewidth_manual(
    values = c("九段线" = 0.6,
               "海岸线" = 0.3,
               "小地图框格" = 0.3)
  ) + 
  annotation_scale(location = "bl", 
                   width_hint = 0.3, 
                   text_family = cnfont) + 
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        legend.position = c(0.1, 0.2),
        legend.title = element_blank(),
        plot.background = element_rect(fill = "white", color = "white")) + 
  annotation_north_arrow(
    location = "tr", 
    which_north = "false",
    pad_y = unit(0.1, "cm"), 
    style = north_arrow_nautical(
      text_family = cnfont
    ) 
  ) -> p1 
 
ggsave("pic1.png", width = 9, height = 8, device = png) 

8.循环|批量制图

 
# 循环所有年份的
dir.create("robotpic")
for (y in 2006:2019) {
  provmap2 %>% 
    filter(year == y | is.na(year)) %>% 
    ggplot() + 
    geom_sf(aes(fill = group), color = "gray", linewidth = 0.01) + 
    geom_sf(data = provlinemap, 
            aes(color = class, linewidth = class),
            show.legend = F) + 
    scale_fill_manual(values = c("gray", "#e0f2f1", "#b2dfdb", "#80cbc4", 
                                 "#4db6ac", "#26a69a", "#009688", 
                                 "#00897b", "#00796b")) + 
    stat_sf_coordinates(data = provmap,
                        geom = "text", color = "#90a4ae",
                        aes(label = label), family = cnfont,
                        fun.geometry = st_point_on_surface,
                        size = 3) + 
    scale_color_manual(
      values = c("九段线" = "#A29AC4",
                 "海岸线" = "#0055AA",
                 "小地图框格" = "black")
    ) + 
    scale_linewidth_manual(
      values = c("九段线" = 0.6,
                 "海岸线" = 0.3,
                 "小地图框格" = 0.3)
    ) + 
    annotation_scale(location = "bl", 
                     width_hint = 0.3, 
                     text_family = cnfont) + 
    guides(fill = guide_legend(ncol = 2, byrow = F)) + 
    theme(axis.title.x = element_blank(),
          axis.title.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          legend.position = c(0.135, 0.2),
          legend.title = element_blank(),
          plot.background = element_rect(fill = "white", color = "white")) + 
    annotation_north_arrow(
      location = "tr", 
      which_north = "false",
      pad_y = unit(0.1, "cm"), 
      style = north_arrow_nautical(
        text_family = cnfont
      ) 
    ) -> ptemp 
  
  ggsave(paste0("robotpic/pic", y, ".png"), width = 9, height = 7, device = png) 
  ggsave(paste0("robotpic/pic", y, ".pdf"), width = 9, height = 7, device = cairo_pdf) 
}