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:
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:
This is not what I want it to do, lol. If I plot only the base layer, it looks like this:
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.
There's a few issues with your approach:
basemap_raster()
data is EPSG:3857, but the CRS for the geoboundaries()
data is EPSG:4326;geom_spatraster()
instead of geom_spatraster_rgb()
for a SpatRaster with RGB values;basemap_raster()
. This is specific to the raster
package and can cause issues if the raster
package is not loaded. For instance, if you try something like terra::project()
on a RGB SpatRaster created with basemap_raster()
without loading the raster
package, it map drop the RGB values. You should instead use basemap_terra()
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: