geospacial state map, remove areas beyond state border

Question

I have the following code which produces a hydrological map of Vermont, but I want everything outside of the state to be white, so that the hydrological data in bordering states don't show up and the Vermont border is clearer. Is there a way to filter the hydological data down to just Vermont or overlay white on everything that isn't within the geometry data for Vermont?

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)

riv_list <- lapply(filenames, st_read)

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()

test<- na_riv_width %>% filter_state("South")

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)  +
    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()

A state identifier isn't in the data. And the geometries span multiple states. So I am looking for a solution as to how to clip bordering states out of the data. — Quinterpret, Mar 04 '23 at 19:12
I've added the following, but get the following error: counties <- counties(state="VT",cb=T,class="sf") %>% sf::st_cast("MULTILINESTRING") st_crop(na_riv, counties) ; Error in geos_op2_geom("intersection", x, y, ...) : st_crs(x) == st_crs(y) is not TRUE — Quinterpret, Mar 04 '23 at 19:27
You need to have the same projection. See here: https://stackoverflow.com/questions/32233576/how-to-clip-a-polygon-shapefile-by-another-polygon-shapefile-in-r — M--, Mar 04 '23 at 19:35
The `st_crop` did limit the data, but I am still seeing rivers branch into bordering states — Quinterpret, Mar 04 '23 at 20:04

score 2 · Accepted Answer · answered Mar 04 '23 at 20:11

You can get the outline of Vermont as an sf object and use it to trim the rivers using st_intersection

vermont <- rnaturalearth::ne_states(country = "United States of America", 
              returnclass = "sf") %>%
           filter(name == "Vermont")

vermont_rivers <- na_riv_width %>% st_intersection(vermont)

ggplot(vermont) +
    geom_sf() +
    geom_sf(data = vermont_rivers,
            aes(color = factor(ORD_FLOW), size = width,
                alpha = factor(ORD_FLOW))) +
    scale_color_manual(NULL, 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_void() +
    theme(legend.position = "none")

geospacial state map, remove areas beyond state border

1 Answers1