R/map.r

In adeckmyn/maps: Draw Geographical Maps

# 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 <- which(is.na(x))
  if (length(breaks) == 0) {
    starts <- 1
    ends <- n
  } else {
    starts <- c(1, breaks + 1)
    ends <- c(breaks - 1, n)
  }
## AD: EXTRA: allow negative values in i reverse direction
  pl <- function(j) if(j>0) x[starts[j]:ends[j]] else x[ends[abs(j)]:starts[abs(j)]]
  as.numeric(unlist(lapply(i,function(j) c(NA,pl(j))))[-1])
}

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

# FUTURE:
#sub.gondata <- function(p, i) {
#  x <- p$gondata
#  n <- length(x)
#  breaks <- which(is.na(x))
#  if (length(breaks) == 0) {
#    starts <- 1
#    ends <- n
#  } else {
#    starts <- c(1, breaks + 1)
#    ends <- c(breaks - 1, n)
#  }
#  lengths <- starts - ends + 1
#  i2 <- unlist(lapply(i,function(j) c(NA,x[starts[j]:ends[j]])))[-1]
#  xy <- sub.polygon(p, i2)
#}

# 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, namefield="name") {
  if (!is.character(database)) {
    if (!as.polygon) stop("map objects require as.polygon=TRUE")
    if (inherits(database,"Spatial")){
      if (inherits(database,"SpatialPolygons")) the.map <- SpatialPolygons2map(database, 
                                                               namefield=namefield)
      else if (inherits(database,"SpatialLines")) the.map <- SpatialLines2map(database, 
                                                             namefield=namefield)
      else stop("database not supported.")
    } else 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 {
## AD TODO: solve the problem for UK vs Ukrain...
## But here you have a simple set of named polygons, not a map database.
## If it comes from anything but "world", we could break something else.
## We just leave it for now.
        regexp <- paste("(^", regions, ")", sep = "", collapse = "|")
## FIXME: find a better solution
# perl regex crashes on >30000 characters, e.g. whole world database
        i <- grep(regexp, the.map$names, ignore.case = TRUE,
                  perl = (length(regions) < 1000) & (nchar(regexp) < 30000) )
      }
      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.range <-
local({
   val <- NULL
   function(new) if(!missing(new)) val <<- new else val
})

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(4.1, 4.1, par("mar")[3], 0.1), myborder = 0.01, 
         namefield="name", lforce = "n", ...)
{
  # 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)
  if (doproj && !requireNamespace("mapproj", quietly=TRUE)) {
    stop("Please install the package 'mapproj' for projections.")
  }
  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"))
  if (length(wrap)>=2 && !doproj && wrap[2] - wrap[1] != 360) 
    stop("The specified longitudes for wrapping are inconsistent, they should be 360 apart.")
  # turn the region names into x and y coordinates
  if (is.character(database)) as.polygon = fill
  else as.polygon <- TRUE
  # if we are going to wrap around anything else than [-180, 180], the simplest way
  # to get all necessary polylines/polygons is to set xlim=NULL *temporarily*
  # (alternatively, the C code must shift longitudes by +/- 360 to fit in boundaries)
  xlim_tmp <- if (length(wrap)>=2) NULL else xlim
  coord <- map.poly(database, regions, exact, xlim_tmp, ylim, 
                    boundary, interior, fill, as.polygon, namefield=namefield)
  if (is.na(coord$x[1])) stop("first coordinate is NA. Bad map data?")
  if (length(wrap)>=2) {
    antarctica <- if (length(wrap) == 2) -89.9 else wrap[3]
    coord <- map.wrap.poly(coord, xlim=wrap[1:2], poly=fill, antarctica=antarctica)
  }
  # we can enforce xlim & ylim exactly
  # this changes the output
  if (lforce=="e") {
    coord <- map.clip.poly(coord, xlim, ylim, poly=fill)
  }
  if (plot) {
    .map.range(coord$range)
  }

  if (doproj) {
    nam <- coord$names
    coord <- mapproj::mapproject(coord, projection = projection,
			parameters = parameters, orientation = orientation)
    coord$projection = projection
    coord$parameters = parameters
    coord$orientation = orientation
    coord$names <- nam
    if (!is.null(xlim) && !is.null(ylim) && lforce %in% c("s","l")) {
      prange <- mapproj::mapproject(x=rep(xlim,2), y=rep(ylim, each=2))
      if (lforce=="s") {
        xlim <- c(max(prange$x[c(1,3)]),min(prange$x[c(2,4)]))
        ylim <- c(max(prange$y[c(1,2)]),min(prange$y[c(3,4)]))
      } else {
        xlim <- c(min(prange$x[c(1,3)]),max(prange$x[c(2,4)]))
        ylim <- c(min(prange$y[c(1,2)]),max(prange$y[c(3,4)]))
      }
    }
    if (plot && coord$error)
      if (all(is.na(coord$x)))
        stop("projection failed for all data")
      else warning("projection failed for some data")
  }
  # AD: we do wrapping first: slightly better than when run after the thinning
  #     also now the output data is also wrapped if plot=FALSE
  if (length(wrap)==1 && wrap) coord <- map.wrap(coord)
  # 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 unless we have pforce=TRUE
      if (is.null(xlim) || (doproj && !(lforce %in% c("s","l")))) xrange <- range(coord$x, na.rm = TRUE)
      else xrange <- xlim
      if (is.null(ylim) || (doproj && !(lforce %in% c("s","l")))) 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 * myborder) * 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 ??!
	par(mar = mar)	# set mai if user-defined mar
	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 * myborder + ((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))
    }
    if (type != "n") {
      # thinning only works correctly if you have polylines from a database
      if (!as.polygon && resolution != 0) {
        pin <- par("pin")
        usr <- par("usr")
        resolution <- resolution * min(diff(usr)[-2]/pin/100)
        coord[c("x", "y")] <- mapthin(coord, resolution)
      }
      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
}