[R-sig-Geo] map projection problem

[Roger Bivand]

On Tue, 1 Apr 2008, Melanie Abecassis wrote:

> 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