

ggplot2优雅的带你绘制中国地图

R语言数据分析指南

2021-12-22

欢迎关注"R数据分析指南"

本节来介绍如何使用sf包来绘制中国地图，本文由阿猹同学完成，我们提供数据可视化的有偿服务，有需要的小伙伴欢迎咨询

加载R包

rm(list = ls())
pacman::p_load(tidyverse, sf, raster, ggspatial, stars)

设置投影

crs_84 <- st_crs("EPSG:4326")  ## WGS 84 大地坐标
crs_al <- st_crs("+proj=aea +lat_1=25 +lat_2=47 +lon_0=105") ## Albers Equal Area Conic投影

地图获取

http://datav.aliyun.com/portal/school/atlas/area_selector#&lat=31.769817845138945&lng=104.29901249999999&zoom=4

china_all <-
    sf::st_read("https://geo.datav.aliyun.com/areas_v3/bound/100000_full.json") %>%
    st_transform(crs_al)
hainan <-
    sf::st_read("https://geo.datav.aliyun.com/areas_v3/bound/460000_full.json") %>%
    st_transform(crs_al)

截取地图

tmp_china <-  # 去除 海南省，九段线
    china_all %>%
    filter(!adcode %in% c("460000", "100000_JD")) %>%
    st_make_valid() %>%
    st_union()

tmp_hainan <-  # 海南省去除 三沙市
    hainan %>%
    filter(!name %in% "三沙市") %>%
    st_make_valid() %>%
    st_union()

组合起来中国大陆边框

china_com <- st_union(tmp_china, tmp_hainan)%>% st_as_sf()

海拔获取

dem <- geoviz::mapbox_dem(
    lat = 35.8617,
    long = 104.1954,
    square_km = 2000,
    api = "pk.eyJ1IjoiYmVueXNmIiwiYSI6ImNrczBtdWE0ajBwNjcydnBqMjRyZDdsOXkifQ.sUcMdooE7b9uQqzfrnWdSQ"
)

china_dem <- dem %>%
    projectRaster(crs = crs_al$wkt) %>% # 修改栅格数据的投影
    aggregate(fact = 3) %>%  ## 降低分辨率，减少运算量
    raster::crop(., raster::extent(china_com)) %>%
    raster::mask(china_com) %>%
    stars::st_as_stars() %>%
    st_as_sf()

颜色配置

colors <- c(
    "#33A02C", "#B2DF8A", "#FDBF6F", "#1F78B4",
    "#999999", "#E31A1C", "#E6E6E6", "#A6CEE3")

绘制大陆区域

p1 <-
    ggplot() +
    geom_sf(aes(fill = layer,  color = layer), data = china_dem) +
    geom_sf(size = .2, fill = "transparent", color = "#060d1b", data = china_all) +
    scale_fill_gradientn(colours = colors) +
    scale_color_gradientn(colours = colors)

截取南海

p2 <-
    p1 +
    coord_sf(crs = crs_84) + ## 将投影坐标转换为大地坐标
    scale_x_continuous(expand = c(0, 0), limits = c(107, 122), breaks = seq(70, 140, 10)) +
    scale_y_continuous(expand = c(0, 0), limits = c(2, 24), breaks = seq(10, 60, 10)) +
    guides(fill = "none", color = "none") +
    theme_bw() +
    theme(
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank()
    )

结合采样点绘图

# 采样点随机生成
set.seed(123)
sites <- data.frame(
    lng = runif(10, 90, 120),
    lat = runif(10, 30, 40),
    lab = letters[1:10])

图形拼接

p1 + coord_sf(crs = crs_al, default_crs = crs_84) +
    geom_point(data = sites, aes(x = lng, y = lat), color = "red",size=2) +  ## 采样点图层
    ggrepel::geom_label_repel(data = sites, aes(x = lng,y=lat,label=lab), family = "serif") + ## 加标签
    annotate(geom = "text", x = 80, y = 18,label="GS(2019)6379",family = "serif", vjust = 0, hjust = 0) +  ##审图号
    scale_x_continuous(expand = c(0, 0),limits=c(72,142),breaks=seq(70, 140, 10)) +
    scale_y_continuous(expand = c(0, 0),limits = c(17,55.5), breaks = seq(10, 60, 10)) +
    labs(fill = "elevant", color = "elevant") +
    theme_bw() +
    theme(
        axis.text = element_text(family ="serif",color="black"),  ## 字体改为新罗马
        axis.title = element_blank(),
        legend.position = c(1,0.8),legend.justification = c(1,1)) +
    annotation_scale(location = "bl") + # 设置距离刻度尺
    annotation_north_arrow(location = "tl", style = north_arrow_nautical(
        fill = c("grey40", "white"), line_col = "grey20")) +
    annotation_custom(ggplotGrob(p2),xmin= 122,xmax = 138,ymin=15,ymax = 29)