Created
May 24, 2021 02:35
-
-
Save bohdanszymanik/22ea0f185af07fdb58db545ab7f0b733 to your computer and use it in GitHub Desktop.
Using ggplot to display a map plot with an underlying tiles
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment