mercator-web-map-tile-conversion.R

# Conversion between Web Map Tiles and Mercator coordinates

require(tibble)

merc_to_wmt = function(x, y, z){
  m_ext = pi * 6378137 # Mercator extent
  res = 2 ^ (1:14)[z]  # grid size
  x_grid = scales::rescale(x, from = c(-m_ext, m_ext), to = c(0, res)) |> floor()
  y_grid = scales::rescale(y, from = c(m_ext, -m_ext), to = c(0, res)) |> floor()
  tibble(x = x_grid |> pmax(0) |> pmin(res-1), y = y_grid |> pmax(0) |> pmin(res-1))
}

wmt_to_merc = function(x, y, z){ 
  m_ext = pi * 6378137 # Mercator extent
  res = 2 ^ (1:14)[z]  # grid size
  x_grid = scales::rescale(x, from = c(0, res), to = c(-m_ext, m_ext))
  y_grid = scales::rescale(y, from = c(0, res), to = c(m_ext, -m_ext))
  tibble(x = x_grid, y = y_grid)
}

# Usage and validation

require(dplyr)

# test data in WGS84
x = runif(n = 1e4, min = -180, max = 180)
y = runif(n = 1e4, min =  -90, max =  90)
z = 5

# project to Mercator
m = tibble(x = x, y = y) %>% 
  sf::sf_project(from = leaflet:::epsg4326, to = leaflet:::epsg3857, pts = .) |> 
  as_tibble() |> `names<-`(c('x', 'y'))

# this function buried in {rosm} converts WGS84 to WMT so we'll use to benchmark
rosm:::tile.xy(180, -90, 5)
rosm:::tile.xy(180, 90, 6)
rosm:::tile.xy(0, 0, 6)

# rosm conversion (very slow)
w = purrr::map2(x, y, rosm:::tile.xy, zoom = z) |> 
  purrr::map_df(~ t(.x) |> as_tibble()) |> `names<-`(c('wx', 'wy'))

# our convert to WMT 
w2 = merc_to_wmt(m$x, m$y, z)

# collate
d = tibble(x = x, y = y) |> cbind(dplyr::rename(m, mx=x, my=y)) |> 
  cbind(w) |> cbind(`names<-`(w2, c('wx2','wy2'))) |> as_tibble()

d |> dplyr::filter(wy2 != wy | wx2 != wx)  # identical results

# WMT to Mercator
wmt_to_merc(d$wx, d$wy, z)