1 2osmtile_bbox = function# compute the bounding box of an OpenStreetmap tile 3### inspired by \code{osmtile} from the package \code{OpenStreetmap} 4### returns the Mercator projection bounding box 5( 6 x=61, ##<< x tile coordinate 7 y=41, ##<< x tile coordinate 8 zoom =7, ##<< zoom level 9 minim = -20037508 ##<< parameter for OSM projection 10){ 11 sc <- abs(minim) * 2 12 13 p1 <- c(x/(2^zoom) * sc + minim, -(y/(2^zoom) * sc + minim))#upper left ?? 14 p2 <- c((x + 1)/(2^zoom) * sc + minim, -((y + 1)/(2^zoom) * 15 sc + minim)) #lowerRight? 16 bbox <- list(upperLeft = p1, lowerRight = p2) 17 return(bbox) 18### bounding box, Mercator projection 19} 20 21 22plotOSM = function#plots OSM map tiles 23### places tiles on plot 24( 25 mt, ##<< list returned by \code{GetMapTiles} 26 upperLeft, ##<< upperLeft corner in lat/lon of the plot region 27 lowerRight, ##<< lowerRight corner in lat/lon of the plot region 28 lat, ##<< latitude values to be overlaid, if any 29 lon, ##<< longitude values to be overlaid, if any 30 add = FALSE, 31 removeMargin = TRUE, 32 verbose=0, ##<< level of verbosity 33 ... ##<< further arguments to be passed to \code{rasterImage} 34) { 35 map <- list(tiles = mt) 36 map$bboxLL = list() 37 38 if (missing(upperLeft) & missing(lowerRight)) { 39 x = range(mt$X) 40 y = range(mt$Y) 41 42 upperLeft = osmtile_bbox(x=x[1], y=y[1], zoom = mt$zoom)$upperLeft 43 lowerRight = osmtile_bbox(x=x[2], y=y[2], zoom = mt$zoom)$lowerRight 44 map$bbox <- list(upperLeft =upperLeft,lowerRight=lowerRight) 45 } else { 46 map$bbox <- list(upperLeft = geosphere_mercator(upperLeft[1L], upperLeft[2L]), 47 lowerRight = geosphere_mercator(lowerRight[1L], lowerRight[2L])) 48 if (verbose) print(map$bbox) 49 map$bbox <- list(p1 = c(x = min(map$bbox$upperLeft[1L], map$bbox$lowerRight[1L]), 50 y = max(map$bbox$upperLeft[2L], map$bbox$lowerRight[2L])), 51 p2 = c(x = max(map$bbox$upperLeft[1L], map$bbox$lowerRight[1L]), 52 y = min(map$bbox$p1[2L], map$bbox$lowerRight[2L]))) 53 } 54 map$bboxLL$upperLeft =geosphere_mercator(map$bbox$upperLeft,inverse = TRUE) 55 map$bboxLL$lowerRight =geosphere_mercator(map$bbox$lowerRight,inverse = TRUE) 56 57 if (verbose) print(map$bbox) 58 59 60 mar <- par("mar") 61 if (add == FALSE) { 62 plot.new() 63 if (removeMargin) 64 par(mar = c(0, 0, 0, 0)) 65 xlim_p = c(map$bbox$upperLeft[1], map$bbox$lowerRight[1]) 66 ylim_p = c(map$bbox$lowerRight[2], map$bbox$upperLeft[2]) 67 plot.window(xlim_p, ylim_p, xaxs = "i", yaxs = "i", asp = T) 68 if (verbose>1) browser() 69 } 70 #if (verbose>1) browser() 71 k=1 72 for (x in mt$X){ 73 for (y in mt$Y){ 74 osmtile=list() 75 osmtile$tile = mt$tiles[[k]] 76 osmtile$x=x;osmtile$y=y;osmtile$zoom=mt$zoom 77 plotOSMtile(osmtile, verbose=verbose, ...) 78 k=k+1 79 } 80 } 81 82 par(mar = mar) 83 invisible(map) 84### returns map object invisibly 85} 86 87plotOSMtile = function# plots a single OSM tile 88### Adds tile to plot 89( 90 osmtile, ##<< tile object 91 zoom, ##<< zoom level 92 add = TRUE, ##<< 93 raster = TRUE, ##<< 94 verbose=0, ##<< level of verbosity 95 ... ##<< further arguments to be passed to \code{rasterImage} 96){ 97 98 xres <- nrow(osmtile$tile) #look up: 256 or 255 ?!! 99 yres <- ncol(osmtile$tile) 100 if (missing(zoom)) zoom = osmtile$zoom 101 102 if (verbose>1) browser() 103 104 bbox=osmtile_bbox(x=osmtile$x, y=osmtile$y, zoom = zoom) 105 bbox$p1 = bbox$upperLeft 106 bbox$p2 = bbox$lowerRight 107 m=0#0.5 #margin 108 xleft = bbox$p1[1] - m * abs(bbox$p1[1] - bbox$p2[1])/yres 109 ybottom = bbox$p2[2] + m * abs(bbox$p1[2] - bbox$p2[2])/xres 110 xright = bbox$p2[1] - m * abs(bbox$p1[1] - bbox$p2[1])/yres 111 ytop = bbox$p1[2] + m * abs(bbox$p1[2] - bbox$p2[2])/xres 112 113 114 if (!raster) 115 image(x = seq(xleft, xright, length = yres), y = seq(ybottom, ytop, length = xres), 116 z = t(matrix(1:(xres * yres), nrow = xres, byrow = TRUE))[, 117 xres:1], col = as.vector(osmtile$tile), add = add, ...) 118 else rasterImage(osmtile$tile, xleft, ybottom, xright, ytop, ...) 119 120 if (verbose) { 121 cat("placing tile at coords:",xleft, ybottom, xright, ytop, "\n") 122 rect(xleft, ybottom, xright, ytop) 123 } 124### returns nothing 125} 126 127geosphere_mercator = function#Transform longitude/latiude points to the Mercator projection. 128### From \code{geosphere::mercator} 129( 130 p, ##<< longitude/latitude of point(s). Can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) 131 inverse = FALSE, ##<< Logical. If TRUE, do the inverse projection (from Mercator to longitude/latitude 132 r = 6378137 ##<< Numeric. Radius of the earth; default = 6378137 m 133) { 134 toRad <- pi/180 135 if (inverse) { 136 p <- .pToMatrix(p) 137 p[, 2] <- pi/2 - 2 * atan(exp(-p[, 2]/r)) 138 p[, 1] <- p[, 1]/r 139 colnames(p) <- c("lon", "lat") 140 return(p/toRad) 141 } 142 else { 143 p <- .pToMatrix(p) * toRad 144 p[, 2] <- log(tan(p[, 2]) + (1/cos(p[, 2]))) 145 p <- p * r 146 colnames(p) <- c("x", "y") 147 return(p) 148 } 149### Mercator projection of lon/lat points 150} 151 152.pToMatrix <- function(p) { 153 154 if (is.data.frame(p)) { 155 p <- as.matrix(p) 156 } else 157 158 if (is.vector(p)){ 159 if (length(p) != 2) { 160 stop('Wrong length for a vector, should be 2') 161 } else { 162 p <- matrix(p, ncol=2) 163 } 164 } else if (is.matrix(p)) { 165 if (ncol(p) != 2) { 166 stop( 'A points matrix should have 2 columns') 167 } 168 cn <- colnames(p) 169 if (length(cn) == 2) { 170 if (toupper(cn[1]) == 'Y' | toupper(cn[2]) == 'X') { 171 warning('Suspect column names (x and y reversed?)') 172 } 173 if (toupper(substr(cn[1],1,3) == 'LAT' | toupper(substr(cn[2],1,3)) == 'LON')) { 174 warning('Suspect column names (longitude and latitude reversed?)') 175 } 176 } 177 } else { 178 stop('points should be vectors of length 2, matrices with 2 columns, or inheriting from a SpatialPoints* object') 179 } 180 181 if (! is.numeric(p) ) { p[] <- as.numeric(p) 182 } 183 184 185 return(p) 186} 187 188 189