Skip to content

Instantly share code, notes, and snippets.

@andrewheiss
Created December 5, 2024 19:30
Show Gist options
  • Save andrewheiss/c8ce0c0989b5a602d70c167c5422cd28 to your computer and use it in GitHub Desktop.
Save andrewheiss/c8ce0c0989b5a602d70c167c5422cd28 to your computer and use it in GitHub Desktop.
library(tidyverse)
library(sf)
#> Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(rnaturalearth)

# All countries
countries <- ne_countries(scale = 50) |>
  filter(iso_a3 != "ATA")

# German states
states_germany <- ne_states(country = "germany")

# Make a little window for zooming
europe_window <- st_sfc(
  st_point(c(3, 46)),  # left (west), bottom (south)
  st_point(c(20, 58)),  # right (east), top (north)
  crs = st_crs("EPSG:4326")   # WGS 84
) %>% 
  st_transform(crs = st_crs("EPSG:3035")) %>%  # ETRS89-extended / LAEA Europe
  st_coordinates()

# This has a double border with German states on the exterior of Germany
ggplot() +
  geom_sf(data = countries, linewidth = 0.5, color = "red") +
  geom_sf(data = states_germany, linewidth = 0.25, color = "blue", fill = "white") +
  coord_sf(
    crs = st_crs("EPSG:3035"),
    xlim = europe_window[, "X"],
    ylim = europe_window[, "Y"],
    expand = FALSE
  ) +
  theme_void()

# Extract the intersections of provinces so that we don't plot the outside
# borders and distort the coastlines
interior_state_borders <- st_intersection(states_germany) |>
  filter(n.overlaps > 1) |> 
  # Remove weird points that st_intersection() adds
  filter(!(st_geometry_type(geometry) %in% c("POINT", "MULTIPOINT")))
#> although coordinates are longitude/latitude, st_intersection assumes that they
#> are planar

# It works!
ggplot() +
  geom_sf(data = states_germany) +
  geom_sf(data = interior_state_borders, color = "red")

# fill no longer works though, since those are just borders, but we can include 
# a border-less version of the German states
ggplot() +
  geom_sf(data = countries, linewidth = 0.5, color = "red") +
  geom_sf(data = states_germany, linewidth = 0, fill = "white") +
  geom_sf(data = interior_state_borders, linewidth = 0.25, color = "blue") +
  coord_sf(
    crs = st_crs("EPSG:3035"),
    xlim = europe_window[, "X"],
    ylim = europe_window[, "Y"],
    expand = FALSE
  ) +
  theme_void()

Created on 2024-12-05 with reprex v2.1.1

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment