Search code examples
rggplot2terratidyterra

Using coord_sf with geom_spatraster changes the resolution

I am trying to plot a subregion of a large raster (e.g., just Arizona out of a raster that covers North America) using tidyterra and ggplot2. I have tried setting the limits of the plot by setting the CRS to lon/lat WGS 84 using coord_sf and then using lim with the limits in lon/lat. When I do this, the resulting plot shows a much lower resolution than what the raster originally has (1 km). I realize I could set the limits in the original CRS of the raster, but I would also like to plot vectors on top in lat/lon WGS 84, rather than the original CRS. Is there a way that I can retain the original resolution while also using lat/long WGS 84 to specify the limits and for additional geoms?

library(tidyverse)
library(ggplot2)
library(terra)
library(tidyterra)
library(maps)

### Example data:

# example raster:
ex <- rast(
  extent = c(-4352000, 3336000, -3341000, 4276000),
  resolution = c(1000, 1000),
  crs = "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
)
ex[] <- runif(n = ncell(ex), min = 0, max = 100)

# example polygon
az <- map_data("state") %>%
  filter(region %in% c("arizona")) %>%
  vect(
    x = cbind(
      as.factor(.$region),    # object (state) id
      .$group,                # polygon part id
      .$long,                 # x
      .$lat,                  # y
      0                             # hole flag (all 0)
    ),
    type = "polygons",
    crs = "+proj=longlat +datum=WGS84"
  )

### Plot:

p <- ggplot() +
  geom_spatraster(data = ex) +
  geom_spatvector(data = az, fill = NA, color = "white") +
  coord_sf(crs = 4326) +
  lims(x = c(-116, -108), y = c(30, 38))
print(p)

enter image description here

Solution

  • To avoid resampling issues, you can terra::crop() your data prior to plotting. Given you are wanting to plot your data in WGS84/EPSG:4326, you have to project your spatraster prior to cropping. Doing so negates the need to define a CRS in coords_sf() as all data already have the same CRS. I have provided two options for defining how to crop your data.

    library(terra)
library(tidyterra)
library(sf)
library(maps)
library(ggplot2)

# Example raster
set.seed(1)
ex <- rast(
  extent = c(-4352000, 3336000, -3341000, 4276000),
  resolution = c(1000, 1000),
  crs = "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
) %>%
  project("EPSG:4326")

ex[] <- runif(n = ncell(ex), min = 0, max = 100)

# Create sf for cropping raster using your extent values
crop_sf <- data.frame(lon = c(-116, -116, -108, -108),
                      lat = c(30, 38, 38, 30)) %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

crop_sf <- do.call(c, st_geometry(crop_sf)) %>%
  st_cast("POLYGON") %>%
  st_sfc() %>%
  st_as_sf(crs = 4326)

# Example polygon
az <- map_data("state") %>%
  filter(region %in% c("arizona")) %>%
  vect(
    x = cbind(
      as.factor(.$region),    # object (state) id
      .$group,                # polygon part id
      .$long,                 # x
      .$lat,                  # y
      0                       # hole flag (all 0)
    ),
    type = "polygons",
    crs = "EPSG:4326"
  )

    Crop ex spatraster using crop_sf:

    ex1 <- crop(ex, crop_sf)

ggplot() +
  geom_spatraster(data = ex1) +
  geom_spatvector(data = az, fill = NA, color = "white") +
  coord_sf(xlim = c(-116, -108), ylim = c(30, 38))

    result

    Crop ex spatraster using az spatvector:

    ex2 <- crop(ex, az, mask = TRUE)

ggplot() +
  geom_spatraster(data = ex2) + # , maxcell = 1e6
  geom_spatvector(data = az, fill = NA, color = "white") +
  scale_fill_continuous(na.value = "transparent") +
  coord_sf(xlim = c(-116, -108), ylim = c(30, 38))

    resul2