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