Overlay sf object on ggmap - not aligned, specific to region

Question

I would like to overlay an sf object (mapping geometry with organization's data) over the google map showing terrain features. This was successful when using other US states, but not Alaska. What do I need to do/change in order to overlay the map more precisely?

ak1<-ggmap(AK.base)+
  geom_sf(data=ak_region_download,mapping=aes(geometry=geometry,fill=n,color=''),inherit.aes = F,alpha=.5)+
  scale_fill_continuous(na.value=NA,low = "darkblue", high = "orange", 
                        name = 'number',
                        breaks=c(1,5,10,20,30,40))+
  scale_color_discrete('')+
  theme_bw()
ak1

For reference I achieved AK.base with a google API

AK.base <- get_map(location = c(lon=-152,lat=63.5), zoom = 4, scale = "auto", 
                   maptype = c("terrain"), 
                   messaging = FALSE, urlonly = FALSE, filename = "ggmapTemp", 
                   crop = FALSE, color = c("color"), source = c("google"))

And the dataframe looks something like this, using the geojson file from the website referenced below.

# https://gis.data.alaska.gov/datasets/DCCED::alaska-borough-and-census-area-boundaries/about
#
# Download shapefile from above link, place ZIP in folder called DataRaw within 
# project directory, unzip
ak_region_download <-
  here("DataRaw/Alaska_Borough_and_Census_Area_Boundaries.shp") |> 
  sf::st_read()
ak_region_download$n <- 1:30

Answer 1

There two main issues you're encountering:

• the get_map() data from Google are in EPSG:3857, but get loaded into R as WGS84/EPSG:4326 without reprojecting, hence the offset
• your AK.base object crosses the anti-meridian, which makes solutions like this one trickier to implement

The fix is a bit convoluted, but I could not work out a more succinct approach.

First, let's load the required packages and examine the extent of the get_map() data:

library(sf)
library(terra)
library(ggmap)
library(tidyterra)

AK.base <- get_map(
  location = c(lon = -152, lat = 63.5),
  zoom = 4,
  scale = "auto",
  maptype = c("terrain"),
  messaging = FALSE,
  urlonly = FALSE,
  filename = "ggmapTemp",
  crop = FALSE,
  color = c("color"),
  source = c("google"))

# View extent of get_map() object
ak_bbox <- setNames(unlist(attr(AK.base, "bb")), 
                    c("ymin", "xmin", "ymax", "xmax"))

ak_bbox[c("xmin", "xmax", "ymin", "ymax")]
#       xmin       xmax       ymin       ymax 
# -180.08105 -123.83105   47.88752   73.58465

Note that the values are WGS84/EPSG:4326 coordinates and the the xmin value is outside the bounding extent of WGS84/EPSG:4326 (-180, 180, -90, 90). To ensure the extent is valid, you can first convert AK.base to a SpatRaster object and crop it using terra::crop():

# Convert get_map() object to SpatRaster
r <- rast(AK.base)

# Create extent object that does not cross anti-meridian
ext_crop <- ext(-180, ext(r)[2:4])

# Crop r
r <- crop(r, ext_crop)

# Get new extent of SpatRaster
r_ext <- ext(r)
setNames(c(r_ext[1:4]), 
         c("xmin", "xmax", "ymin", "ymax"))
#       xmin       xmax       ymin       ymax 
# -179.99316 -123.83105   47.88752   73.58465

The xmin is now within the permitted bounding extent of WGS84/EPSG:4326. However, these values should be EPSG:3857, which is Google's default CRS. To correct this:

# Convert r extent coordinates to EPSG:3857
ext3857 <- ext(project(vect(ext(r), crs = crs(r)), "EPSG:3857"))

# Change (and not project) CRS of r
crs(r) <- "EPSG:3857"

# Use ext3857 coordinates to modify extent values of r
ext(r) <- ext3857

r
# class       : SpatRaster 
# dimensions  : 1280, 1278, 3  (nrow, ncol, nlyr)
# resolution  : 4891.97, 4891.97  (x, y)
# extent      : -20036747, -13784809, 6088163, 12349884  (xmin, xmax, ymin, ymax)
# coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857) 
# source(s)   : memory
# colors RGB  : 1, 2, 3 
# names       : red, green, blue 
# min values  :  12,     8,    8 
# max values  : 254,   254,  254

As shown, the data now have the correct CRS and extent.

Finally, plot your data. As r is a SpatRaster, you can't use ggmap() so use tidyterra::geom_spatraster_rgb() instead:

# Project ak_region_download to same CRS as SpatRaster
ak_region_download <- st_transform(ak_region_download, st_crs(r))

# Plot
ggplot() +
  geom_spatraster_rgb(data = r, maxcell = Inf) +
  geom_sf(data = ak_region_download,
          aes(fill = n, color = NA),
          inherit.aes = F,
          alpha = .5) +
  scale_fill_continuous(na.value = NA, low = "darkblue", high = "orange",
                        name = "number",
                        breaks = c(1,5,10,20,30,40)) +
  coord_sf(xlim = ext(r)[c(1,2)],
           expand = FALSE) +
  scale_color_discrete("")+
  theme_bw()

If you examine the result, you will see that the Canada/US border does not align with your sf data. I checked the Google data against:

ak_sf <- tigris::states() |>
  filter(NAME == "Alaska") |>
  st_transform(st_crs(r))

and the same issue occurred. I can't say definitively whether this a projection issue or an issue with the data, but if you want to avoid the misaligned borders, you can use ggmap::get_googlemap(). It allows you to modify the default output using the style = parameter. The example below sets the country borders to "off":

AK.base <- get_googlemap(center = c(lon = -152, lat = 63.5),
                         source = "google",
                         maptype = "terrain",
                         scale = 2,
                         zoom = 4,
                         style = c(feature = "administrative.country",
                                   element = "geometry.stroke",
                                   visibility = "off")
                         )

Applying the rest of the code results in: