1map.where <- function(database = "world", x, y, ...)
2{
3  if(missing(y)) {
4    if(is.matrix(x)) { y <- x[, 2]; x <- x[, 1] }
5    else if(is.list(x) && !is.null(x$y)) { y <- x$y; x <- x$x }
6    else { y <- x[[2]]; x <- x[[1]] }
7  }
8  if(is.character(database)) {
9    mapbase <- mapenvir(database)
10    gon <- .C(C_map_where,
11       as.character(mapbase),
12       as.double(x),
13       as.double(y),
14       as.integer(length(x)),
15       integer(length(x)))[[5]]
16    # this must be database, not mapbase
17    nam <- mapname(database, ".")
18    gon[gon == 0] = NA
19    # names(nam)[gon] #AD: this only works if database.N is perfectly ordered
20    names(nam)[match(gon,nam)]
21  }
22  else {
23    if (inherits(database,"SpatialPolygons")) {
24      database <- SpatialPolygons2map(database, ...)
25    }
26    if(num.polygons(database) != length(database$names))
27      stop("map object must have polygons (fill=TRUE)")
28    n = length(database$x)
29    result = .C(C_map_in_polygon,
30       as.double(database$x), as.double(database$y), as.integer(n),
31       as.double(x), as.double(y), as.integer(length(x)),
32       integer(length(x)), NAOK = TRUE)[[7]]
33    result[result == 0] = NA
34    database$names[result]
35  }
36}
37as.matrix.polygon <- function(x, ...) {
38  p = x
39  if(is.null(p)) return(p)
40  if(is.list(p) && !is.data.frame(p)) p <- cbind(p$x, p$y)
41  p
42}
43in.one.polygon <- function(p, x) {
44  # returns a logical vector, whose length is nrow(x)
45  if(is.null(p)) return(NA)
46  p <- as.matrix.polygon(p)
47  if(is.list(x) && !is.data.frame(x)) x <- cbind(x$x, x$y)
48  if(is.vector(x)) dim(x) <- c(1, 2)
49  # p and x are matrices
50  .C(C_map_in_one_polygon,
51     as.double(p[, 1]), as.double(p[, 2]), as.integer(nrow(p)),
52     as.double(x[, 1]), as.double(x[, 2]), as.integer(nrow(x)),
53     logical(nrow(x)), as.integer(TRUE))[[7]]
54}
55in.polygon <- function(p, x) {
56  # returns a logical vector, whose length is nrow(x)
57  if(is.null(p)) return(NA)
58  p <- as.matrix.polygon(p)
59  if(is.list(x) && !is.data.frame(x)) x <- cbind(x$x, x$y)
60  if(is.vector(x)) dim(x) <- c(1, 2)
61  # p and x are matrices
62  .C(C_map_in_polygon,
63     as.double(p[, 1]), as.double(p[, 2]), as.integer(nrow(p)),
64     as.double(x[, 1]), as.double(x[, 2]), as.integer(nrow(x)),
65     logical(nrow(x)), NAOK = TRUE)[[7]] > 0
66}
67
68# polygon is not assumed closed
69area.polygon <- function(p) {
70  if(is.null(p)) return(NA)
71  p <- as.matrix.polygon(p)
72  n <- nrow(p)
73  x1 <- p[, 1]
74  i2 <- c(n, 1:(n - 1))
75  x2 <- p[i2, 1]
76  y1 <- p[, 2]
77  y2 <- p[i2, 2]
78  0.5*abs(sum(x1*y2 - x2*y1))
79}
80
81centroid.polygon <- function(p) {
82  if(is.null(p)) return(c(NA, NA))
83  p <- as.matrix.polygon(p)
84  n <- nrow(p)
85  x1 <- p[, 1]
86  i2 <- c(n, 1:(n - 1))
87  x2 <- p[i2, 1]
88  y1 <- p[, 2]
89  y2 <- p[i2, 2]
90  a <- x1*y2 - x2*y1
91  s <- sum(a)*3
92  if(s == 0) c(mean(x1), mean(y1))
93  else c(sum((x1 + x2)*a)/s, sum((y1 + y2)*a)/s)
94}
95
96# applies fun to all sub-polygons of p
97apply.polygon <- function(p, fun, names. = NULL) {
98  if(is.null(p)) return(p)
99  if(is.null(names.) && !is.null(p$names)) names. = p$names
100  p = as.matrix.polygon(p)
101  n <- nrow(p)
102  breaks <- (1:n)[is.na(p[, 1])]
103  starts <- c(1, breaks + 1)
104  ends <- c(breaks - 1, n)
105  m <- length(starts)
106  result <- list()
107  for(i in 1:m) {
108    this.p = if(ends[i] >= starts[i]) p[starts[i]:ends[i], ] else NULL
109    result[[i]] <- fun(this.p)
110  }
111  names(result) = names.
112  result
113}
114
115num.polygons <- function(p) {
116  if(is.list(p)) 1 + sum(is.na(p$x))
117  else 1 + sum(is.na(p[, 1]))
118}
119
120# range.polygon <- function(..., na.rm = FALSE) {
121#   p <- as.list.polygon(...)
122#   lapply(p[c("x", "y")], range, na.rm = na.rm)
123# }
124
125map.text <- function(database, regions = ".", exact = FALSE, labels,
126    cex = 0.75, add = FALSE, move = FALSE, ...) {
127  if(!add) map(database=database, regions=regions, exact=exact, ...)
128  # get polygons
129  cc = match.call(expand.dots=TRUE)
130  cc[[1]] = as.name("map")
131  cc$fill = TRUE
132  cc$plot = FALSE
133  cc$regions = regions
134  cc$exact = exact
135  cc$move = cc$add = cc$cex = cc$labels = NULL
136  cc$resolution = 0
137  m = eval(cc)
138  if(missing(labels)) {
139    labels = gsub(".*,", "", m$names)
140  }
141  if(num.polygons(m) != length(labels))
142    stop("map object must have polygons (fill=TRUE) and equal number of labels")
143  x = apply.polygon(m, centroid.polygon)
144  # convert m into a matrix
145  x <- t(array(unlist(x), c(2, length(x))))
146  if(move) {
147# AD: this option should probably be removed, as the code is not available
148    # require(mining)
149    move.collisions2 <- get("move.collisions2")	# to prevent check NOTE
150    w = strwidth(labels, units = "inches", cex = cex)
151    h = strheight(labels, units = "inches", cex = cex)
152    x = move.collisions2(x[, 1], x[, 2], w, h)
153  }
154  # want to omit map-specific options here (like "exact")
155  text(x, labels, cex = cex, ...)
156  invisible(m)
157}
158
159identify.map <- function(x, n = 1, index = FALSE, ...) {
160  # identify polygons in a map
161  # must click near the center of the polygon
162  m = x
163  if(!is.list(m)) stop("must provide a map object")
164  if(num.polygons(m) != length(m$names))
165    stop("map object must have polygons (fill=TRUE)")
166  x = apply.polygon(m, centroid.polygon)
167  x <- t(array(unlist(x), c(2, length(x))))
168  i = identify(x[, 1], x[, 2], labels = m$names, n = n, ...)
169  if(index) i else m$names[i]
170}
171
172area.map <- function(m, regions = ".", sqmi=TRUE, ...) {
173  # returns the areas of given regions,
174  # combining the areas of all regions which match.
175  if(!is.list(m)) stop("must provide a map object")
176  if(num.polygons(m) != length(m$names))
177    stop("map object must have polygons (fill=TRUE)")
178  proj = m$projection
179  m = map.poly(m,regions,as.polygon=TRUE,...)
180  area = unlist(apply.polygon(m, area.polygon))
181  merge <- regions[match.map(m, regions, ...)]
182  names(merge) <- m$names
183  merge = factor(merge, levels = regions)
184  area = drop(indicators.factor(merge) %*% area)
185  areaSqMiles <- function(proj) {
186    # returns a factor f such that f*area.map() is in square miles.
187    if(is.null(proj)) proj = "no projection"
188    if(proj %in% c("mollweide","azequalarea","aitoff")) 2*15732635
189    else if(proj %in% c("sinusoidal","bonne","cylequalarea","albers")) 15732635
190    else if(proj == "sp_albers") 15745196
191    else {
192      warning(paste("sq.mile correction unavailable for",proj))
193      1
194    }
195  }
196  if(sqmi) area*areaSqMiles(proj) else area
197}
198indicators.factor <- function(y) {
199  # convert a factor into a matrix of indicators
200  # result is level by case
201  # works if y contains NAs
202  r <- array(0, c(length(levels(y)), length(y)), list(levels(y), names(y)))
203  for(lev in levels(y)) r[lev, y == lev] <- 1
204  r
205}
206
207