[R-sig-Geo] How to overlay world map using non-standard projection

[Roger Bivand]

On Thu, 20 Jan 2011, Robert Hijmans wrote:

>
>> Problem is that I've been unable to overlay a map of the world. I've tried
> to use the map() function, and >I've generated the map I want (right
> projection and centered correctly) with the call:
>
>> library(maps)
>> world <-map(projection = "cylequalarea", parameter=30, interior=FALSE,
> wrap=TRUE, ylim=c(->60,60),xlim=c(-180,180), orientation=c(90,-205,0))
>
>> But when I try to overlay this map on the plot of my transformed
> SpatialGrid I find that the scale of both >maps is completely different.
> Their coordinates don't match:
>
> Perhaps you need to provide more parameters to map() (ellipsoid, datum?), as
> it appears that it now uses something equivalent to an ellipsoid with radii
> of about ~1 meter

When the map() function in the maps package calls functions in the mapproj 
package to project geographical coordinates to planar, the output 
coordinates are planar but in a nonstandard metric not documented in that 
package. If you would like control over the coordinates, convert the 
object returned by map() with map2SpatialLines() or map2SpatialPolygons() 
in maptools, and project with spTransform() methods in rgdal instead:

library(maps)
wrld <- map(plot=FALSE, interior=FALSE, wrap=TRUE,
   ylim=c(-60,60),xlim=c(-179.9, 179.9))
# shaving off edges to avoid horizontal line artefact
library(maptools)
wrld_sp <- map2SpatialLines(wrld)
proj4string(wrld_sp) <- CRS("+proj=longlat")
library(rgdal)
cea_wrld_sp <- spTransform(wrld_sp, CRS("+proj=cea +lat_ts=30"))
plot(cea_wrld_sp, axes=TRUE)

I'm assuming that you meant CRS("+proj=cea +lat_ts=30"), where the second 
value is the standard parallel.

Roger

>
> Robert
>

-- 
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