Skip to content

map projection problem

2 messages · Melanie Abecassis, Roger Bivand

Original
Melanie Abecassis
# 
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.


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
Roger Bivand
#
On Tue, 1 Apr 2008, Melanie Abecassis wrote:

            

      
      
      
    

        
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