[R-sig-Geo] seamless maps

[Roger Bivand]

On Thu, 8 Jan 2009, Matt Oliver wrote:

> Hans,
>
> Take a look at the proj4 package, specifically project()

Or rather rgdal, which is a fully featured implementation with a proper 
interface to the PROJ.4 library. Unfortunately, there does not seem to be 
a way of shifting the prime meridian in a helpful way - and the original 
question was not about projection, but about joining +180 to -180, which 
is obvious in spherical but not in planar terms. rgdal provides methods 
for projection and datum transformation for objects defined in the sp 
package, as well as for projection of matrices of coordinates.

Roger

PS:

This gives something, but at the cost of meaningless coordinates and 
artefacts for Greenland and Iceland, which could be subsetted out:

wrld_simpl_nA <- wrld_simpl[wrld_simpl$NAME != "Antarctica",]
library(rgdal)
wrld_simpl_nA_22k <- spTransform(wrld_simpl_nA,
   CRS("+proj=merc +datum=WGS84 +x_0=22000000"))
proj4string(wrld_simpl_nA_22k) <- CRS("+proj=merc +datum=WGS84")
wrld_simpl_nA_22k_ll <- spTransform(wrld_simpl_nA_22k,
   CRS("+proj=longlat"))
plot(wrld_simpl_nA_22k_ll, axes=TRUE)

The alternative is as given before - dividing the polygons into ones to be 
recentered and ones to be left as is.

>
> see http://trac.osgeo.org/proj/
>
> There are many projections that allow you to recenter a map on a specific
> longitude
>
> For example
>
>
> x <- seq(-180, 180)
> y <- seq(-90, 90)
>
> xy <- expand.grid(x, y)
>
> plot(xy)
>
> mer <- cbind(0, y)  #prime meridian
> points(mer[,1], mer[,2], col = "white")
>
>
> proj.xy <- project(xy, "+proj=eck4") ###project grid into eckert IV
> proj.mer <- project(mer, "+proj=eck4") ###project meridian into eckert IV
>
> plot(proj.xy)
> points(proj.mer[,1], proj.mer[,2], col = "white")
>
>
> proj.xy <- project(xy, "+proj=eck4 +lon_0=90W") ###project grid into eckert
> IV with central meridian at 90W
> proj.mer <- project(mer, "+proj=eck4 +lon_0=90W") ###project meridian into
> eckert IV with central meridian at 90W
>
>
> plot(proj.xy)
> points(proj.mer[,1], proj.mer[,2], col = "white")
>
>
>
> On Thu, Jan 8, 2009 at 4:07 AM, Hans-J?rg Bibiko <bibiko at eva.mpg.de> wrote:
>
>> Hi,
>>
>> I wonder if there's a workaround to generate seamless maps based on
>> SpatialPolygons objects. E.g. to produce a map showing Australia at the left
>> edge and America at the right one ( xlim := min: 110?E , max: 30?W ).
>> Or the way around, is there a function to re-center a given map by defining
>> the median longitude? (I know the function 'recenter' which produce a map
>> from 0? to 360?)
>>
>> Kind regards,
>>
>> --Hans
>> _______________________________________________
>> 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