Search code examples
rr-sf

R - how to produce a distance map with the sf package?


Consider the following code:

library(sf)
library(tidyverse)

nc <- st_read(system.file("shape/nc.shp", package="sf")) %>% 
  st_union()

centroid <- st_centroid(nc)

ggplot(nc) +
  geom_sf() +
  geom_sf(data = centroid)

enter image description here

I now want to create a distance map which shows how far away the area is from the centroid. I'm looking for something that is in the style of third image from the bottom on this here link.

How can I create the map? Thanks.


Solution

  • You could create a distance raster, need to transform to projected CRS first and the relevant ggplot layer is from tidyterra. EPSG:6542 happens to be metric, so you might want to pick something else or just convert units.

    library(sf)
    library(terra)
    library(tidyterra)
    library(ggplot2)
    
    # transform to projected CRS for distance calculation, units are meters
    # NAD83(2011) / North Carolina (EPSG:6542) 
    nc <- st_read(system.file("shape/nc.shp", package="sf")) |>
      st_union() |>
      st_transform(nc, crs = 6542)
    
    centroid <- st_centroid(nc)
    
    nc_spatvect <- vect(nc)
    c_spatvect  <- vect(centroid)
    dist_rast <- distance(rast(nc_spatvect, resolution = 1000 ),c_spatvect) |> mask(nc_spatvect)
    
    # meters to km
    dist_rast <- dist_rast / 1000
    
    ggplot(nc) +
      geom_spatraster(data = dist_rast, aes(fill = lyr.1)) +
      scale_fill_viridis_c(na.value = NA, option = "plasma", direction = -1, labels = scales::label_number(suffix = "km")) +
      geom_sf(fill = NA) +
      geom_sf(data = centroid) +
      theme_bw()
    

    Or if you just need a simple distance plot, maybe drop sf, ggplot2 and tidyterra:

    
    nc_spatvect <- vect(system.file("shape/nc.shp", package="sf")) |> 
      aggregate() |> 
      project("epsg:6542")
    
    c_spatvect <- centroids(nc_spatvect)
    
    distance(rast(nc_spatvect, resolution = 1000), c_spatvect) |>
      mask(nc_spatvect) |> 
      plot(col = viridis::plasma(100, direction = -1))
    plot(c_spatvect, add = TRUE)
    

    Created on 2023-04-24 with reprex v2.0.2