如何在R中优化佛蒙特州水文地图，隐藏州界外数据？

方法一：过滤佛蒙特州范围内的水文数据 #

直接保留州界内的河流数据，从根源上避免显示邻州内容，是更高效的方案。

关键步骤： #

• 获取佛蒙特州的边界矢量数据，并转换为与河流数据一致的坐标系
• 通过空间交集操作，只保留完全或部分位于州界内的河段

修改后的完整代码： #

library(osmdata)
library(ggplot2)
library(tidyverse)
require(sf)
require(tigris)
library(httr)

bbx <- getbb("VT")

url <- "https://data.hydrosheds.org/file/HydroRIVERS/HydroRIVERS_v10_na_shp.zip"

res <- GET(url,
           write_disk("na_rivers.zip"),
           progress())
unzip("na_rivers.zip")
filenames <- list.files("HydroRIVERS_v10_na_shp", pattern="*.shp", full.names=T)

list_riv <- lapply(filenames, sf::st_read)
na_riv <- list_riv[[1]] %>%
  sf::st_cast("MULTILINESTRING")

na_riv_width <- na_riv %>%
  dplyr::mutate(
    width = as.numeric(ORD_FLOW),
    width = dplyr::case_when(
      width == 3 ~ 1,
      width == 4 ~ 0.8,
      width == 5 ~ 0.6,
      width == 6 ~ 0.4,
      width == 7 ~ 0.2,
      width == 8 ~ 0.2,
      width == 9 ~ 0.1,
      width == 10 ~ 0.1,
      TRUE ~ 0
    )
  ) %>%
  sf::st_as_sf()

# 获取佛蒙特州边界并转换坐标系
vt_boundary <- states(cb = TRUE) %>% 
  filter(NAME == "Vermont") %>% 
  st_transform(st_crs(na_riv_width))

# 过滤出佛蒙特州范围内的河流
vt_riv_width <- na_riv_width %>% 
  st_intersection(vt_boundary)

get_river_map <- function() {
  p <- 
    ggplot() +
    geom_sf(
      data = vt_riv_width,  # 使用过滤后的河流数据
      aes(
        color = factor(ORD_FLOW), size = width,
        alpha = factor(ORD_FLOW)
      )
    ) +
    coord_sf(xlim = c(min(bbx[1,]), max(bbx[1,])),
             ylim = c(min(bbx[2,]), max(bbx[2,])),
             expand = FALSE)  +
    labs(
      y = "", subtitle = "",
      x = "",
      title = "",
      caption = ""
    ) +
    scale_color_manual(
      name = "",
      values = c(
        "#08306b", "#08519c", "#2171b5",
                 "#4292c6", "#6baed6", "#9ecae1",
                 "#c6dbef", "#deebf7", "white"
      )
    ) +
    scale_size(range = c(0, .3)) +
    scale_alpha_manual(values = c(
      "3" = 1, "4" = 1, "5" = .7, "6" = .6,
      "7" = .4, "8" = .3, "9" = .2, "10" = .1
    )) +
    theme_minimal() +
    theme(
      panel.background = element_blank(),
      legend.background = element_blank(),
      legend.position = "none",
      panel.border = element_blank(),
      panel.grid.minor = element_blank(),
      panel.grid.major = element_blank(),
      plot.title = element_text(
        size = 40, color = "#2171b5", hjust = 0.5, vjust = 0
      ),
      plot.subtitle = element_text(
        size = 14, color = "#ac63a0", hjust = 0.5, vjust = 0
      ),
      plot.caption = element_text(
        size = 10, color = "grey60", hjust = 0.5, vjust = 10
      ),
      axis.title.x = element_text(
        size = 10, color = "grey20", hjust = 0.5, vjust = -6
      ),
      legend.text = element_text(
        size = 9, color = "grey20"
      ),
      legend.title = element_text(size = 10, color = "grey20"),
      strip.text = element_text(size = 12),
      plot.margin = unit(c(t = 1, r = -2, b = -1, l = -2), "lines"),
      axis.title.y = element_blank(),
      axis.ticks = element_blank(),
      axis.text.x = element_blank(),
      axis.text.y = element_blank()
    )
  
  return(p)
}

p1 <- get_river_map()

方法二：叠加白色遮罩层 #

如果不想修改原始河流数据集，可以在地图最上层添加一个白色多边形，覆盖佛蒙特州边界外的区域，实现“遮罩”效果。

关键步骤： #

• 创建覆盖当前绘图范围的矩形多边形
• 用佛蒙特州边界“挖空”该矩形，得到边界外的区域
• 将这个区域填充白色，放在所有图层的最上方

修改后的完整代码： #

library(osmdata)
library(ggplot2)
library(tidyverse)
require(sf)
require(tigris)
library(httr)

bbx <- getbb("VT")

url <- "https://data.hydrosheds.org/file/HydroRIVERS/HydroRIVERS_v10_na_shp.zip"

res <- GET(url,
           write_disk("na_rivers.zip"),
           progress())
unzip("na_rivers.zip")
filenames <- list.files("HydroRIVERS_v10_na_shp", pattern="*.shp", full.names=T)

list_riv <- lapply(filenames, sf::st_read)
na_riv <- list_riv[[1]] %>%
  sf::st_cast("MULTILINESTRING")

na_riv_width <- na_riv %>%
  dplyr::mutate(
    width = as.numeric(ORD_FLOW),
    width = dplyr::case_when(
      width == 3 ~ 1,
      width == 4 ~ 0.8,
      width == 5 ~ 0.6,
      width == 6 ~ 0.4,
      width == 7 ~ 0.2,
      width == 8 ~ 0.2,
      width == 9 ~ 0.1,
      width == 10 ~ 0.1,
      TRUE ~ 0
    )
  ) %>%
  sf::st_as_sf()

# 获取佛蒙特州边界
vt_boundary <- states(cb = TRUE) %>% 
  filter(NAME == "Vermont") %>% 
  st_transform(st_crs(na_riv_width))

# 创建绘图范围的矩形多边形
bbox_poly <- st_bbox(c(xmin = min(bbx[1,]), xmax = max(bbx[1,]),
                       ymin = min(bbx[2,]), ymax = max(bbx[2,]))) %>% 
  st_as_sfc() %>% 
  st_transform(st_crs(vt_boundary))

# 生成遮罩：矩形减去佛蒙特州区域
mask <- st_difference(bbox_poly, vt_boundary)

get_river_map <- function() {
  p <- 
    ggplot() +
    geom_sf(
      data = na_riv_width,
      aes(
        color = factor(ORD_FLOW), size = width,
        alpha = factor(ORD_FLOW)
      )
    ) +
    coord_sf(xlim = c(min(bbx[1,]), max(bbx[1,])),
             ylim = c(min(bbx[2,]), max(bbx[2,])),
             expand = FALSE)  +
    # 添加白色遮罩层，放在最上层
    geom_sf(data = mask, fill = "white", color = NA) +
    labs(
      y = "", subtitle = "",
      x = "",
      title = "",
      caption = ""
    ) +
    scale_color_manual(
      name = "",
      values = c(
        "#08306b", "#08519c", "#2171b5",
                 "#4292c6", "#6baed6", "#9ecae1",
                 "#c6dbef", "#deebf7", "white"
      )
    ) +
    scale_size(range = c(0, .3)) +
    scale_alpha_manual(values = c(
      "3" = 1, "4" = 1, "5" = .7, "6" = .6,
      "7" = .4, "8" = .3, "9" = .2, "10" = .1
    )) +
    theme_minimal() +
    theme(
      panel.background = element_blank(),
      legend.background = element_blank(),
      legend.position = "none",
      panel.border = element_blank(),
      panel.grid.minor = element_blank(),
      panel.grid.major = element_blank(),
      plot.title = element_text(
        size = 40, color = "#2171b5", hjust = 0.5, vjust = 0
      ),
      plot.subtitle = element_text(
        size = 14, color = "#ac63a0", hjust = 0.5, vjust = 0
      ),
      plot.caption = element_text(
        size = 10, color = "grey60", hjust = 0.5, vjust = 10
      ),
      axis.title.x = element_text(
        size = 10, color = "grey20", hjust = 0.5, vjust = -6
      ),
      legend.text = element_text(
        size = 9, color = "grey20"
      ),
      legend.title = element_text(size = 10, color = "grey20"),
      strip.text = element_text(size = 12),
      plot.margin = unit(c(t = 1, r = -2, b = -1, l = -2), "lines"),
      axis.title.y = element_blank(),
      axis.ticks = element_blank(),
      axis.text.x = element_blank(),
      axis.text.y = element_blank()
    )
  
  return(p)
}

p1 <- get_river_map()

内容的提问来源于stack exchange，提问作者Quinterpret