Using ggplot to display a map plot with an underlying tiles

mappingWithGgplot.R

library(tidyverse)
library(readr)
library(ggmap)
library(ggplot2)

library(readxl)
library(plotly)
library(here)
library(httr)

# approach based entirely on with small changes
# https://yutani.rbind.io/post/2018-06-09-plot-osm-tiles/
  

deg2num<-function(lat_deg, lon_deg, zoom){
  lat_rad <- lat_deg * pi /180
  n <- 2.0 ^ zoom
  xtile <- floor((lon_deg + 180.0) / 360.0 * n)
  ytile = floor((1.0 - log(tan(lat_rad) + (1 / cos(lat_rad))) / pi) / 2.0 * n)
  return( c(xtile, ytile))
  #  return(paste(paste("https://a.tile.openstreetmap.org", zoom, xtile, ytile, sep="/"),".png",sep=""))
}

# Returns data frame containing detailed info for all zooms
deg2num.all<-function(lat_deg, lon_deg){
  nums <- as.data.frame(matrix(ncol=6,nrow=21))
  colnames(nums) <- c('zoom', 'x', 'y', 'mapquest_osm', 'mapquest_aerial', 'osm')
  rownames(nums) <- 0:20
  for (zoom in 0:20) {
    num <- deg2num(lat_deg, lon_deg, zoom)
    nums[1+zoom,'zoom'] <- zoom
    nums[1+zoom,'x'] <- num[1]
    nums[1+zoom,'y'] <- num[2]
    nums[1+zoom,'mapquest_osm'] <- paste('http://otile1.mqcdn.com/tiles/1.0.0/map/', zoom, '/', num[1], '/', num[2], '.jpg', sep='')
    nums[1+zoom,'mapquest_aerial'] <- paste('http://otile1.mqcdn.com/tiles/1.0.0/sat/', zoom, '/', num[1], '/', num[2], '.jpg', sep='')
    nums[1+zoom,'osm'] <- paste('https://a.tile.openstreetmap.org/', zoom, '/', num[1], '/', num[2], '.png', sep='')
  }
  return(nums)
}

deg2num(-36.84, 174.44, 6) # Auckland at zoom 6 - this is correct

# read some point/line/polygon data eg
nc <- sf::read_sf(system.file("C:/some/shapefile.shp"))
nc1 <- st_read("C:/someother/shapefile.shp")

# for the logic below to work we need to be in lat lon degrees using eps:4326 (not a projection to eg metres from a point like NZTM)
(nc1_WGS84 <- st_transform(nc1, 4326))

# get the bbox
(b <- sf::st_bbox(nc1_WGS84))
# this gives
# xmin      ymin      xmax      ymax 
# 167.15089 -46.90227 178.49183 -34.45623
# in comparison the original NZTM gives
# xmin    ymin    xmax    ymax 
# 1142700 4794128 2084136 6187249

class(b)
st_crs(b)

# calculate the lengths of x and y of the bbox
x_len <- b["xmax"] - b["xmin"]
y_len <- b["ymax"] - b["ymin"]

# calculate the minimum zoom level that is smaller than the lengths
x_zoom <- sum(x_len < 360 / 2^(0:19)) - 1
y_zoom <- sum(y_len < 170.1022 / 2^(0:19)) - 1

zoom <- min(x_zoom, y_zoom)
zoom

sec <- function(x) {
  1 / cos(x)
}

lonlat2xy <- function(lat_deg, lon_deg, zoom) {
  n <- 2^zoom
  
  x <- (n * (lat_deg + 180)) %/% 360
  lon_rad <- lon_deg * pi / 180
  y <- (n * (1 - log(tan(lon_rad) + sec(lon_rad)) / pi)) %/% 2
  
  list(x = x, y = y)
}

p <- ggplot(nc1) +
  geom_sf() +
  annotate("rect", xmin = b["xmin"], xmax = b["xmax"], ymin = b["ymin"], ymax = b["ymax"],
           colour = alpha("red", 0.4), fill = "transparent", linetype = "dashed", size = 1.2) +
  coord_sf(
    crs = 4326,
    expand = FALSE) 

p

corners <- expand.grid(x = b[c(1, 3)], y = b[c(2, 4)])

p +
  geom_point(aes(x, y), corners[2:3,], colour = "red", size = 5)

xy <- lonlat2xy(b[c("xmin", "xmax")], b[c("ymin", "ymax")], zoom)

tiles <- expand.grid(x = seq(xy$x["xmin"], xy$x["xmax"]),
                     y = seq(xy$y["ymin"], xy$y["ymax"]))

tiles

urls <- sprintf("https://a.tile.openstreetmap.org/%d/%d/%d.png", zoom, tiles$x, tiles$y)

xy2lonlat <- function(x, y, zoom) {
  n <- 2^zoom
  
  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(lon_deg = lon_deg, lat_deg = lat_deg)
}

library(purrr)
library(dplyr)

nw_corners <- pmap_dfr(tiles, xy2lonlat, zoom = zoom)
# add 1 to x and y to get the south-east corners
se_corners <- pmap_dfr(mutate_all(tiles, `+`, 1), xy2lonlat, zoom = zoom)

names(nw_corners) <- c("xmin", "ymax")
names(se_corners) <- c("xmax", "ymin")

tile_positions <- bind_cols(nw_corners, se_corners)
tile_positions
zoom

p +
  geom_point(aes(x, y), corners[2:3,], colour = "red", size = 5) +
  geom_rect(data = tile_positions,
            aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
            colour = "blue", fill = "transparent")

# see below for ideas on managing secrets
# https://cran.r-project.org/web/packages/httr/vignettes/secrets.html
username <- rstudioapi::showPrompt(title="username", message="Enter username", default = "")
password <- rstudioapi::askForPassword()

get_tile <- function(url) {
  # build a local path
  path <- stringr::str_extract(url, "/\\d+/\\d+/\\d+.png")
  local_png <- here::here(file.path("data", "osm-tiles", path))

  if (!file.exists(local_png)) {
    dir.create(dirname(local_png), showWarnings = FALSE, recursive = TRUE)
    
    r <- GET(
          url,
          use_proxy("someproxyserver.nz", port=8080, username=username, password=password, auth="ntlm"),
          add_headers(`User-Agent` = "From R code"),
          write_disk(local_png)
    )
  }
  png::readPNG(local_png)
}

get_tile("https://a.tile.openstreetmap.org/3/7/5.png")

pngs <- map(urls, get_tile)

args <- tile_positions %>%
  mutate(raster = pngs)

args

ggplot(nc1_WGS84) +
  pmap(args, annotation_raster, interpolate = TRUE) +
  geom_sf(fill = alpha("red", 0.3)) +
  # don't forget the license notice!
  labs(caption = "\U00a9 OpenStreetMap contributors") +
  theme_minimal()


# hmm the output doesn't quite line up - my coordinate reference systems aren't right - close but not close enough