map projection problem

Hi,
I am unsuccessfully trying to change the map projection  in map function
of the maps library.

Here is my code :

library(maps)
library(mapdata)
library(mapproj)

map(database = "world", fill = TRUE, col = 1, plot = TRUE,add=F,
xlim = c(-200,-120), ylim = c(-5,45))
map.axes()

This gives me a map centered on Hawaii.


I now want to have *the same map but in the mercator projection* :

map(database = "world", fill = TRUE, col = 1, plot = TRUE,add=F,
projection="mercator", xlim = c(-200,-120), ylim = c(-5,45))
map.axes()

And I get a weird map of the US  and parts of Canada only.

When fill=TRUE, map() actually returns the whole polygon, so you are 
seeing what you asked map() for, that's just the way it works. If you 
don't need fill, it dangles a bit, but only gives you a little more than 
you asked for.

There are two alternatives, one using map which leaves you with the 
continental US and Canada (I'm going to -180 rather than -200):

library(maps)
library(mapproj)
m0 <- map(database = "world", plot = FALSE, fill=TRUE,
  xlim = c(-200,-120), ylim = c(-5,45))
IDs <- sapply(strsplit(m0$names, ":"), function(x) x[1])
library(maptools)
m1 <- map2SpatialPolygons(m0, IDs=IDs, proj4string = CRS("+proj=longlat"))
summary(m1)
library(rgdal)
m2 <- spTransform(m1, CRS("+proj=merc"))
bb <- cbind(c(-180,-120), c(-5,45))
bbSP <- SpatialPoints(bb, proj4string = CRS("+proj=longlat"))
bbmerc <- spTransform(bbSP, CRS("+proj=merc"))
plot(m2, col=1, xlim=bbc[,1], ylim=bbc[,2])
summary(m2)

which gets you there but isn't clipped, or:

library(maptools)
gshhs.c.b <- system.file("share/gshhs_c.b", package="maptools")
g0 <- Rgshhs(gshhs.c.b, xlim=c(180, 240), ylim=c(-5, 45), level=1)
summary(g0$SP)
library(rgdal)
g2 <- spTransform(g0$SP, CRS("+proj=merc"))
plot(g2, col=1)
summary(g2)

using the coarse version of the GSHHS shorelines shipped with maptools. 
Both shrink the bounding box of the object to the actual shorelines.

In both of these cases, the output Mercator coordinates are the true 
coordinates, not the display coordinates used by mapproj.

Hope this helps,

Roger

Then I tried this :

A=mapproject(c(-200,-120), c(-5,45), projection="mercator")

map(database = "world", fill = TRUE, col = 1, plot = TRUE,add=F,
projection="mercator", xlim = A$x, ylim = A$y)
map.axes()

which should be the correct way to do this I think, but I get this :
Error in map.poly(database, regions, exact, xlim, ylim, boundary,
interior,  :
 nothing to draw: all regions out of bounds

Could anyone help me with this ?
Thanks a lot,
Melanie

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no