1"mapthin" <- 2function(xy, delta, symmetric = TRUE) 3{ 4 x <- xy$x 5 y <- xy$y 6 xy <- .C(C_map_thin, 7 x = as.double(x), 8 y = as.double(y), 9 n = as.integer(length(x)), 10 as.double(delta), 11 as.integer(symmetric), 12 NAOK = TRUE)[c("x", "y", "n")] 13 length(xy$x) <- xy$n 14 length(xy$y) <- xy$n 15 xy[c("x", "y")] 16} 17 18# add axes to a map 19"map.axes" <- 20function(...) 21{ 22 axis(1, ...) 23 axis(2, ...) 24 box(...) 25 invisible() 26} 27 28"map.cities" <- 29function (x = world.cities, country = "", label = NULL, minpop = 0, 30 maxpop = Inf, capitals = 0, cex = par("cex"), projection = FALSE, 31 parameters = NULL, orientation = NULL, pch = 1, ...) 32{ 33 if (missing(x) && !exists("world.cities")) { 34 world.cities <- maps::world.cities 35 } 36 usr <- par("usr") 37 if (!missing(projection) && projection != FALSE) { 38 if (requireNamespace("mapproj", quietly = TRUE)) { 39 if (is.character(projection)) { 40 projx <- mapproj::mapproject(x$long, x$lat, projection = projection, 41 parameters = parameters, orientation = orientation) 42 } else { 43 if (nchar(mapproj::.Last.projection()$projection) > 0) { 44 projx <- mapproj::mapproject(x$long, x$lat) 45 } else stop("No projection defined\n") 46 } 47 x$long <- projx$x 48 x$lat <- projx$y 49 } else stop("mapproj package not available\n") 50 } else { 51 if (usr[2] > (180 + 0.04*(usr[2] - usr[1]))) 52 x$long[x$long < 0] <- 360 + x$long[x$long < 0] 53 } 54 selection <- x$long >= usr[1] & x$long <= usr[2] & x$lat >= usr[3] & 55 x$lat <= usr[4] & (x$pop >= minpop & x$pop <= maxpop) & ((capitals == 0) | 56 (x$capital >= 1)) 57 if (country != "") 58 selection <- selection & x$country.etc == country 59 selection0 <- selection & (x$capital == 0) & (capitals == 0) 60 selection01 <- selection & (x$capital <= 1) & (capitals <= 1) 61 selection1 <- selection & (x$capital == 1) & (capitals == 1) 62 selection2 <- selection & (x$capital == 2) & (capitals == 2) 63 selection3 <- selection & (x$capital == 3) & (capitals == 3) 64 if (is.null(label)) 65 label <- sum(selection) < 20 66 cxy <- par("cxy") 67 if (sum(selection01) > 0) 68 points(x$long[selection01], x$lat[selection01], pch = pch, 69 cex = cex * 0.6, ...) 70 if (sum(selection0) > 0) 71 if (label) 72 text(x$long[selection0], x$lat[selection0] + cxy[2] * cex * 0.7, 73 paste(" ", x$name[selection0], sep = ""), cex = cex * 0.7, ...) 74 if (sum(selection1) > 0) { 75 points(x$long[selection1], x$lat[selection1], pch = pch, cex = cex, ...) 76 if (label) { 77 text(x$long[selection1], x$lat[selection1] + cxy[2] * cex, 78 paste(" ", x$name[selection1], sep = ""), cex = cex * 1.2, ...) 79 } 80 } 81 if (sum(selection2) > 0) { 82 points(x$long[selection2], x$lat[selection2], pch = pch, cex = cex, ...) 83 if (label) { 84 text(x$long[selection2], x$lat[selection2] + cxy[2] * cex * 1.1, 85 paste(" ", x$name[selection2], sep = ""), cex = cex * 1.1, ...) 86 } 87 } 88 if (sum(selection3) > 0) { 89 points(x$long[selection3], x$lat[selection3], pch = pch, cex = cex, ...) 90 if (label) { 91 text(x$long[selection3], x$lat[selection3] + cxy[2] * cex * 0.9, 92 paste(" ", x$name[selection3], sep = ""), cex = cex * 0.9, ...) 93 } 94 } 95 invisible() 96} 97 98# draw a scale bar on a map 99"map.scale" <- 100function (x, y, relwidth = 0.15, metric = TRUE, ratio = TRUE, ...) 101{ 102 # old version 103 # format.pretty <- function(x) { 104 # as.character(pretty(x * c(0.99, 1.01), n = 2)[2]) 105 # } 106 # minka: new version 107 format.pretty <- function(x, digits = 2) { 108 x = signif(x, 2) 109 prettyNum(formatC(x, format = "fg", digits = digits), big.mark = ",") 110 } 111 usr <- par("usr") 112 if (missing(y)) 113 y <- (9 * usr[3] + usr[4])/10 114 if (abs(y) >= 90) 115 warning("location of scale out of this world!") 116 if (missing(x)) 117 #x <- (0.9 - relwidth) * usr[2] + (0.1 + relwidth) * usr[1] 118 x <- (9 * usr[1] + usr[2])/10 119 cosy <- cos((2 * pi * y)/360) 120 perdeg <- (2 * pi * (6356.78 + 21.38 * cosy) * cosy)/360 121 scale <- (perdeg * 100000)/(2.54 * (par("pin")/diff(par("usr"))[-2])[1]) 122 if (metric) 123 unit <- "km" 124 else { 125 perdeg <- perdeg * 0.6213712 126 unit <- "mi" 127 } 128 len <- perdeg * relwidth * (usr[2] - usr[1]) 129 ats <- pretty(c(0, len), n = 2) 130 nats <- length(ats) 131 labs <- as.character(ats) 132 labs[nats] <- paste(labs[nats], unit) 133 linexy <- matrix(NA, ncol = 2, nrow = 3 * nats) 134 colnames(linexy) <- c("x", "y") 135 cxy <- par("cxy") 136 dy <- cxy[2] * par("tcl") 137 dx <- ats[nats]/perdeg/(nats - 1) 138 linexy[1, ] <- c(x, y) 139 linexy[2, ] <- c(x, y + dy) 140 for (i in 1:(nats - 1)) { 141 linexy[3 * i, ] <- c(x + (i - 1) * dx, y) 142 linexy[3 * i + 1, ] <- c(x + i * dx, y) 143 linexy[3 * i + 2, ] <- c(x + i * dx, y + dy) 144 } 145 lines(linexy) 146 # minka: this is broken 147 text(x + ats/perdeg, y + dy - 0.5 * cxy[2], labs, adj = c(0.4, 0.5), ...) 148 # minka: added ratio option 149 if(ratio) 150 text(x, y + 0.5 * cxy[2], 151 paste("scale approx 1:", format.pretty(scale), sep = ""), 152 adj = 0, ...) 153 invisible(scale) 154} 155 156map.wrap <- function(p, xlim=NULL) { 157 # new version Alex Deckmyn 158 # faster, a bit more robust (check data range), bugs fixed, xlim option added 159 # insert NAs to break lines that wrap around the globe. 160 # does not work properly with polygons. 161 # p is list of x and y vectors 162 # xc is the central value of the longitude co-ordinate: usually 0, but world2 has [0,360], so 180 163 if (is.null(xlim)) { 164 xc <- 0 165 xr <- diff(range(p$x, na.rm=TRUE)) 166 } else { 167 xc <- mean(xlim) 168 xr <- diff(xlim) 169 } 170 dx <- abs(diff(p$x - xc)) 171 dax <- abs(diff(abs(p$x - xc))) 172 j <- which(dx/dax > 50 & dx > (xr * 0.8) ) 173 if (length(j)==0) return(p) 174 j <- c(j, length(p$x)) 175 ind <- seq_along(p$x) 176 177 index <- c(ind[1:j[1]], 178 unlist(lapply(1:(length(j)-1), function(k) c(NA, ind[(j[k]+1):j[(k + 1)]]) )) ) 179# notice that p$name is not returned! 180# so if wrapping is applied, you loose polygon names 181# but they would be wrong, anyway... 182 list(x = p$x[index], y = p$y[index]) 183} 184 185