SUMMARY: code to modify world map split

["Peter Ward" <ward@mathstat.dal.ca>]

I received many helpful suggestions on how to modify the world map so that
it is split at another longitude. I eventually used the following code,
which simply adds the lower longitudes to the upper longitudes. Of course,
the new lines are  not continuous for each land mass, so you cannot plot the
map as polygons (only plot as lines).

Many thanks to the S-NEWS list for the helpful advice.



Peter Ward
Dept. of Biology
Dalhousie University   ph.      (902) 494 3910
Halifax B3H 4J1         fax       (902) 494 3736
CANADA                 e-mail: ward@mathstat.dal.ca
You might try something like this:

> mymap <- map("world",xlim=c(20,180))
> mymap2  <- map("world",xlim=c(-180,19.99))
> mymap2$x <- mymap2$x+360
> mymap$x <- c(mymap$x,NA,mymap2$x)
> mymap$y <- c(mymap$y,NA,mymap2$y)
> mymap$range <- c(20,360+20)
> plot(mymap,type="l",xlim=c(20,380),axes="n")

 [mailto:thacker@aoml.noaa.gov]

********************************************************
The way I did it, to create the world2 databases, was to go back to the
original data, and modify each longitude point.  So instead of the range
being (-pi, pi), I changed it to (0, 2pi).  This isn't just a matter of
adding pi (otherwise the whole thing will still look the same!).  You
would have to add 2pi to all values less than 20 degrees.  Then of
course you have to recreate the binary form of the database, and then
generate the S-Plus format of this (no trivial task, but I have a
program to do this).

OR ... you can emulate this by something like:

library(maps)
world <- map("world", plot=F)
world$x[!is.na(world$x) & world$x < 20] <- world$x[!is.na(world$x) & world$x
< 20] + 360
motif()
plot(world, type="n")
lines(world)

And then you realise that you have to take care of the new breaks in the
data!.

A quick and dirty fix is:

world$y[!is.na(world$x) & abs(diff(world$x)) > 300] <- NA
world$x[!is.na(world$x) & abs(diff(world$x)) > 300] <- NA
plot(world, type="n")
lines(world)

ray@mcs.vuw.ac.nz
********************************************************
I remember the problem now:  the lines across the cut point are not clipped
properly, they go from edge to edge.  One way around this is to get rid of
all long lines between points.  A quick and dirty way to do this is:

map.omit.wraparound <-
function(xy, threshhold)
{
        # find places where line from xy[i] to xy[i+1] appears to
        # wrap around page and insert NA there.
        d <- abs(diff(xy$x))
        i <- !is.na(d) & d > threshhold
        i <- cumsum(c(1, i + 1))
        val <- list(x = rep(as.numeric(NA), length(i)), y = rep(as.numeric(
                NA), length(i)))
        val$x[i] <- xy$x
        val$y[i] <- xy$y
        val
}

and edit map to call it, if you wish, with the following change.  It helps
with the example
        map("state",proj="mercator",orientation=c(0,90,0))

*** /tmp/map.old        Wed Nov 20 12:59:23 2002
--- /tmp/map.new        Wed Nov 20 12:59:15 2002
***************
*** 2,8 ****
  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", ...)
  {
        # parameter checks
        if(!missing(resolution) && !plot) stop(
--- 2,9 ----
  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",
!       omit.wraparound = T, ...)
  {
        # parameter checks
        if(!missing(resolution) && !plot) stop(
***************
*** 84,89 ****
--- 85,98 ----
                                2]/uin[2])
                        coord[c("x", "y")] <- mapthin(coord, resolution)
                }
+               if(omit.wraparound) {
+                       # experimental: get rid of wraparound.  Assume any
line
+                       # covering more than 90% of width of plot should be
going around
+                       # back of map.  The 90% should be another argument
to map.
+                       coord[c("x", "y")] <- map.omit.wraparound(coord,
+                               threshhold = diff(range(coord$x, na.rm = T))
*
+                               0.9)
+               }
                # suppress warnings about clipping
                if(type != "n") {
                        oerr <- par(err = -1)


bill@insightful.com
********************************************************
You can use the coastline tables that I put in statlib long ago, and plot
the map any way you want, divided anywhere you want.  But then there is
the matter of projections....  There are no geopolitical boundaries.

wofsy@fas.harvard.edu