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