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