get_openstreetmap_tile_url.R

library(sf)

# load data
nc <- read_sf(system.file("shape/nc.shp", package="sf"))

# get bbox 
b <- st_bbox(nc)

# calculate tile number ---------------------------

width_or_height <- max(b["xmax"] - b["xmin"], b["ymax"] - b["ymin"])
zoom <- sum(width_or_height < 360 / 2^(0:19)) - 1 + 2
#                                                 ^^^
#                                                 Add more resolution

n <- 2^zoom

lat_deg <- b[c("xmin", "xmax")] + 180
x <- (n * lat_deg) %/% 360
lon_rad <- b[c("ymin", "ymax")] * pi / 180
sec <- function(x) 1/cos(x)
y <- (n * (1 - log(tan(lon_rad) + sec(lon_rad)) / pi)) %/% 2

xy <- unique(expand.grid(x = seq(x[1], x[2]), y = seq(y[1], y[2])))

# get tile data -----------------------

urls <- purrr::map2(xy$x, xy$y, ~glue::glue("https://a.tile.openstreetmap.org/{zoom}/{.x}/{.y}.png"))

pngs <- purrr::map(urls, httr::GET)
pngs <- purrr::map(pngs, httr::content)
pngs_g <- purrr::map(pngs, grid::rasterGrob, interpolate = TRUE)

# calculate the tile position from tile number ------------

xy_to_lonlat_ <- function(x, y, n) {
  lon_deg = x / n * 360.0 - 180.0
  lat_rad = atan(sinh(pi * (1 - 2 * y / n)))
  lat_deg = lat_rad * 180.0 / pi
  list(x = lon_deg, y = lat_deg)
}

xy_to_lonlat <- function(x, y, n) {
  nw_corner <- xy_to_lonlat_(x, y, n)
  se_corner <- xy_to_lonlat_(x + 1, y + 1, n)
  list(xmin = nw_corner$x, ymin = se_corner$y, xmax = se_corner$x, ymax = nw_corner$y)
}

pos <- purrr::map2(xy$x, xy$y, xy_to_lonlat, n = n)

# plot -------------------

ggplot(nc) +
  purrr::map2(pngs, pos, ~ annotation_raster(.x, xmin = .y$xmin, ymin = .y$ymin, xmax = .y$xmax, ymax = .y$ymax)) +
  geom_sf(fill = alpha("red", 0.3)) +
  theme_minimal()