空间数据绘图with R

本文由SCY编写，转载请注明出处。

这篇文章致力于更新R绘制空间数据的方法。

前提条件

加载需要的包

library(terra)        # 空间数据处理
library(scico)        # 绘图
library(tidyverse)    # 数据处理
library(tidyterra)    # 空间数据制图
library(ggspatial)    # 空间数据制图
library(patchwork)    # 拼图

导入数据：（1）DEM；（2）fvc_slope、fvc_pvalue & fvc_zvalue

setwd("D:/Data")                    # 设置工作空间
DEM <- rast("dem.tif")              # 读取DEM
StudyRegion <- vect("yjq.shp")      # 读取研究区边界

绘图要求：边框，经纬度，内部轴刻度线，指北针，比例尺，图例分段显示

绘制DEM

ggplot() +
  geom_spatraster(data = DEM) +     # 绘制DEM
  geom_spatvector(
    data = StudyRegion,             # 绘制研究区边界
    fill = NA,                      # 无填充
    color = "gray50"                # 灰色边框
    ) +
  scale_fill_distiller(
    palette = "Greens",             # 调色板
    name = "DEM (m)",               # 图例标题        
    na.value = NA,                  # 无数据值不显示
    direction = 1,                  # 色带方向
    breaks = seq(0, 7000, 1000)     # 图例分段
  ) +
  guides(fill = guide_coloursteps(direction = "horizontal")) +       # 图例形状，显示方向
  annotation_north_arrow(
    height = unit(1, "cm"),         # 指北针高度
    width = unit(0.6, "cm"),        # 指北针宽度
    location = "tl") +              # 指北针位置
  annotation_scale(
    style = "ticks",                # 比例尺样式
    location = "bl",                # 比例尺位置
    width_hint = 0.5) +             # 比例尺宽度
  theme(
    legend.title.position = "top",                  # 图例标题位置（在上方）
    legend.position = "bottom",                     # 图例位置（在底部）
    legend.key.width = unit(2, "cm"),               # 图例宽度
    legend.key.height = unit(0.3, "cm"),            # 图例高度
    legend.text = element_text(vjust = 0),          # 图例文本位置（垂直居左）
    panel.grid = element_line(
      color = "gray80",                             # 网格线颜色
      linetype = 2),                                # 网格线类型（虚线）
    panel.background = element_blank(),             # 绘图区背景（空白）
    panel.border = element_rect(
      color = "black",                              # 绘图区边框颜色（黑色）
      fill = NA,                                    # 无填充
      linewidth = 0.5),                             # 边框线宽
    axis.ticks.length = unit(-0.1, "cm"),           # 轴刻度线长度（负号表示内部显示）
    axis.text = element_text(color = "black")       # 轴文本颜色
  )

成图效果：

绘制fvc_slope、fvc_pvalue & fvc_zvalue

slope <- rast("fvc_slope.tif")
pvalue <- rast("fvc_MKtest.tif")
zvalue <- rast("fvc_Zs.tif")

根据p值标记通过95%检验区域（函数）

stippling <- function(x){
  grid <- aggregate(x, fact = 20, fun = min)        # 降采样
  significant <- grid < 0.05                        # 通过95%检验
  significant_coords <- xyFromCell(significant, which(values(significant)))       # 提取坐标
  return(as.data.frame(significant_coords))         # 返回数据框
}

p1 <- ggplot() +
  geom_spatraster(data = slope) +                   # 绘制斜率
  scale_fill_scico(                      
    palette = "vik",                                # 调色板
    na.value = NA,
    limits = c(-0.03, 0.03),                        # 控制斜率范围
    breaks = seq(-0.03, 0.03, 0.01)                 # 控制斜率分段
  ) +
  guides(fill = guide_colorbar(title = "Slope")) +  # 图例形状，标题
  annotate(
    "text",                                         # 添加文本（panel 编号）
    x = 75,                                         # x坐标
    y = Inf,                                        # y坐标
    label = "(a)",                                  # 文本内容
    hjust = 1.1,                                    # 水平位置
    vjust = 1.1,                                    # 垂直位置
    size = 5                                        # 字体大小
  ) +
  geom_point(
    data = stippling(pvalue),                       # 通过95%检验区域的坐标
    aes(x = x, y = y),                              # x，y坐标
    color = "black",                                
    size = 0.1
  ) +
  theme(
    legend.title.position = "top",                  # 图例标题位置（在上方）       
    legend.title = element_text(face = "italic", margin = margin(b = 20)),    # 图例标题样式
    legend.key.width = unit(0.25, "cm"),            # 图例宽度
    legend.key.height = unit(0.8, "cm"),            # 图例高度
    legend.ticks = element_blank(),                 # 图例刻度线
    legend.text = element_text(vjust = 0),          # 图例文本位置（垂直居左）
    panel.grid = element_line(color = "gray80", linetype = 2),            # 网格线
    panel.background = element_blank(),             # 绘图区背景（空白）
    panel.border = element_rect(                    # 绘图区边框
      color = "black",                              # 边框颜色
      fill = NA,                                    # 无填充
      linewidth = 0.5                               # 边框线宽
    ),
    axis.title = element_blank(),                   # 轴标题
    axis.ticks.length = unit(-0.1, "cm"),           # 轴刻度线长度（负号表示内部显示）
    axis.text = element_text(color = "black"),      # 轴文本颜色
    axis.text.y = element_text(angle = 90, hjust = 0.5)        # y轴文本角度（垂直显示）
  )

p2 <- ggplot() +
  geom_spatraster(data = zvalue) +
  scale_fill_distiller(
    palette = "RdBu",
    na.value = NA,
    limits = c(-7.5, 7.5),
    breaks = seq(-7.5, 7.5, 2.5)
  ) +
  guides(fill = guide_colorbar(title = "Z Value")) +
  annotate(
    "text",
    x = 75,
    y = Inf,
    label = "(b)",
    hjust = 1.1,
    vjust = 1.1,
    size = 5
  ) +
  theme(
    legend.title.position = "top",
    legend.title = element_text(face = "italic", margin = margin(b = 20)),
    legend.key.width = unit(0.25, "cm"),
    legend.key.height = unit(0.8, "cm"),
    legend.ticks = element_blank(),
    legend.text = element_text(vjust = 0),
    panel.grid = element_line(color = "gray80", linetype = 2),
    panel.background = element_blank(),
    panel.border = element_rect(
      color = "black",
      fill = NA,
      linewidth = 0.5
    ),
    axis.title = element_blank(),
    axis.ticks.length = unit(-0.1, "cm"),
    axis.text = element_text(color = "black"),
    axis.text.y = element_text(angle = 90, hjust = 0.5)
  )

p1 + p2 + plot_layout(ncol = 1)      # 拼图（一列分布）

成图效果：
![fvc变化趋势](https://scy11.oss-cn-beijing.aliyuncs.com/img/202408192210342.png）