rggplot2mapsrastertidyterra

Basemap in basemaps_gglayer() not working

I'm trying to add a basemap to some of my figures. It seems like the basemaps package is the best way to go, but I haven't been able to get it to work.

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

When I try adding a basemap with 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 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"

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 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") as 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.

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.

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