.packageName <- "maps"
"mapgetg" <-
function(database = "world", gons, fill = FALSE, xlim = c(-1e30, 1e30),
	ylim = c(-1e30, 1e30))
{
	ngon <- length(gons)
	gnames <- names(gons)
	dbname <- paste(database, "MapEnv", sep = "")
	data(list = dbname)
	mapbase <- paste(Sys.getenv(get(dbname)), database, sep = "")
	z <- .C("mapgetg", PACKAGE="maps",
		as.character(mapbase),
		gons = as.integer(gons),
		as.integer(ngon),
		sizes = integer(ngon),
		error = as.integer(0),
		as.double(c(xlim, ylim)),
		as.integer(fill))[c("gons", "sizes", "error")]
	gons <- z$gons
	sizes <- z$sizes
	if(z$error < 0)
		stop("error in reading polygon headers")
	z <- .C("mapgetg", PACKAGE="maps",
		as.character(mapbase),
		as.integer(gons),
		as.integer(ngon),
		lines = integer(sum(sizes)),
		error = as.integer(1),
		as.double(c(xlim, ylim)),
		as.integer(fill))[c("lines", "error")]
	if(z$error < 0)
		stop("error in reading polyline numbers")
	lines <- z$lines
	ok <- sizes > 0
	list(number = lines, size = sizes[ok], name = gnames[ok])
}

"mapgetl" <-
function(database = "world", lines, xlim = c(-1e30, 1e30), ylim = c(-1e30,
	1e30), fill = FALSE)
{
	nline <- as.integer(length(lines))
	if(nline == 0)
		return(integer(0))
	dbname <- paste(database, "MapEnv", sep = "")
	data(list = dbname)
	mapbase <- paste(Sys.getenv(get(dbname)), database, sep = "")
	z <- .C("mapgetl", PACKAGE="maps",
		as.character(mapbase),
		linesize = as.integer(lines),
		error = as.integer(nline),
		as.integer(0),
		as.double(0),
		as.double(0),
		as.double(c(xlim, ylim)),
		as.integer(fill))[c("linesize", "error")]
	if(z$error < 0)
		return(integer(0))
	ok <- z$linesize != 0
	lines <- lines[ok]
	nline <- length(lines)
	if(nline == 0)
		return(integer(0))
	linesize <- z$linesize[ok]
	N <- sum(linesize) + nline - 1
	xy <- .C("mapgetl", PACKAGE="maps",
		as.character(mapbase),
		as.integer(lines),
		as.integer(nline),
		as.integer(1),
		x = double(N),
		y = double(N),
		range = double(4),
		as.integer(fill))[c("x", "y", "range")]
	flip <- c(1, -1, 1, -1)
	box <- flip * pmax(flip * c(xlim, ylim), flip * xy$range)
	if(any(diff(box)[-2] <= 0)) {
		if(missing(xlim) || missing(ylim))
			stop("nothing to draw: data range and limits don't intersect"
				)
		else box <- c(xlim, ylim)
	}
	xy$range <- box
	xy
}

"mapname" <-
function(database = "world", patterns, exact = FALSE)
{
  dbname <- paste(database, "MapEnv", sep = "")
  data(list = dbname)
  mapbase <- paste(Sys.getenv(get(dbname)), database, sep = "")
  # rewritten by Tom Minka
  fname <- paste(sep = "", mapbase, ".N")
  cnames <- read.delim(fname, as.is = TRUE, header = FALSE)
  nam <- as.character(cnames[[1]])
  if(exact) {
    i = match(patterns, nam)
    if(any(is.na(i))) i = NULL
  } else {
    regexp <- paste("(^", patterns, ")", sep = "", collapse = "|")
    i <- grep(regexp, nam, ignore.case = TRUE)
  }
  if(length(i) == 0) return(NULL)
  r <- cnames[i, 2]
  names(r) <- nam[i]
  return(r)
}

"maptype" <-
function(database = "world")
{
  if(is.character(database)) {
	dbname <- paste(database, "MapEnv", sep = "")
	data(list = dbname)
	mapbase <- paste(Sys.getenv(get(dbname)), database, sep = "")
        # minka: maptypes are now 1,2 instead of 0,1
	switch(.C("maptype", PACKAGE="maps",
		as.character(mapbase),
		integer(1))[[2]] + 2, "unknown", "spherical", "planar", "spherical")
  } else {
    # map object
    "spherical"
  }
}

char.to.ascii <- function(s) {
  # returns the ascii code for a character (0 for an empty string)
  n = length(s)
  .C("char_to_ascii", PACKAGE="maps",
    as.integer(n), as.character(s), integer(n))[[2]]
}
is.regexp <- function(s) {
  # test: is.regexp(c("too", ".*dak", "kj[t]"))
  pattern = ".*[\.\?\*\[].*"
  result = logical(length(s))
  result[grep(pattern, s)] = TRUE
  result
}
match.map <- function(database, regions, exact = FALSE, warn = TRUE) {
  invperm <- function(p) {
    # returns the inverse permutation, so that invperm(p)[p] = 1:n
    ip = p
    ip[p] = 1:length(p)
    ip
  }
  if(is.character(database)) {
    dbname <- paste(database, "MapEnv", sep = "")
    data(list = dbname)
    mapbase <- paste(Sys.getenv(get(dbname)), database, sep = "")
    fname <- paste(sep = "", mapbase, ".N")
    x <- read.delim(fname, header = FALSE)
    nam <- as.character(x[[1]])
  }
  else {
    nam <- database$names
  }
  nam = tolower(nam)
  regions = tolower(regions)
  if(!exact && any(is.regexp(regions))) {
    match.map.grep(nam, regions, warn)
  } else {
    # sort regions and table
    ord.nam = order(nam)
    nam = nam[ord.nam]
    ord.regions = order(regions)
    regions = regions[ord.regions]
    result = .C("map_match", PACKAGE="maps",
      as.integer(length(nam)), as.character(nam),
      as.integer(length(regions)), as.character(regions),
      result = integer(length(nam)), as.integer(exact))[["result"]]
    # 0 -> NA
    is.na(result[result == 0]) = TRUE
    # back to original order
    ord.regions[result[invperm(ord.nam)]]
    #result
  }
}
match.map.slow <- function(nam, regions, warn = FALSE) {
  result <- rep(NA, length = length(nam))
  # use hashing to prune the potential prefix matches
  # doesn't work for regexp
  let = min(nchar(regions))
  names(nam) = NULL
  hash = factor(char.to.ascii(substr(nam, let, let)))
  nam.bin = split(nam, hash)
  index.bin = split(1:length(nam), hash)
  for(i in 1:length(regions)) {
    pattern = regions[i]
    regexp <- paste("(^", pattern, ")", sep = "", collapse = "|")
    #r <- grep(regexp, nam)
    region.hash = as.character(char.to.ascii(substr(pattern, let, let)))
    r <- grep(regexp, nam.bin[[region.hash]])
    if(length(r) > 0) {
      r = index.bin[[region.hash]][r]
      result[r] <- i
    } else if(warn) warning(paste(pattern, "is not in the map"))
  }  
  result
}
match.map.grep <- function(nam, regions, warn = FALSE) {
  result <- rep(NA, length = length(nam))
  for(i in 1:length(regions)) {
    pattern = regions[i]
    regexp <- paste("(^", pattern, ")", sep = "", collapse = "|")
    r <- grep(regexp, nam)
    if(length(r) > 0) {
      result[r] <- i
    } else if(warn) warning(paste(pattern, "is not in the map"))
  }  
  result
}
# R data structure for maps is list(x, y, names)
# names is character vector naming the polygons
# (contiguous groups of non-NA coordinates)

# returns the jth contiguous section of non-NA values in vector x
# should consecutive NA's count as one?
subgroup <- function(x, i) {
  n <- length(x)
  breaks <- (1:n)[is.na(x)]
  if (length(breaks) == 0) {
    starts <- 1; ends <- n
  } else {
    starts <- c(1, breaks + 1)
    ends <- c(breaks - 1, n)
  }
  result <- numeric(0)
  for(j in i) {
    p <- x[starts[j]:ends[j]]
    if (length(result) == 0) result <- p
    else result <- c(result, NA, p)
  }
  result
}

sub.polygon <- function(p, i) {
  lapply(p[c("x", "y")], function(x) subgroup(x, i))
}

# returns a sub-map of the named map corresponding to the given regions
# regions is a vector of regular expressions to match to the names in the map
# regions outside of xlim, ylim may be omitted
map.poly <- function(database, regions = ".", exact = FALSE,
                     xlim = NULL, ylim = NULL, boundary = TRUE,
		     interior = TRUE, fill = FALSE, as.polygon = FALSE) {
  if (!is.character(database)) {
    if (!as.polygon) stop("map objects require as.polygon=TRUE")
    the.map <- database
    if (identical(regions,".")) {
      # speed up the common case
      nam = the.map$names
      coord <- the.map[c("x", "y")]
    } else {
      # same as mapname()
      if (exact) {
        i = match(regions, the.map$names)
        if (any(is.na(i))) i = NULL
      } else {
      regexp <- paste("(^", regions, ")", sep = "", collapse = "|")
        i <- grep(regexp, the.map$names, ignore.case = TRUE)
      }
      if (length(i) == 0) stop("no recognized region names")
      nam <- the.map$names[i]
      coord <- sub.polygon(the.map, i)
    }
    coord$range <- c(range(coord$x, na.rm = TRUE), range(coord$y, na.rm = TRUE))
  } else {
    # turn the polygon numbers into a list of polyline numbers
    gon <- mapname(database, regions, exact)
    n <- length(gon)
    if (n == 0) stop("no recognized region names")
    if (is.null(xlim)) xlim <- c(-1e+30, 1e+30)
    if (is.null(ylim)) ylim <- c(-1e+30, 1e+30)
    # turn the polygon numbers into a list of polyline numbers
    line <- mapgetg(database, gon, as.polygon, xlim, ylim)
    if (length(line$number) == 0)
            stop("nothing to draw: all regions out of bounds")
    # turn the polyline numbers into x and y coordinates
    if (as.polygon) {
      coord <- mapgetl(database, line$number, xlim, ylim, fill) 
      # assemble lines into polygons
      gonsize <- line$size
      keep <- rep(TRUE, length(gonsize))
      coord[c("x", "y")] <- makepoly(coord, gonsize, keep)
    }
    else {
      l <- abs(line$number)
      if (boundary && interior) l <- unique(l)
      else if (boundary) l <- l[!match(l, l[duplicated(l)], FALSE)]
      else l <- l[duplicated(l)]
      coord <- mapgetl(database, l, xlim, ylim, fill)
      if (length(coord) == 0)
              stop("all data out of bounds")
    }
    nam <- line$name
  }
  list(x = coord$x, y = coord$y, range = coord$range, names = nam)
}

map <-
function(database = "world", regions = ".", exact = FALSE,
	 boundary = TRUE, interior = TRUE, projection = "",
	 parameters = NULL, orientation = NULL, fill = FALSE,
	 col = 1, plot = TRUE, add = FALSE, namesonly = FALSE, 
         xlim = NULL, ylim = NULL, wrap = FALSE,
         resolution = if (plot) 1 else 0, type = "l", bg = par("bg"),
         mar = c(0, 0, par("mar")[3], 0.1), border = 0.01, ...)
{
  # parameter checks
  if (resolution>0 && !plot) 
    stop("must have plot=TRUE if resolution is given")
  if (!fill && !boundary && !interior)
    stop("one of boundary and interior must be TRUE")
  doproj <- !missing(projection) || !missing(parameters) || !missing(
          orientation)
  coordtype <- maptype(database)
  if (coordtype == "unknown") 
     stop("missing database or unknown coordinate type")
  if (doproj && coordtype != "spherical") 
    stop(paste(database, "database is not spherical; projections not allowed"))
  # turn the region names into x and y coordinates
  if (is.character(database)) as.polygon = fill
  else as.polygon = TRUE
  coord <- map.poly(database, regions, exact, xlim, ylim, 
                    boundary, interior, fill, as.polygon)
  if (is.na(coord$x[1])) stop("first coordinate is NA.  bad map data?")
  if (plot) {
    assign(".map.range", coord$range, envir = globalenv())
  }
  if (doproj) {
    nam <- coord$names
    library(mapproj)
    coord <- mapproject(coord, pr = projection, pa = parameters,
                        or = orientation)
    coord$projection = projection
    coord$parameters = parameters
    coord$orientation = orientation
    if (plot && coord$error)
      if (all(is.na(coord$x)))
        stop("projection failed for all data")
      else warning("projection failed for some data")
    coord$names <- nam
  }
  # do the plotting, if requested
  if (plot) {
    # for new plots, set up the coordinate system;
    # if a projection was done, set the aspect ratio
    # to 1, else set it so that a long-lat square appears
    # square in the middle of the plot
    if (!add) {
      opar = par(bg = bg)
      if (!par("new")) plot.new()
      # xlim, ylim apply before projection
      if (is.null(xlim) || doproj) xrange <- range(coord$x, na.rm = TRUE)
      else xrange <- xlim
      if (is.null(ylim) || doproj) yrange <- range(coord$y, na.rm = TRUE)
      else yrange <- ylim
      if (coordtype != "spherical" || doproj) {
	aspect <- c(1, 1) 
      } else
        aspect <- c(cos((mean(yrange) * pi)/180), 1)
      d <- c(diff(xrange), diff(yrange)) * (1 + 2 * border) * aspect
      if (coordtype != "spherical" || doproj) {
        plot.window(xrange, yrange, asp = 1/aspect[1])
      } else {
        # must have par(xpd = FALSE) for limits to have an effect ??!
	p <- par("fin") -
	  as.vector(matrix(c(0, 1, 1, 0, 0, 1, 1, 0), nrow = 2) %*% par("mai"))
	par(pin = p)
        p <- par("pin")
        p <- d * min(p/d)
        par(pin = p)
        d <- d * border + ((p/min(p/d) - d)/2)/aspect
        usr <- c(xrange, yrange) + rep(c(-1, 1), 2) * rep(d, c(2, 2))
        par(usr = usr)
      }
      on.exit(par(opar))
    }
    # thinning only works if you have polylines from a database
    if (is.character(database) && resolution != 0 && type != "n") {
      pin <- par("pin")
      usr <- par("usr")
      resolution <- resolution * min(diff(usr)[-2]/pin/100)
      coord[c("x", "y")] <- mapthin(coord, resolution)
    }
    if (type != "n") {
      if (wrap) coord = map.wrap(coord)
      if (fill) polygon(coord, col = col, ...)
      else lines(coord, col = col, type = type, ...)
    }
  }
  # return value is names or coords, but not both
  class(coord) = "map"
  value <- if (namesonly) coord$names else coord
  if (plot) invisible(value)
  else value
}

"makepoly" <-
function(xy, gonsize, keep)
{
  # remove NAs and duplicate points so that a set of polylines becomes a set
  # of polygons.
  # xy is a set of polylines, separated by NAs.
  # gonsize is a vector, giving the number of lines to put into each polygon
  # note that a polyline may consist of a single point
  x <- xy$x
  y <- xy$y
  n <- length(x)
  gonsize <- gonsize[ - length(gonsize)]
  discard <- seq(length(x))[is.na(x)]
  if (length(discard) > 0) {
    # locations of (possible) duplicate points
    dups = c(discard - 1, n)
    # only polylines with > 1 point have duplicates
    i = which(diff(c(0, dups)) > 2);
    discard <- c(discard, dups[i]);
  }
  # first part of discard is the NAs, second part is duplicates
  # gonsize tells us which NAs to preserve
  if (length(gonsize) > 0)
    discard <- discard[ - cumsum(gonsize)]
  if (length(discard) > 0) {
    x <- x[ - discard]
    y <- y[ - discard]
  }
  keep <- rep(keep, diff(c(0, seq(length(x))[is.na(x)], length(x))))
  closed.polygon(list(x = x[keep], y = y[keep]))
}
closed.polygon <- function(p) {
  # p is a set of polylines, separated by NAs
  # for each one, the first point is copied to the end, giving a closed polygon
  x = p$x
  y = p$y
  n = length(x)
  breaks <- seq(length(x))[is.na(x)]
  starts <- c(1, breaks + 1)
  ends <- c(breaks - 1, n)
  x[ends + 1] = x[starts]
  y[ends + 1] = y[starts]
  x = insert(x, breaks + 1)
  y = insert(y, breaks + 1)
  list(x = x, y = y)
}
insert <- function(x, i, v = NA) {
  # insert v into an array x, at positions i
  # e.g. insert(1:7, c(2, 5, 8))
  n = length(x)
  new.n = n + length(i)
  m = logical(new.n)
  i.new = i - 1 + seq(length(i))
  m[i.new] = TRUE
  x = x[(1:new.n) - cumsum(m)]
  x[i.new] = v
  x
}
"mapthin" <-
function(xy, delta, symmetric = TRUE)
{
	x <- xy$x
	y <- xy$y
	xy <- .C("mapthin", PACKAGE="maps",
		x = as.double(x),
		y = as.double(y),
		n = as.integer(length(x)),
		as.double(delta),
		as.integer(symmetric),
		NAOK = TRUE)[c("x", "y", "n")]
	length(xy$x) <- xy$n
	length(xy$y) <- xy$n
	xy[c("x", "y")]
}

# add axes to a map
"map.axes" <-
function()
{
	axis(1)
	axis(2)
	box()
	invisible()
}

"map.cities" <-
function(x = world.cities, country = "", label = NULL, minpop = 0, maxpop = Inf, 
	capitals = 0, cex = par("cex"), ...)
{
        if(missing(x)) data(world.cities)	# Using lazy evaluation
	usr <- par("usr")
	if(usr[2] > 180)
		x$long[x$long < 0] <- 360 + x$long[x$long < 0]
	selection <- x$long >= usr[1] & x$long <= usr[2] & x$lat >= usr[3] & x$
		lat <= usr[4] & (x$pop >= minpop & x$pop <= maxpop) & ((
		capitals == 0) | (x$capital >= 1))
	if(country != "")
		selection <- selection & x$country.etc == country
	selection0 <- selection & (x$capital == 0) & (capitals == 0)
	selection01 <- selection & (x$capital <= 1) & (capitals <= 1)
	selection1 <- selection & (x$capital == 1) & (capitals == 1)
	selection2 <- selection & (x$capital == 2) & (capitals == 2)
	selection3 <- selection & (x$capital == 3) & (capitals == 3)
	if(is.null(label))
		label <- sum(selection) < 20
	cxy <- par("cxy")
	if(sum(selection01) > 0)
		points(x$long[selection01], x$lat[selection01], pch = 1, cex = 
			cex * 0.6, ...)
	if(sum(selection0) > 0)
		if(label)
			text(x$long[selection0], x$lat[selection0] + cxy[
				2] * cex * 0.7, paste(" ", x$name[selection0], 
				sep = ""), cex = cex * 0.7, ...)
	if(sum(selection1) > 0) {
		points(x$long[selection1], x$lat[selection1], pch = 1, cex = 
			cex, ...)
		text(x$long[selection1], x$lat[selection1] + cxy[2] * cex,
			paste(" ", x$name[selection1], sep = ""), cex = cex * 
			1.2, ...)
	}
	if(sum(selection2) > 0) {
		points(x$long[selection2], x$lat[selection2], pch = 1, cex = 
			cex, ...)
		text(x$long[selection2], x$lat[selection2] + cxy[2] * cex *
			1.1, paste(" ", x$name[selection2], sep = ""), cex = 
			cex * 1.1, ...)
	}
	if(sum(selection3) > 0) {
		points(x$long[selection3], x$lat[selection3], pch = 1, cex = 
			cex, ...)
		text(x$long[selection3], x$lat[selection3] + cxy[2] * cex *
			0.9, paste(" ", x$name[selection3], sep = ""), cex = 
			cex * 0.9, ...)
	}
	invisible()
}

# draw a scale bar on a map
"map.scale" <-
function (x, y, relwidth = 0.15, metric = TRUE, ratio = TRUE, ...) 
{
  # old version
  format.pretty <- function(x) {
    as.character(pretty(x * c(0.99, 1.01), n = 2)[2])
  }
  # minka: new version
  format.pretty <- function(x, digits = 2) {
    x = signif(x, 2)
    prettyNum(formatC(x, format = "fg", digits = digits), big.mark = ",")
  }
  usr <- par("usr")
  if (missing(y)) 
    y <- (9 * usr[3] + usr[4])/10
  if (abs(y) >= 90) 
    warning("location of scale out of this world!")
  if (missing(x)) 
    #x <- (0.9 - relwidth) * usr[2] + (0.1 + relwidth) * usr[1]
    x <- (9 * usr[1] + usr[2])/10
  cosy <- cos((2 * pi * y)/360)
  perdeg <- (2 * pi * (6356.78 + 21.38 * cosy) * cosy)/360
  scale <- (perdeg * 100000)/(2.54 * (par("pin")/diff(par("usr"))[-2])[1])
  if (metric) 
    unit <- "km"
  else {
    perdeg <- perdeg * 0.6213712
    unit <- "mi"
  }
  len <- perdeg * relwidth * (usr[2] - usr[1])
  ats <- pretty(c(0, len), n = 2)
  nats <- length(ats)
  labs <- as.character(ats)
  labs[nats] <- paste(labs[nats], unit)
  linexy <- matrix(NA, ncol = 2, nrow = 3 * nats)
  colnames(linexy) <- c("x", "y")
  cxy <- par("cxy")
  dy <- cxy[2] * par("tcl")
  dx <- ats[nats]/perdeg/(nats - 1)
  linexy[1, ] <- c(x, y)
  linexy[2, ] <- c(x, y + dy)
  for (i in 1:(nats - 1)) {
    linexy[3 * i, ] <- c(x + (i - 1) * dx, y)
    linexy[3 * i + 1, ] <- c(x + i * dx, y)
    linexy[3 * i + 2, ] <- c(x + i * dx, y + dy)
  }
  lines(linexy)
  # minka: this is broken
  text(x + ats/perdeg, y + dy - 0.5 * cxy[2], labs, adj = c(0.4, 0.5), ...)
  # minka: added ratio option
  if(ratio)
    text(x, y + 0.5 * cxy[2],
         paste("scale approx 1:", format.pretty(scale), sep = ""),
         adj = 0, ...)
  invisible(scale)
}

map.wrap <- function(p) {
  # insert NAs to break lines that wrap around the globe.
  # does not work properly with polygons.
  # p is list of x and y vectors
  dx = abs(diff(p$x))
  dax = abs(diff(abs(p$x)))
  j = which(dx/dax > 50)
  j = c(j, length(p$x))
  start = 1
  x = c()
  y = c()
  for(i in j) {
    if(length(x) > 0) {
      x = c(x, NA)
      y = c(y, NA)
    }
    x = c(x, p$x[start:i])
    y = c(y, p$y[start:i])
    start = i + 1
  }
  list(x = x[2:length(x)], y = y[2:length(y)])
}
"map.old" <-
function (database = "state", regions = ".", exact = F, boundary = T, 
    interior = T, fill = F, projection = "", parameters = NULL, 
    orientation = rep(NA, 3), color = 1, add = F, plot = T, namesonly = F, 
    xlim = c(-1e+30, 1e+30), ylim = c(-1e+30, 1e+30), resolution = 1, 
    type = "l", ...) 
{
  "makepoly" <- function(xy, gonsize, keep) {
    x <- xy$x
    y <- xy$y
    n <- length(x)
    gonsize <- gonsize[ - length(gonsize)]
    discard <- seq(x)[is.na(x)]
    if(length(discard) > 0)
      discard <- c(discard, discard - 1, n)
    if(length(gonsize) > 0)
      discard <- discard[ - cumsum(gonsize)]
    if(length(discard) > 0) {
      x <- x[ - discard]
      y <- y[ - discard]
    }
    keep <- rep(keep, diff(c(0, seq(x)[is.na(x)], length(x))))
    list(x = x[keep], y = y[keep])
  }

    if (!missing(resolution) && !plot) 
        stop("must have plot=T if resolution is given")
    if (!fill && !boundary && !interior) 
        stop("one of boundary and interior must be TRUE")
    doproj <- !missing(projection) || !missing(parameters) || 
        !missing(orientation)
    coordtype <- maptype(database)
    if (coordtype == "unknown") 
        stop("missing database or unknown coordinate type")
    if (doproj && coordtype != "spherical") 
        stop(paste(database, "database is not spherical; projections not allowed"))
    gon <- mapname(database, regions, exact)
    n <- length(gon)
    if (n == 0) 
        stop("nothing to draw: no recognized region names")
    line <- mapgetg(database, gon, fill, xlim, ylim)
    if (length(line$number) == 0) {
        if (missing(xlim) || missing(ylim)) 
            stop("nothing to draw: all regions out of bounds")
        else coord <- list(x = c(xlim[1], NA, xlim[2]), y = c(ylim[1], 
            NA, ylim[2]), range = c(xlim, ylim))
        if (fill) 
            stop("cannot fill: all data out of bounds")
    }
    else {
        if (fill) 
            coord <- mapgetl(database, line$number, xlim, ylim)
        else {
            l <- abs(line$number)
            if (boundary && interior) 
                l <- unique(l)
            else if (boundary) 
                l <- l[!match(l, l[duplicated(l)], F)]
            else l <- l[duplicated(l)]
            coord <- mapgetl(database, l, xlim, ylim)
            if (length(coord) == 0) 
                if (missing(xlim) || missing(ylim)) 
                  stop("all data out of bounds")
                else coord <- list(x = c(xlim[1], NA, xlim[2]), 
                  y = c(ylim[1], NA, ylim[2]), range = c(xlim, 
                    ylim))
        }
    }
    if (doproj) {
        coord <- mapproject(coord, pr = projection, pa = parameters, 
            or = orientation)
        if (plot && coord$error) 
            if (all(is.na(coord$x))) 
                stop("projection failed for all data")
            else warning("projection failed for some data")
    }
    if (fill) {
        gonsize <- line$size
        color <- rep(color, length = length(gonsize))
        keep <- !is.na(color)
        coord[c("x", "y")] <- makepoly(coord, gonsize, keep)
        color <- color[keep]
    }
    if (plot) {
        if (!add) {
	    # remove margins
            par(pin = par("fin"))
            plot.new()
            #xrange <- coord$range[1:2]
            #yrange <- coord$range[3:4]
            xrange <- range(coord$x, na.rm = T)
            yrange <- range(coord$y, na.rm = T)
	    border <- c(0.01, 0.01)
            if (!missing(xlim)) {
                xrange[1] <- max(xrange[1], xlim[1])
                xrange[2] <- min(xrange[2], xlim[2])
		border[1] <- 0
            }
            if (!missing(ylim)) {
                yrange[1] <- max(yrange[1], ylim[1])
                yrange[2] <- min(yrange[2], ylim[2])
		border[2] <- 0
            }
            aspect <- if (coordtype != "spherical" || doproj) 
                c(1, 1)
            else c(cos((mean(yrange) * pi)/180), 1)
            d <- c(diff(xrange), diff(yrange)) * aspect
            p <- par("pin")
            # assumes coordinates are already clipped
            p <- d*min(p/d)
	    par(pin = p)
            d <- d*border + ((p/min(p/d) - d)/2)/aspect
            usr <- c(xrange, yrange) + rep(c(-1, 1), 2) * rep(d, c(2, 2))
            par(usr = usr)
        }
        if (resolution != 0 && type != "n") {
            pin <- par("pin")
            usr <- par("usr")
            resolution <- resolution * min(diff(usr)[-2]/pin/100)
	    coord[c("x", "y")] <- mapthin(coord, resolution)
        }
        if (type != "n") {
            #oerr <- par(err = -1)
            #on.exit(par(oerr))
            if (fill) 
                polygon(coord, col = color, ...)
            else lines(coord, col = color, type = type, ...)
        }
    }
    value <- if (namesonly) 
        line$name
    else coord[c("x", "y", "range")]
    if (plot) 
        invisible(value)
    else value
}
map.where <- function(database = "world", x, y)
{
  if(missing(y)) {
    if(is.matrix(x)) { y <- x[, 2]; x <- x[, 1] }
    else if(is.list(x) && !is.null(x$y)) { y <- x$y; x <- x$x }
    else { y <- x[[2]]; x <- x[[1]] }
  }
  if(is.character(database)) {
    dbname <- paste(database, "MapEnv", sep = "")
    data(list = dbname)
    mapbase <- paste(Sys.getenv(get(dbname)), database, sep = "")
    gon <- .C("map_where", PACKAGE="maps",
       as.character(mapbase),
       as.double(x),
       as.double(y),
       as.integer(length(x)),
       integer(length(x)))[[5]]
    # this must be database, not mapbase
    nam <- mapname(database, ".")
    gon[gon == 0] = NA
    names(nam)[gon]
  }
  else {
    if(num.polygons(database) != length(database$names))
      stop("map object must have polygons (fill=TRUE)")
    n = length(database$x)
    result = .C("map_in_polygon", PACKAGE="maps",
       as.double(database$x), as.double(database$y), as.integer(n),
       as.double(x), as.double(y), as.integer(length(x)),
       integer(length(x)), NAOK = TRUE)[[7]]
    result[result == 0] = NA
    database$names[result]
  }
}
as.matrix.polygon <- function(x) {
  p = x
  if(is.null(p)) return(p)
  if(is.list(p) && !is.data.frame(p)) p <- cbind(p$x, p$y)
  p
}
in.one.polygon <- function(p, x) {
  # returns a logical vector, whose length is nrow(x)
  if(is.null(p)) return(NA)
  p <- as.matrix.polygon(p)
  if(is.list(x) && !is.data.frame(x)) x <- cbind(x$x, x$y)
  if(is.vector(x)) dim(x) <- c(1, 2)
  # p and x are matrices
  .C("map_in_one_polygon", PACKAGE="maps",
     as.double(p[, 1]), as.double(p[, 2]), as.integer(nrow(p)),
     as.double(x[, 1]), as.double(x[, 2]), as.integer(nrow(x)),
     logical(nrow(x)), as.integer(TRUE))[[7]]
}
in.polygon <- function(p, x) {
  # returns a logical vector, whose length is nrow(x)
  if(is.null(p)) return(NA)
  p <- as.matrix.polygon(p)
  if(is.list(x) && !is.data.frame(x)) x <- cbind(x$x, x$y)
  if(is.vector(x)) dim(x) <- c(1, 2)
  # p and x are matrices
  .C("map_in_polygon", PACKAGE="maps",
     as.double(p[, 1]), as.double(p[, 2]), as.integer(nrow(p)),
     as.double(x[, 1]), as.double(x[, 2]), as.integer(nrow(x)),
     logical(nrow(x)), NAOK = TRUE)[[7]] > 0
}

# polygon is not assumed closed
area.polygon <- function(p) {
  if(is.null(p)) return(NA)
  p <- as.matrix.polygon(p)
  n <- nrow(p)
  x1 <- p[, 1]
  i2 <- c(n, 1:(n - 1))
  x2 <- p[i2, 1]
  y1 <- p[, 2]
  y2 <- p[i2, 2]
  0.5*abs(sum(x1*y2 - x2*y1))
}

centroid.polygon <- function(p) {
  if(is.null(p)) return(c(NA, NA))
  p <- as.matrix.polygon(p)
  n <- nrow(p)
  x1 <- p[, 1]
  i2 <- c(n, 1:(n - 1))
  x2 <- p[i2, 1]
  y1 <- p[, 2]
  y2 <- p[i2, 2]
  a <- x1*y2 - x2*y1
  s <- sum(a)*3
  if(s == 0) c(mean(x1), mean(y1))
  else c(sum((x1 + x2)*a)/s, sum((y1 + y2)*a)/s)
}

# applies fun to all sub-polygons of p
apply.polygon <- function(p, fun, names. = NULL) {
  if(is.null(p)) return(p)
  if(is.null(names.) && !is.null(p$names)) names. = p$names
  p = as.matrix.polygon(p)
  n <- nrow(p)
  breaks <- (1:n)[is.na(p[, 1])]
  starts <- c(1, breaks + 1)
  ends <- c(breaks - 1, n)
  m <- length(starts)
  result <- list()
  for(i in 1:m) {
    this.p = if(ends[i] >= starts[i]) p[starts[i]:ends[i], ] else NULL
    result[[i]] <- fun(this.p)
  }
  names(result) = names.
  result
}

num.polygons <- function(p) {
  if(is.list(p)) 1 + sum(is.na(p$x))
  else 1 + sum(is.na(p[, 1]))
}

range.polygon <- function(..., na.rm = FALSE) {
  p <- as.list.polygon(...)
  lapply(p[c("x", "y")], range, na.rm = na.rm)
}

map.text <- function(database, regions = ".", labels, cex = 0.75, add = FALSE,
                     move = FALSE, ...) {
  if(!add) map(database, regions, ...)
  # get polygons
  cc = match.call(expand.dots=TRUE)
  cc[[1]] = as.name("map")
  cc$fill = TRUE
  cc$plot = FALSE
  cc$move = cc$add = cc$cex = cc$labels = NULL
  cc$resolution = 0
  m = eval(cc)
  if(missing(labels)) {
    labels = gsub(".*,", "", m$names)
  }
  if(num.polygons(m) != length(labels))
    stop("map object must have polygons (fill=TRUE)")
  x = apply.polygon(m, centroid.polygon)
  # convert m into a matrix
  x <- t(array(unlist(x), c(2, length(x))))
  if(move) {
    library(mining)
    w = strwidth(labels, units = "inches", cex = cex)
    h = strheight(labels, units = "inches", cex = cex)
    x = move.collisions2(x[, 1], x[, 2], w, h)
  }
  # want to omit map-specific options here (like "exact")
  text(x, labels, cex = cex, ...)
  invisible(m)
}

identify.map <- function(x, n = 1, index = FALSE, ...) {
  # identify polygons in a map
  # must click near the center of the polygon
  m = x
  if(!is.list(m)) stop("must provide a map object")
  if(num.polygons(m) != length(m$names))
    stop("map object must have polygons (fill=TRUE)")
  x = apply.polygon(m, centroid.polygon)
  x <- t(array(unlist(x), c(2, length(x))))
  i = identify(x[, 1], x[, 2], labels = m$names, n = n, ...)
  if(index) i else m$names[i]
}

area.map <- function(m, regions = ".", sqmi=TRUE, ...) {
  # returns the areas of given regions,
  # combining the areas of all regions which match.
  if(!is.list(m)) stop("must provide a map object")
  if(num.polygons(m) != length(m$names))
    stop("map object must have polygons (fill=TRUE)")
  proj = m$projection
  m = map.poly(m,regions,as.polygon=TRUE,...)
  area = unlist(apply.polygon(m, area.polygon))
  merge <- regions[match.map(m, regions, ...)]
  names(merge) <- m$names
  merge = factor(merge, levels = regions)
  area = drop(indicators.factor(merge) %*% area)
  areaSqMiles <- function(proj) {
    # returns a factor f such that f*area.map() is in square miles.
    if(is.null(proj)) proj = "no projection"
    if(proj %in% c("mollweide","azequalarea","aitoff")) 2*15732635
    else if(proj %in% c("sinusoidal","bonne","cylequalarea","albers")) 15732635
    else if(proj == "sp_albers") 15745196
    else {
      warning(paste("sq.mile correction unavailable for",proj))
      1
    }
  }
  if(sqmi) area*areaSqMiles(proj) else area
}
indicators.factor <- function(y) {
  # convert a factor into a matrix of indicators
  # result is level by case
  # works if y contains NAs
  r <- array(0, c(length(levels(y)), length(y)), list(levels(y), names(y)))
  for(lev in levels(y)) r[lev, y == lev] <- 1
  r
}

# Functions for smoothing or dis-aggregating data over map regions.

# m is a map object
# z is a named vector
# res is resolution of sampling grid
# span is kernel parameter (larger = smoother)
#   span = Inf is a special case which invokes cubic spline kernel.
#   span is scaled by the map size, and is independent of res.
# result is a frame
smooth.map <- function(m, z, res = 50, span = 1/10, averages = FALSE,
                       type = c("smooth", "interp"), merge = FALSE) {
  #if(is.data.frame(z)) z = as.named.vector(z)
  if(averages) {
    # turn averages into sums
    z = z * area.map(m, names(z), sqmi=FALSE)
  }
  # sampling grid
  xlim <- range(m$x, na.rm = TRUE)
  ylim <- range(m$y, na.rm = TRUE)
  midpoints <- function(start, end, n) {
    inc <- (end - start)/n
    seq(start + inc/2, end - inc/2, len = n)
  }
  # 2*res is an assumption about aspect ratio (usually true)
  if(length(res) == 1) res = c(2*res, res)
  xs <- midpoints(xlim[1], xlim[2], res[1])
  ys <- midpoints(ylim[1], ylim[2], res[2])
  x <- expand.grid(x = xs, y = ys)
  if(FALSE) {
    # add centroids to the sample points
    xc = apply.polygon(m[c("x", "y")], centroid.polygon)
    # convert m into a matrix
    xc <- t(array(unlist(xc), c(2, length(xc))))
    xc = data.frame(x = xc[, 1], y = xc[, 2])
    x = rbind(x, xc)
  }
  radius = sqrt(diff(xlim)^2 + diff(ylim)^2)/2
  lambda = 1/(span*radius)^2
  #cat("lambda = ", lambda, "\n")
  cell.area = diff(xs[1:2])*diff(ys[1:2])

  r <- factor(map.where(m, x))
  if(merge) {
    # merge regions
    # merge[r] is the parent of region r
    # regions with merge[r] = NA are considered absent from the map
    # (no sample points will be placed there)
    # this can be slow on complex maps
    merge <- names(z)
    merge <- merge[match.map(m, merge)]
    names(merge) <- m$names
    levels(r) <- merge[levels(r)]
  }
  # remove points not on the map
  i <- !is.na(r)
  x <- x[i, ]
  r <- r[i]
  xo = x
  if(TRUE) {
    # kludge - drop regions with no samples
    n = table(r)
    bad = (n == 0)
    newlevels = levels(r)
    newlevels[bad] = NA
    levels(r) = newlevels
  }
  # put z in canonical order, and drop values which are not in the map
  z = z[levels(r)]
  # remove regions not named in z, or where z is NA
  bad = is.na(z)
  z = z[!bad]
  newlevels = levels(r)
  newlevels[bad] = NA
  levels(r) = newlevels
  i <- !is.na(r)
  x <- x[i, ]
  r <- r[i]
  # do all regions have sample points?
  n = table(r)
  if(any(n == 0)) stop(paste(paste(names(n)[n == 0], collapse = ", "), "have no sample points"))
  type <- match.arg(type)
  if(FALSE) {
    # code for these is in 315/util.r
    # most time is spent here
    w <- switch(type,
                mass = gp.smoother(x, x, r, lambda),
                smooth = kr.smoother(x, x, r, lambda))
    #list(x = x, r = r, w = w)
    z = drop(z %*% w)
    cbind(x, z = z)
  } else {
    if(type == "smooth") {
      z = kernel.smooth(x, z, xo, lambda, r)
    } else {
      z = gp.smooth(x, z, xo, lambda, r)
    }
    z = z/cell.area
    cbind(xo, z = z)
  }
}

gp.smooth <- function(x, z, xo, lambda, r) {
  # predict a function measured at locations x to new locations xo
  krr = kernel.region.region(x, r, lambda)
  white.z = solve(krr, z)
  kernel.smooth(x, white.z, xo, lambda, r, normalize = FALSE)
}

kernel.smooth <- function(x, z, xo, lambda, region = NULL, normalize = TRUE) {
  # predict a function measured at locations x to new locations xo
  if(!is.matrix(x)) dim(x) <- c(length(x), 1)
  if(!is.matrix(xo)) dim(xo) <- c(length(xo), 1)
  n = nrow(x)
  if(is.null(region)) region = 1:n
  if(length(region) < n) stop("region must have same length as x")
  region = as.integer(region)
  if(any(is.na(region))) stop("region has NAs")
  if(max(region) > length(z)) stop("not enough measurements for the number of regions")
  no = nrow(xo)
  if(normalize) {
    # divide by region sizes
    z = as.double(z/as.numeric(table(region)))
  }
  .C("kernel_smooth", PACKAGE="maps",
     as.integer(n), as.integer(ncol(x)),
     as.double(t(x)), z, as.integer(region),
     as.integer(no), as.double(t(xo)), zo = double(no),
     as.double(lambda), as.integer(normalize))$zo
}

kernel.region.region <- function(x, region, lambda) {
  if(!is.matrix(x)) dim(x) <- c(length(x), 1)
  region = as.integer(region)
  nr = max(region)
  krr = .C("kernel_region_region", PACKAGE="maps",
    as.integer(nrow(x)), as.integer(ncol(x)),
    as.double(t(x)),
    region, as.double(lambda), as.integer(nr), krr = double(nr*nr))$krr
  dim(krr) = c(nr, nr)
  krr
}
kernel.region.x <- function(x, region, z, lambda) {
  if(!is.matrix(x)) dim(x) <- c(length(x), 1)
  if(!is.matrix(z)) dim(z) <- c(length(z), 1)
  region = as.integer(region)
  nr = max(region)
  no = nrow(z)
  krx = .C("kernel_region_x", PACKAGE="maps",
    as.integer(nrow(x)), as.integer(ncol(x)),
    as.double(t(x)), region, as.integer(no), as.double(t(z)),
    as.double(lambda), as.integer(nr), krx = double(nr*no))$krx
  dim(krx) = c(nr, no)
  krx
}
.First.lib <- function(lib, pkg) {
  # minka: only do this if R_MAP_DATA_DIR doesn't exist?
  if (Sys.getenv("R_MAP_DATA_DIR") == "")
    Sys.putenv("R_MAP_DATA_DIR"=paste(lib, pkg, "mapdata/", sep="/"))
  library.dynam("maps", pkg, lib)
}
