1# R data structure for maps is list(x, y, names)
2# names is character vector naming the polygons
3# (contiguous groups of non-NA coordinates)
4
5# returns the jth contiguous section of non-NA values in vector x
6# should consecutive NA's count as one?
7subgroup <- function(x, i) {
8  n <- length(x)
9  breaks <- which(is.na(x))
10  if (length(breaks) == 0) {
11    starts <- 1
12    ends <- n
13  } else {
14    starts <- c(1, breaks + 1)
15    ends <- c(breaks - 1, n)
16  }
17## AD: EXTRA: allow negative values in i reverse direction
18  pl <- function(j) if(j>0) x[starts[j]:ends[j]] else x[ends[abs(j)]:starts[abs(j)]]
19  as.numeric(unlist(lapply(i,function(j) c(NA,pl(j))))[-1])
20}
21
22sub.polygon <- function(p, i) {
23  lapply(p[c("x", "y")], function(x) subgroup(x, i))
24}
25
26# FUTURE:
27#sub.gondata <- function(p, i) {
28#  x <- p$gondata
29#  n <- length(x)
30#  breaks <- which(is.na(x))
31#  if (length(breaks) == 0) {
32#    starts <- 1
33#    ends <- n
34#  } else {
35#    starts <- c(1, breaks + 1)
36#    ends <- c(breaks - 1, n)
37#  }
38#  lengths <- starts - ends + 1
39#  i2 <- unlist(lapply(i,function(j) c(NA,x[starts[j]:ends[j]])))[-1]
40#  xy <- sub.polygon(p, i2)
41#}
42
43# returns a sub-map of the named map corresponding to the given regions
44# regions is a vector of regular expressions to match to the names in the map
45# regions outside of xlim, ylim may be omitted
46map.poly <- function(database, regions = ".", exact = FALSE,
47                     xlim = NULL, ylim = NULL, boundary = TRUE,
48		     interior = TRUE, fill = FALSE, as.polygon = FALSE, namefield="name") {
49  if (!is.character(database)) {
50    if (!as.polygon) stop("map objects require as.polygon=TRUE")
51    if (inherits(database,"Spatial")){
52      if (inherits(database,"SpatialPolygons")) the.map <- SpatialPolygons2map(database,
53                                                               namefield=namefield)
54      else if (inherits(database,"SpatialLines")) the.map <- SpatialLines2map(database,
55                                                             namefield=namefield)
56      else stop("database not supported.")
57    } else the.map <- database
58    if (identical(regions,".")) {
59      # speed up the common case
60      nam = the.map$names
61      coord <- the.map[c("x", "y")]
62    } else {
63      # same as mapname()
64      if (exact) {
65        i = match(regions, the.map$names)
66        if (any(is.na(i))) i = NULL
67      } else {
68## AD TODO: solve the problem for UK vs Ukrain...
69## But here you have a simple set of named polygons, not a map database.
70## If it comes from anything but "world", we could break something else.
71## We just leave it for now.
72        regexp <- paste("(^", regions, ")", sep = "", collapse = "|")
73## FIXME: find a better solution
74# perl regex crashes on >30000 characters, e.g. whole world database
75        i <- grep(regexp, the.map$names, ignore.case = TRUE,
76                  perl = (length(regions) < 1000) & (nchar(regexp) < 30000) )
77      }
78      if (length(i) == 0) stop("no recognized region names")
79      nam <- the.map$names[i]
80      coord <- sub.polygon(the.map, i)
81    }
82    coord$range <- c(range(coord$x, na.rm = TRUE), range(coord$y, na.rm = TRUE))
83  } else {
84    # turn the polygon numbers into a list of polyline numbers
85    gon <- mapname(database, regions, exact)
86    n <- length(gon)
87    if (n == 0) stop("no recognized region names")
88    if (is.null(xlim)) xlim <- c(-1e+30, 1e+30)
89    if (is.null(ylim)) ylim <- c(-1e+30, 1e+30)
90    # turn the polygon numbers into a list of polyline numbers
91    line <- mapgetg(database, gon, as.polygon, xlim, ylim)
92    if (length(line$number) == 0)
93            stop("nothing to draw: all regions out of bounds")
94    # turn the polyline numbers into x and y coordinates
95    if (as.polygon) {
96      coord <- mapgetl(database, line$number, xlim, ylim, fill)
97      # assemble lines into polygons
98      gonsize <- line$size
99      keep <- rep(TRUE, length(gonsize))
100      coord[c("x", "y")] <- makepoly(coord, gonsize, keep)
101    }
102    else {
103      l <- abs(line$number)
104      if (boundary && interior) l <- unique(l)
105      else if (boundary) l <- l[!match(l, l[duplicated(l)], FALSE)]
106      else l <- l[duplicated(l)]
107      coord <- mapgetl(database, l, xlim, ylim, fill)
108      if (length(coord) == 0)
109              stop("all data out of bounds")
110    }
111    nam <- line$name
112  }
113  list(x = coord$x, y = coord$y, range = coord$range, names = nam)
114}
115
116.map.range <-
117local({
118   val <- NULL
119   function(new) if(!missing(new)) val <<- new else val
120})
121
122map <-
123function(database = "world", regions = ".", exact = FALSE,
124	 boundary = TRUE, interior = TRUE, projection = "",
125	 parameters = NULL, orientation = NULL, fill = FALSE,
126	 col = 1, plot = TRUE, add = FALSE, namesonly = FALSE,
127         xlim = NULL, ylim = NULL, wrap = FALSE,
128         resolution = if (plot) 1 else 0, type = "l", bg = par("bg"),
129         mar = c(4.1, 4.1, par("mar")[3], 0.1), myborder = 0.01,
130         namefield="name", lforce = "n", ...)
131{
132  # parameter checks
133  if (resolution>0 && !plot)
134    stop("must have plot=TRUE if resolution is given")
135  if (!fill && !boundary && !interior)
136    stop("one of boundary and interior must be TRUE")
137  doproj <- !missing(projection) || !missing(parameters) || !missing(
138          orientation)
139  if (doproj && !requireNamespace("mapproj", quietly=TRUE)) {
140    stop("Please install the package 'mapproj' for projections.")
141  }
142  coordtype <- maptype(database)
143  if (coordtype == "unknown")
144     stop("missing database or unknown coordinate type")
145  if (doproj && coordtype != "spherical")
146    stop(paste(database, "database is not spherical; projections not allowed"))
147  if (length(wrap)>=2 && !doproj && wrap[2] - wrap[1] != 360)
148    stop("The specified longitudes for wrapping are inconsistent, they should be 360 apart.")
149  # turn the region names into x and y coordinates
150  if (is.character(database)) as.polygon = fill
151  else as.polygon <- TRUE
152  # if we are going to wrap around anything else than [-180, 180], the simplest way
153  # to get all necessary polylines/polygons is to set xlim=NULL *temporarily*
154  # (alternatively, the C code must shift longitudes by +/- 360 to fit in boundaries)
155  xlim_tmp <- if (length(wrap)>=2) NULL else xlim
156  coord <- map.poly(database, regions, exact, xlim_tmp, ylim,
157                    boundary, interior, fill, as.polygon, namefield=namefield)
158  if (is.na(coord$x[1])) stop("first coordinate is NA. Bad map data?")
159  if (length(wrap)>=2) {
160    antarctica <- if (length(wrap) == 2) -89.9 else wrap[3]
161    coord <- map.wrap.poly(coord, xlim=wrap[1:2], poly=fill, antarctica=antarctica)
162  }
163  # we can enforce xlim & ylim exactly
164  # this changes the output
165  if (lforce=="e") {
166    coord <- map.clip.poly(coord, xlim, ylim, poly=fill)
167  }
168  if (plot) {
169    .map.range(coord$range)
170  }
171
172  if (doproj) {
173    nam <- coord$names
174    coord <- mapproj::mapproject(coord, projection = projection,
175			parameters = parameters, orientation = orientation)
176    coord$projection = projection
177    coord$parameters = parameters
178    coord$orientation = orientation
179    coord$names <- nam
180    if (!is.null(xlim) && !is.null(ylim) && lforce %in% c("s","l")) {
181      prange <- mapproj::mapproject(x=rep(xlim,2), y=rep(ylim, each=2))
182      if (lforce=="s") {
183        xlim <- c(max(prange$x[c(1,3)]),min(prange$x[c(2,4)]))
184        ylim <- c(max(prange$y[c(1,2)]),min(prange$y[c(3,4)]))
185      } else {
186        xlim <- c(min(prange$x[c(1,3)]),max(prange$x[c(2,4)]))
187        ylim <- c(min(prange$y[c(1,2)]),max(prange$y[c(3,4)]))
188      }
189    }
190    if (plot && coord$error)
191      if (all(is.na(coord$x)))
192        stop("projection failed for all data")
193      else warning("projection failed for some data")
194  }
195  # AD: we do wrapping first: slightly better than when run after the thinning
196  #     also now the output data is also wrapped if plot=FALSE
197  if (length(wrap)==1 && wrap) coord <- map.wrap(coord)
198  # do the plotting, if requested
199  if (plot) {
200    # for new plots, set up the coordinate system;
201    # if a projection was done, set the aspect ratio
202    # to 1, else set it so that a long-lat square appears
203    # square in the middle of the plot
204    if (!add) {
205      opar = par(bg = bg)
206      if (!par("new")) plot.new()
207      # xlim, ylim apply before projection unless we have pforce=TRUE
208      if (is.null(xlim) || (doproj && !(lforce %in% c("s","l")))) xrange <- range(coord$x, na.rm = TRUE)
209      else xrange <- xlim
210      if (is.null(ylim) || (doproj && !(lforce %in% c("s","l")))) yrange <- range(coord$y, na.rm = TRUE)
211      else yrange <- ylim
212      if (coordtype != "spherical" || doproj) {
213	aspect <- c(1, 1)
214      } else
215        aspect <- c(cos((mean(yrange) * pi)/180), 1)
216      d <- c(diff(xrange), diff(yrange)) * (1 + 2 * myborder) * aspect
217      if (coordtype != "spherical" || doproj) {
218        plot.window(xrange, yrange, asp = 1/aspect[1])
219      } else {
220        # must have par(xpd = FALSE) for limits to have an effect ??!
221	par(mar = mar)	# set mai if user-defined mar
222	p <- par("fin") -
223	  as.vector(matrix(c(0, 1, 1, 0, 0, 1, 1, 0), nrow = 2) %*% par("mai"))
224	par(pin = p)
225        p <- par("pin")
226        p <- d * min(p/d)
227        par(pin = p)
228        d <- d * myborder + ((p/min(p/d) - d)/2)/aspect
229        usr <- c(xrange, yrange) + rep(c(-1, 1), 2) * rep(d, c(2, 2))
230        par(usr = usr)
231      }
232      on.exit(par(opar))
233    }
234    if (type != "n") {
235      # thinning only works correctly if you have polylines from a database
236      if (!as.polygon && resolution != 0) {
237        pin <- par("pin")
238        usr <- par("usr")
239        resolution <- resolution * min(diff(usr)[-2]/pin/100)
240        coord[c("x", "y")] <- mapthin(coord, resolution)
241      }
242      if (fill) polygon(coord, col = col, ...)
243      else lines(coord, col = col, type = type, ...)
244    }
245  }
246  # return value is names or coords, but not both
247  class(coord) = "map"
248  value <- if (namesonly) coord$names else coord
249  if (plot) invisible(value)
250  else value
251}
252
253"makepoly" <-
254function(xy, gonsize, keep)
255{
256  # remove NAs and duplicate points so that a set of polylines becomes a set
257  # of polygons.
258  # xy is a set of polylines, separated by NAs.
259  # gonsize is a vector, giving the number of lines to put into each polygon
260  # note that a polyline may consist of a single point
261  x <- xy$x
262  y <- xy$y
263  n <- length(x)
264  gonsize <- gonsize[ - length(gonsize)]
265  discard <- seq(length(x))[is.na(x)]
266  if (length(discard) > 0) {
267    # locations of (possible) duplicate points
268    dups = c(discard - 1, n)
269    # only polylines with > 1 point have duplicates
270    i = which(diff(c(0, dups)) > 2);
271    discard <- c(discard, dups[i]);
272  }
273  # first part of discard is the NAs, second part is duplicates
274  # gonsize tells us which NAs to preserve
275  if (length(gonsize) > 0)
276    discard <- discard[ - cumsum(gonsize)]
277  if (length(discard) > 0) {
278    x <- x[ - discard]
279    y <- y[ - discard]
280  }
281  keep <- rep(keep, diff(c(0, seq(length(x))[is.na(x)], length(x))))
282  closed.polygon(list(x = x[keep], y = y[keep]))
283}
284closed.polygon <- function(p) {
285  # p is a set of polylines, separated by NAs
286  # for each one, the first point is copied to the end, giving a closed polygon
287  x = p$x
288  y = p$y
289  n = length(x)
290  breaks <- seq(length(x))[is.na(x)]
291  starts <- c(1, breaks + 1)
292  ends <- c(breaks - 1, n)
293  x[ends + 1] = x[starts]
294  y[ends + 1] = y[starts]
295  x = insert(x, breaks + 1)
296  y = insert(y, breaks + 1)
297  list(x = x, y = y)
298}
299insert <- function(x, i, v = NA) {
300  # insert v into an array x, at positions i
301  # e.g. insert(1:7, c(2, 5, 8))
302  n = length(x)
303  new.n = n + length(i)
304  m = logical(new.n)
305  i.new = i - 1 + seq(length(i))
306  m[i.new] = TRUE
307  x = x[(1:new.n) - cumsum(m)]
308  x[i.new] = v
309  x
310}
311