rggplot2mapsrastertidyterra

plot basemap in ggplot2 (basemaps_gglayer() not working)

i'm trying to add a basemap to some of my figures. from reading online, it seems like the basemaps package is the best way to go, but i haven't been able to get it to work, and have come across several issues.

here's the original code for the figure, which runs perfectly fine:

ggplot() +
  geom_raster(subset(results, !is.na(var)), mapping = aes(long, lat, fill = var), interpolate = TRUE) +
  coord_sf(datum = "ESRI:102001") +
  scale_fill_gradientn('Variance', colours = colors) + #limits = c(0, 1)) +
  geom_sf(data = lg.parks, fill = NA, color = "black") + 
  geom_sf(data = canada, fill = NA, color = "black") + 
  annotation_scale(location = "bl", width_hint = 0.5) +
  annotation_north_arrow(location = "bl", which_north = "true", 
                         pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
                         style = north_arrow_fancy_orienteering) +
  theme_void() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank())

here is what my figure looks like now:

a map of canada with a legend reading 'variance', which shows values from 0 to 0.2

i would like to add a basemap to this figure.

when i try adding the line basemap_gglayer(ext = ext(canada)) + , it returns the error

Error in st_transform.sfc(st_as_sfc(st_bbox(x)), crs = crs_webmerc) : 
  cannot transform sfc object with missing crs

i seem to be able to work around this by specifying ext = st_bbox(canada). however, if i try plotting my figure with the line basemap_gglayer(ext = st_bbox(canada)), i get another error:

Error in `scale_fill_gradientn()`:
! Discrete values supplied to continuous scale.
ℹ Example values: "#91DED8", "#91DED8", "#91DED8", "#91DED8", and "#91DED8"

so that isn't working. i tried then saving the basemap as a separate raster, then plotting it using the geom_spatraster function from the tidyterra package.

bm <- basemaps::basemap_raster(ext = st_bbox(canada))
bm <- as(bm, "SpatRaster")

ggplot() +
  geom_spatraster(data = bm) +
  geom_raster(subset(results, !is.na(var)), mapping = aes(long, lat, fill = var), interpolate = TRUE) +
  coord_sf(datum = "ESRI:102001") +
  scale_fill_gradientn('Variance', colours = colors) + #limits = c(0, 1)) +
  geom_sf(data = lg.parks, fill = NA, color = "black") + 
  geom_sf(data = canada, fill = NA, color = "black") + 
  annotation_scale(location = "bl", width_hint = 0.5) +
  annotation_north_arrow(location = "bl", which_north = "true", 
                         pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
                         style = north_arrow_fancy_orienteering) +
  theme_void() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank())

and here is the output:

a dark blue map of canada with park shapefiles on top of it

this is really not what i want it to do, lol. if i plot only the base layer, it looks like this:

map of canada in white with surrounding ocean in various shades of teal

ideally, my final figure looks like ^ this, with my original figure just layered on top.

i have my defaults set to basemaps::set_defaults(map_service = "esri", map_type = "world_terrain_base"), since i want to avoid the google/stadia maps which require an API key. i'm running R 4.3.3 on linux ubuntu v 22.04.4.

EDIT: here is what my results dataframe looks like and the raster of the dataframe:

> head(results)
     lat     long  park ecozone        var    mean_res      mean       cv
1 69.425 -125.025 FALSE     ARC 0.02040808 -0.15004042 0.2764668 0.682502
2 72.125 -124.125  TRUE     ARC 0.03036255 -0.07685251 0.1612495 1.433145
3 73.025 -123.675  TRUE     ARC 0.02992938 -0.07491194 0.1648732 1.381233
4 73.025 -123.225  TRUE     ARC 0.02981253 -0.07340767 0.1551282 1.444258
5 72.575 -122.775  TRUE     ARC 0.02772352 -0.06835325 0.1393356 1.590863
6 72.125 -122.325  TRUE     ARC 0.02867879 -0.05654028 0.1109921 1.988154

> results.rast
class       : SpatRaster 
dimensions  : 91, 196, 1  (nrow, ncol, nlyr)
resolution  : 0.45, 0.45  (x, y)
extent      : -141, -52.8, 42.2, 83.15  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857) 
source(s)   : memory
name        :      var 
min value   : 0.000000 
max value   : 0.177362

i have generated some reproducible code below:

library(ggplot2)
library(rgeoboundaries) #for canada shapefile
library(basemaps)
library(tidyterra)
library(terra) #for crop + mask functions

#canada shapefile
canada <- geoboundaries("Canada")

basemaps::set_defaults(map_service = "esri", map_type = "world_terrain_base")

#generate a raster
x <- raster(ncol=50, nrow=50, xmn=-141, xmx=-52, ymn=41, ymx=83)
y <- crop(x, canada)
values(y) <- 1:ncell(y)
y <- mask(y, canada)
z <- as.data.frame(y, xy = TRUE) #convert to dataframe for ggplot

ggplot() +
  #basemap_gglayer(ext = st_bbox(canada)) + 
geom_raster(subset(z, !is.na(layer)), mapping = aes(x, y, fill = layer), interpolate = TRUE) +
  coord_sf(datum = "ESRI:102001") +
  geom_sf(data = canada, fill = NA, color = "black") + 
  theme_void()

i'm not sure if this is an issue with the package, or an issue on my end. or maybe there's just a better way to get a basemap! any suggestions are appreciated

EDIT: note that the results.rast raster in my original code was created from the results dataframe, rather than the other way around (as i did in the reproducible code) - however, i assume it should work the same way. the dataframe i'm working with was created from rasters many, many steps ago (so if the issue is coming from these original rasters, i will have to re-run months of work....at which point i'd rather give up on the basemaps alltogether)

Solution

  • There's a few issues with your approach:

    The following repex addresses all of these issues by setting a consistent CRS for all data, by expanding the extent of x, and by using the correct plotting functions. Note that you did not provide copies of your colors and lg.parks objects so these have been replaced or omitted.

    Edit responding to updated question and comments A common error when working with spatial data is to incorrectly assign a CRS to the wrong coordinate values - in this case, assigning EPSG:3857 to WGS84/EPSG:4326 values. It is important to know where you are starting from in order to avoid projection issues.

    Load packages and create data:

    library(rgeoboundaries) # For Canada shapefile
library(basemaps)
library(tidyterra)
library(sf) # For sf objects, e.g. st_transform()
library(terra)
library(ggplot2)
library(ggspatial)

# Canada shapefile, project to match default basemap_raster crs
canada <- geoboundaries("Canada") %>%
  st_transform(3857)

# Basemap
set_defaults(map_service = "esri", map_type = "world_terrain_base")
bm <- basemap_terra(ext = st_bbox(canada)) %>%
  as("SpatRaster")

# Get extent of Canada sf for min/max values of example x raster
st_bbox(canada)
#      xmin      ymin      xmax      ymax 
# -15696344   5113338  -5857561  17955343 

# Generate an example raster using canada min/max values
x <- rast(ncols = 50, nrows = 50, nlyrs = 1, 
          xmin = -15696344, xmax = -5857561, 
          ymin = 5113338, ymax = 17955343, 
          names = c("rast_val"), 
          crs = "EPSG:3857") %>%
  init("cell")

y <- crop(x, canada)
values(y) <- 1:ncell(y)
y <- mask(y, canada)
z <- as.data.frame(y, xy = TRUE) # Convert to dataframe for ggplot

    Plot data:

    ggplot() +
  geom_spatraster_rgb(data = bm) +
  geom_raster(subset(z, !is.na(rast_val)),
              mapping = aes(x, y, fill = rast_val, alpha = 0.5),
              interpolate = TRUE) +
  geom_sf(data = canada, fill = NA, color = "black") +
  scale_fill_gradientn("Variance", colours = terrain.colors(5)) +
  # geom_sf(data = lg.parks, fill = NA, color = "black") + 
  annotation_scale(location = "bl", 
                   width_hint = 0.5,
                   pad_x = unit(0.18, "in")) +
  annotation_north_arrow(location = "bl",
                         which_north = "true", 
                         pad_x = unit(0.2, "in"), 
                         pad_y = unit(4, "in"),
                         style = north_arrow_fancy_orienteering) +
  guides(alpha = "none") +
  coord_sf(datum = "ESRI:102001") +
  theme_void()

    The following warning will appear:

    Scale on map varies by more than 10%, scale bar may be inaccurate

    Result:

    result