rmappingmapboxprojection

Why does mapboxapi return a skewed map with apparently incorrect latitude and misplaced relative to data?

I am trying to create several maps across wide latitudinal extents. I am using the Mapbox maps through mapboxapi in R. I would like to know how to do this correctly so that both the map and the data I put on it are accurately placed. Do I need to re-project?

The map apparently places longitude fairly correctly, but the Straits of Gibraltar should lie at ~36 degrees of latitude.

basemap <- layer_static_mapbox(location = c(-20,10,20,55),
                               # style_id = "streets-v12",
                               style_url = "mapbox://styles/mapbox/streets-v12",
                               username = "User",
                               access_token = "Token",
                               attribution = T,
                               logo = T)

basemap

The situation gets worse when I add the data (in "EPSG:4326"): map+data

Whereas the quick-and-dirty borders mapping accurately places the data and the world:

ggplot(data, aes(lon, lat)) +
  borders(xlim = c(-10, 10), ylim = c(10, 55)) +
  geom_hline(yintercept = 36, lty = 2, color = "cornflowerblue") +
  geom_point()

fastmap

Solution

  • Your issue arises from something to do with how the location = coordinates are translated. From the documentation regarding location:

    While the input sf object can be in an arbitrary coordinate reference system, if a length-4 bounding box vector is supplied instead it must represent WGS84 longitude/latitude coordinates and be in the order c(xmin, ymin, xmax, ymax)

    For some reason, using the c(xmin, ymin, xmax, ymax) option is not correctly being interpreted. I suspect that as the CRS for the mapbox data is WGS84 Pseudo-Mercator/EPSG:3857, something is not correctly being translated. There may be another setting required but I could not determine how to get that option to work.

    However, here is a reprex that will achieve your desired outcome. It involves using an sf object based on your coordinates that has been projected to mapbox's WGS84 Pseudo-Mercator/EPSG:3857. I have extend the coordinates slightly to reflect your example images. I have used static_mapbox() instead of layer_static_mapbox() but this approach should work regardless.

    library(mapboxapi)
library(magick)
library(terra)
library(sf)
library(ggplot2)
library(tidyterra)

# Add API key (may be prompted to restart or to run readRenviron("~/.Renviron") )
mapboxapi::mb_access_token("pk.XXXX", install = TRUE, overwrite = TRUE)
readRenviron("~/.Renviron")

# Check API key is correct
Sys.getenv("MAPBOX_PUBLIC_TOKEN")
# "pk.XXXX"

# Create sf coordinates for basemap extent (slightly larger than your example)
# and project to WGS84 Pseudo-Mercator/EPSG:3857
map_ext <- as.polygons(ext(-20, 25, 10, 55), crs = "EPSG:4326") |>
  st_as_sf() |>
  st_transform(3857)

# Create static basemap using map_ext for extent
basemap <- static_mapbox(location = map_ext,
                         buffer_dist = 0,
                         style_url = "mapbox://styles/mapbox/streets-v12",
                         username = "mapbox",
                         attribution = TRUE,
                         logo = TRUE)

# Convert basemap to raw
img_raw <- image_data(basemap, channels = "rgb")

# Convert raw to a raster array
img_array <- as.integer(img_raw)
dim(img_array) <- c(image_info(basemap)$height, image_info(basemap)$width, 3)

# Create a SpatRaster object from array
r_base <- rast(img_array)

# Assign extent and CRS
ext(r_base) <- ext(map_ext)
crs(r_base) <- "EPSG:3857"

# Example point data (if not already, you will need to project your data to "EPSG:3857")
data <- data.frame(lon = seq(-15, 10, 1),
                   lat = seq(15, 50, length.out = 26)) |>
  st_as_sf(coords = c("lon", "lat"), crs = 4326) |>
  st_transform(3857)

ggplot() +
  geom_spatraster_rgb(data = r_base) +
  geom_sf(data = data) +
  coord_sf(expand = FALSE)

    1