R-script voor het downloaden en georefereren van een OpenStreetMap tile

tile.R

setwd("D:/temp")

# Functie om op basis van latitude, longitude en zoomniveau de bijbehorende OpenStreetMap tile te berekenen
# Zie http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames

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))
}

# Martinitoren:
# latitude = 53.21921
# longitude = 6.568
# zoomniveau = 18
lat <- 53.21921
lon <- 6.568
zoom <- 18

tile <- deg2num(lat, lon, zoom)
x <- tile[1]
y <- tile[2]

# Download tile
# URL template: http://[abc].tile.openstreetmap.org/zoom/x/y.png
url <- paste0(paste("http://a.tile.openstreetmap.org", zoom, x, y, sep="/"), ".png")
download.file(url, destfile = basename(url), mode = "wb")

# Genereer een World-file (*.PGW) voor het zojuist gedownloade PNG-bestand
# Bereken eerst de bounding box van de tile
tile2lon <- function(xtile, zoom) {
    return (xtile / 2.0 ^ zoom * 360.0 - 180.0);
}

tile2lat <- function(ytile, zoom) {
    n <- pi -  ytile / 2.0 ^ zoom * 2 * pi;
    return (atan(sinh(n)) * 180 / pi);
}

boundingBox <- function(xtile, ytile, zoom){
    north <- tile2lat(ytile, zoom);
    south <- tile2lat(ytile + 1, zoom);
    west <- tile2lon(xtile, zoom);
    east <- tile2lon(xtile + 1, zoom);
    return( c(north, south, east, west))
}

bb <- boundingBox(x, y, zoom)

# Maak het PGW-bestand aan
TileWidthDeg <- bb[3] - bb[4]
TileHeightDeg <- bb[1] - bb[2]
fileConn<-file(paste0(y, ".pgw"))
writeLines(c(as.character(TileWidthDeg / 256), "0.0", "0.0", as.character(-1 * TileHeightDeg / 256), as.character(bb[4]), as.character(bb[1])), fileConn)
close(fileConn)