使用 R 语言绘制填充地图、栅格地图 + 等降水量线

最近有个小伙伴问到了关于在栅格地图上绘制等高线的问题。今天我们就以降水量栅格数据为例讲解如何在 R 语言中绘制这种地图。

在学习本课程之前，需要预先学习之前的课程：

使用 R 语言绘制历年中国省市区县地图（小地图版本+长版）：https://rstata.duanshu.com/#/brief/course/379d19770956478ebc4d919910fca74c

附件中提供了 2021 年各月的降水量栅格数据：pre_2021.nc，数据来源于 中国1km分辨率逐月降水量数据集（1901-2023）：https://data.tpdc.ac.cn/zh-hans/data/faae7605-a0f2-4d18-b28f-5cee413766a2 。

填充地图

首先我们来绘制各省份的年均降水量填充地图。

读取栅格数据和省级矢量数据：

library(tidyverse)
library(sf)
library(terra)

# 读取降水量栅格数据
rast("pre_2021.nc") -> rst
rst

#> class       : SpatRaster
#> dimensions  : 5146, 7849, 12  (nrow, ncol, nlyr)
#> resolution  : 0.008333333, 0.008333333  (x, y)
#> extent      : 71.28589, 136.6942, 15.75212, 58.63546  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84
#> source      : pre_2021.nc
#> varname     : pre (monthly precipitation)
#> names       : pre_1, pre_2, pre_3, pre_4, pre_5, pre_6, ...

# 年均
app(rst, mean, na.rm = T) -> rst
rst

#> class       : SpatRaster
#> dimensions  : 5146, 7849, 1  (nrow, ncol, nlyr)
#> resolution  : 0.008333333, 0.008333333  (x, y)
#> extent      : 71.28589, 136.6942, 15.75212, 58.63546  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84
#> source(s)   : memory
#> name        :        mean
#> min value   :    4.166667
#> max value   : 6252.083333

# 提取中国范围的
read_sf("2021行政区划/省.shp") %>%
  st_union() %>%
  nngeo::st_remove_holes() %>%
  st_sfc() -> cn

plot(cn)

rst %>%
  terra::crop(vect(cn)) %>%
  terra::mask(vect(cn)) -> rst

# 读取省级行政区划
read_sf("2021行政区划/省.shp") %>%
  select(-contains("类型")) -> prov
prov

#> Simple feature collection with 34 features and 2 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 73.50235 ymin: 3.83703 xmax: 135.0957 ymax: 53.56362
#> Geodetic CRS:  WGS 84
#> # A tibble: 34 × 3
#>    省             省代码                                                geometry
#>    <chr>           <dbl>                                      <MULTIPOLYGON [°]>
#>  1 安徽省         340000 (((116.426 34.64413, 116.4263 34.64, 116.4261 34.63816…
#>  2 澳门特别行政区 820000 (((113.5815 22.19861, 113.5812 22.19779, 113.5789 22.1…
#>  3 北京市         110000 (((116.6255 41.05897, 116.6258 41.05887, 116.6264 41.0…
#>  4 福建省         350000 (((117.689 23.53157, 117.6888 23.53144, 117.6875 23.53…
#>  5 甘肃省         620000 (((106.0718 35.44995, 106.0713 35.44955, 106.0707 35.4…
#>  6 广东省         440000 (((110.5657 20.33342, 110.5604 20.32783, 110.5565 20.3…
#>  7 广西壮族自治区 450000 (((109.2167 20.90438, 109.2144 20.90184, 109.2108 20.9…
#>  8 贵州省         520000 (((105.0948 24.92524, 105.0947 24.9252, 105.0946 24.92…
#>  9 海南省         460000 (((112.0522 3.838239, 112.0483 3.837329, 112.0181 3.83…
#> 10 河北省         130000 (((118.2715 38.98381, 118.2715 38.982, 118.2707 38.979…
#> # ℹ 24 more rows

分区域计算平均降水量：

terra::extract(rst, vect(prov), fun = mean, na.rm = T) %>%
  as_tibble() %>%
  bind_cols(st_drop_geometry(prov)) %>%
  select(-ID) %>%
  rename(年度降水量_mm = mean) -> df
df

#> # A tibble: 34 × 3
#>    年度降水量_mm 省             省代码
#>            <dbl> <chr>           <dbl>
#>  1         1119. 安徽省         340000
#>  2         1597. 澳门特别行政区 820000
#>  3          519. 北京市         110000
#>  4         1299. 福建省         350000
#>  5          286. 甘肃省         620000
#>  6         1271. 广东省         440000
#>  7         1227. 广西壮族自治区 450000
#>  8         1071. 贵州省         520000
#>  9         1279. 海南省         460000
#> 10          485. 河北省         130000
#> # ℹ 24 more rows

然后就可以绘制填充地图了，这部分代码类似之前的课程：

# 绘制填充地图
# 中国地图通常使用这样的坐标系
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"

# 读取小地图版本的中国省级地图
read_sf("chinaprov2019mini/chinaprov2019mini.shp") %>%
  filter(!is.na(省代码)) -> provmap

# 线条
read_sf("chinaprov2019mini/chinaprov2019mini_line.shp") %>%
  filter(class %in% c("九段线", "海岸线", "小地图框格")) %>%
  select(class) -> provlinemap

# 小地图的范围
small_bbox <- st_bbox(c(xmin = 120000,
                        xmax = 1766004.1,
                        ymax = 2557786.0,
                        ymin = 320000),
                      crs = st_crs(mycrs)) %>%
  st_as_sfc()

# 绘图
library(ggspatial)
library(ggnewscale)

# 绘制填充地图

# 和地图数据合并
provmap %>%
  left_join(df) -> provmap2

hist(provmap2$年度降水量_mm)

# 绘图
provmap2 %>%
  mutate(年度降水量_mm = if_else(is.na(年度降水量_mm), -400,
                             年度降水量_mm),
         group = cut(年度降水量_mm,
                     breaks = seq(-400, 2400, by = 400), include.lowest = T,
                     labels = c("No data", "0 ~ 400", "400 ~ 800", "800 ~ 1200",
                                "1200 ~ 1600", "1600 ~ 2000", "2000 ~ 2400"))) %>%
  # st_drop_geometry() %>%
  # count(group)
  ggplot() +
  geom_sf(aes(fill = group),
          color = "gray", linewidth = 0.1) +
  geom_sf(data = provlinemap,
          aes(color = class, linewidth = class),
          show.legend = F) +
  scale_color_manual(
    values = c("九段线" = "#A29AC4",
               "海岸线" = "#0055AA",
               "小地图框格" = "black",
               "省份" = "black")
  ) +
  scale_linewidth_manual(
    values = c("九段线" = 0.6,
               "海岸线" = 0.3,
               "小地图框格" = 0.3,
               "省份" = 0.3)
  ) +
  scale_fill_manual(values = c("white", c("#d53e4f", "#fc8d59",
                                          "#fee08b", "#e6f598",
                                          "#99d594", "#3288bd"))) +
  new_scale_color() +
  stat_sf_coordinates(data = provmap,
                      geom = "text", color = "gray",
                      aes(label = 省), family = cnfont,
                      fun.geometry = st_point_on_surface,
                      size = 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(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        legend.position = c(0.12, 0.2),
        plot.background = element_rect(fill = "white",
                                       color = "white")) +
  guides(fill = guide_legend(ncol = 2)) +
  labs(fill = "Precipitation\n(thounds mm)") +
  annotation_north_arrow(
    location = "tr",
    which_north = "false",
    pad_y = unit(0.1, "cm"),
    style = north_arrow_fancy_orienteering(
      text_family = cnfont
    )
  ) -> p1

ggsave("pic1.png", width = 10, height = 8.5, device = png)

栅格地图

栅格地图的绘制方法也和之前课程中讲解的一样。栅格数据绘制地图需要转换成点或者多边形，所以这里我们改用 raster 包处理栅格数据：

# 栅格数据
library(raster)
raster(rst) -> rst

# 计算等降水量线
rasterToContour(rst, levels = seq(0, 2400, by = 400)) %>%
  st_as_sf() -> line
line

#> Simple feature collection with 6 features and 1 field
#> Geometry type: MULTILINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 78.54866 ymin: 18.43559 xmax: 132.7324 ymax: 49.96448
#> Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs
#>     level                       geometry
#> C_1   400 MULTILINESTRING ((78.70524 ...
#> C_2   800 MULTILINESTRING ((79.12219 ...
#> C_3  1200 MULTILINESTRING ((89.0151 2...
#> C_4  1600 MULTILINESTRING ((92.68309 ...
#> C_5  2000 MULTILINESTRING ((117.7477 ...
#> C_6  2400 MULTILINESTRING ((120.5423 ...

plot(line)

# 栅格数据转点（高分辨率）或者多边形（低分辨率），这里使用多边形试试：
rst %>%
  aggregate(fact = 10) %>% # 降低分辨率，提升绘图速度
  rasterToPolygons() %>%
  st_as_sf() -> rstpoints

rstpoints %>%
  st_transform(mycrs) -> rstpoints

# 提取小地图内的点
rstpoints %>%
  st_intersection(small_bbox) -> rstpoints_small

# 主图中的点
st_bbox(c(xmin = -2725586,
          xmax = 2982768,
          ymax = 6000000,
          ymin = 1800655),
        crs = st_crs(mycrs)) %>%
  st_as_sfc() -> plotbbox
rstpoints %>%
  st_intersection(plotbbox) -> rstpoints_main

# 移动小地图内的点到恰当位置：
rstpoints_small %>%
  mutate(geometry = geometry * 0.5 + c(2100000, 1665139)) %>%
  sf::st_set_crs(mycrs) -> rstpoints_small

# 合并两个部分
bind_rows(rstpoints_main, rstpoints_small) -> rstpoints_all

line %>%
  st_transform(mycrs) -> line

# 绘图
ggplot(provmap) +
  geom_sf(data = rstpoints_all, aes(fill = mean),
          linewidth = 0.001,
          color = NA) +
  scale_fill_gradientn(colours = c("#9e0142", "#d53e4f",
                                   "#f46d43", "#fdae61",
                                   "#fee08b", "#ffffbf",
                                   "#e6f598", "#abdda4",
                                   "#66c2a5", "#3288bd",
                                   "#5e4fa2"),
                       name = "Precipitation\n(thounds mm)") +
  new_scale_color() +
  geom_sf(fill = NA, color = "gray",
          linewidth = 0.03) +
  geom_sf(data = provlinemap,
          aes(color = class, linewidth = class),
          show.legend = F) +
  geom_sf(data = line, color = "white", linewidth = 0.5) +
  stat_sf_coordinates(data = provmap,
                      geom = "text", color = "gray",
                      aes(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(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        legend.position = c(0.09, 0.2),
        legend.direction = "horizontal",
        legend.title.position = "top") +
  annotation_north_arrow(
    location = "tr",
    which_north = "false",
    pad_y = unit(0.1, "cm"),
    style = north_arrow_fancy_orienteering(
      text_family = cnfont
    )
  ) -> p2

ggsave("pic2.png", width = 10, height = 8.5, device = png)

这样我们就把等降水线添加到了地图上。

添加等降水线文本标签

如果想把等降水量线的文本标签标注到线上，可以借助 geomtextpath 包。

安装：

remotes::install_github("AllanCameron/geomtextpath")

然后就可以绘制了：

library(geomtextpath)
line %>%
  st_cast("LINESTRING") %>%
  mutate(level = paste0(level, " mm")) -> line2

ggplot(provmap) +
  geom_sf(data = rstpoints_all, aes(fill = mean),
          linewidth = 0.001,
          color = NA) +
  scale_fill_gradientn(colours = c("#9e0142", "#d53e4f",
                                   "#f46d43", "#fdae61",
                                   "#fee08b", "#ffffbf",
                                   "#e6f598", "#abdda4",
                                   "#66c2a5", "#3288bd",
                                   "#5e4fa2"),
                       name = "Precipitation\n(thounds mm)") +
  new_scale_color() +
  geom_sf(fill = NA, color = "gray",
          linewidth = 0.03) +
  geom_sf(data = provlinemap,
          aes(color = class, linewidth = class),
          show.legend = F) +
  geom_textsf(data = line2, aes(label = level),
              family = cnfont,
              linecolour = "white", color = "gray30",
              text_smoothing = 50, size = 3) +
  stat_sf_coordinates(data = provmap,
                      geom = "text", color = "gray",
                      aes(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(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        legend.position = c(0.09, 0.2),
        legend.direction = "horizontal",
        legend.title.position = "top") +
  annotation_north_arrow(
    location = "tr",
    which_north = "false",
    pad_y = unit(0.1, "cm"),
    style = north_arrow_fancy_orienteering(
      text_family = cnfont
    )
  ) -> p3

ggsave("pic3.png", width = 10, height = 8.5, device = png)

点击这里跳转到 RStata 短书平台获取附件：使用 R 语言绘制填充地图、栅格地图 + 等降水量线