[6317] in s-news-athena
Re: maps
daemon@ATHENA.MIT.EDU (John Wallace)
Mon Jan 30 17:40:02 1995
Date: Mon, 30 Jan 1995 14:02:57 -0800 (PST)
From: John Wallace <jrw@fish.washington.edu>
To: Greg Unrau <gregu@ocgy.ubc.ca>
Cc: spluslist <S-news@stat.wisc.edu>
In-Reply-To: <9501192121.AA00813@penguin.ocgy.ubc.ca>
On Tue, 17 Jan 1995, Greg Unrau wrote:
> Well, strange is an understatement. I am using the maps library available
> from statlib and am wanting to extract a portion of Canada from the 'canada'
> map. I gave lat and lon limits to the map argument and it seems to do ok
> (i.e no smoke or keyboard jammings). The issue is what is plotted. Instead of
> a nicely filled map I Have straight lines running through the middle of the
> region and on one side it's filled and on the other it is blank. Any suggestions?
>
> Greg
Sorry for the slow reply, I have been on medical leave with a new baby boy.
The following is a repost from 24 June '94.
-----------------------------------------------
I have added a option 'east.to.west' to the map function, when set to
TRUE this option subtracts 360 from all the eastern longitudes. This
changes for example, 175E to -185, so that it maps on the left side
correctly. The program automatically sets 'east.to.west' under the
following conditions:
if((database == "world" | database == "world.thin") & (regions ==
"usa" | regions == "*alaska")) east.to.west <- T
For those interested the following did the trick:
if(east.to.west)
coord$x[is.finite(coord$x) & coord$x > 0] <-
coord$x[is.finite(coord$x) & coord$x > 0] - 360
I have appended a dump of the changed map function below.
Cheers,
--
-John Wallace jw@u.washingtion.edu
-Biostatistician, Splus Guru Level of Enlightenment, 7+ (log scale)
-Fisheries Research Institute, University of Washington
-phone: (206) 543-1513 fax: (206) 685-7471
On Wed, 22 Jun 1994, Jeff Simonoff wrote:
>
> About 9 months ago a question was raised regarding the apparent inability
> to construct a map of the world (using the "Maps in S" functions map()
> and mapproject()) centered at the 180 degree longitudinal meridian (the
> Pacific Ocean). As I recall, no solution to this problem was given.
> Has anyone come up with a solution, either in changes to the functions
> or to the database?
>
> Thanks for any help anyone can give ...
>
> Jeff Simonoff
> NYU Dept. of Statistics & OR
> jsimonoff@stern.nyu.edu
>
"map"<-
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",
east.to.west = F, ...)
{
# REVSION: Option 'east.to.west' Added
# REVISED BY: John R. Wallace (jw@u.washington.edu)
# LAST REVISED: 2 Feb '94
# parameter checks
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"))
# turn the region names into a list of polygon numbers
gon <- mapname(database, regions, exact)
n <- length(gon)
if(n == 0) stop("nothing to draw: no recognized region names")
# turn the polygon numbers into a list of polyline numbers
line <- mapgetg(database, gon, fill, xlim, ylim)
if(length(line$number) == 0)
stop("nothing to draw: all regions out of bounds")
if((database == "world" | database == "world.thin") & (regions ==
"usa" | regions == "*alaska")) east.to.west <- T
# turn the polyline numbers into x and y coordinates
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)
stop("all data out of bounds")
}
if(east.to.west)
coord$x[is.finite(coord$x) & coord$x > 0] <- coord$x[is.finite(
coord$x) & coord$x > 0] - 360
coord$range[1:2] <- range(coord$x, na.rm = T)
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")
}
# for filled regions, turn NA breaks at polylines into
# NA breaks at polygons, deleting polygons for which
# there is a corresponding NA color
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]
}
# 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 sqaure appears
# square in the middle of the plot
if(!add) {
xrange <- coord$range[1:2]
yrange <- coord$range[3:4]
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")
d <- d/100 + ((p/min(p/d) - d)/2)/aspect
par(usr = c(xrange, yrange) + rep(c(-1, 1), 2) * rep(d,
c(2, 2)))
if(!par("new"))
frame()
}
# do thinning
if(resolution != 0 && type != "n") {
uin <- par("uin")
rsz <- par("rsz")
resolution <- resolution * min(rsz[1]/uin[1], rsz[2]/
uin[2])
coord[c("x", "y")] <- mapthin(coord, resolution)
}
# suppress warnings about clipping
if(type != "n") {
oerr <- par(err = -1)
on.exit(par(oerr))
if(fill)
polygon(coord, col = color, ...)
else lines(coord, col = color, type = type, ...)
}
}
# return value is names or coords, but not both
value <- if(namesonly) line$name else coord[c("x", "y", "range")]
if(plot)
invisible(value)
else value
}
