Skip to content

seamless maps

8 messages · Matt Oliver, Roger Bivand, Hans-Jörg Bibiko

Original
Hans-Jörg Bibiko
# 
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
Roger Bivand
#
On Thu, 8 Jan 2009, Hans-J?rg Bibiko wrote:

            

      
      
      
    

        
Not without recentering, I think. As you've probably found out:

library(maptools)
load(url("http://spatial.nhh.no/R/etc/TM_WORLD_BORDERS_SIMPL-0.2.RData"))
rc <- recenter(world_simp)

gives nasty artefacts where polygons cross the Prime Meridian, so subsets 
of countries crossing 180? and west of but not crossing the Prime Meridian 
need to be recentered, and added back into the object. That is a good deal 
of work. Do you need the country polygons, or can you manage with just the 
shoreline (as polygons):

gshhs.c.b <- system.file("share/gshhs_c.b", package="maptools")
Px <- c(110, 330)
Py <- c(-62, 90)
P <- Rgshhs(gshhs.c.b, xlim=Px, ylim=Py, level=1)
plot(P$SP)

A coarse borders file is also available, but only as lines, and varies in 
levels of detail between countries see ?Rgshhs.

If you do need the country polygons, we'd need to list the blocks 
of polygons needing recentering. The maps package also has Pacific view 
databases for shorelines and some borders.

Hope this helps,

Roger

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Hans-Jörg Bibiko
    

    

      
      #
    

  


  

      
On 08.01.2009, at 14:27, Roger Bivand wrote:

        

            

      
      
      
    

        
Dear Roger and Matt,

many thanks for the hints. I'll have a look at these.


Best,

--Hans

**********************************************************
Hans-Joerg Bibiko
Max Planck Institute for Evolutionary Anthropology
Department of Linguistics
Deutscher Platz 6     phone:   +49 (0) 341 3550 341
D-04103 Leipzig       fax:     +49 (0) 341 3550 333
Germany               e-mail:  bibiko at eva.mpg.de

  
  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Thu, 8 Jan 2009, Matt Oliver wrote:

        

            

      
      
      
    

        
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.

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Thu, 8 Jan 2009, Hans-J?rg Bibiko wrote:

        

            

      
      
      
    

        
A further idea if you don't need the polygons is:

wrld_simpl_SL <- as(as(wrld_simpl, "SpatialPolygons"), "SpatialLines")
wrld_simpl_SL_0_360 <- recenter(nowrapSpatialLines(wrld_simpl_SL))
plot(wrld_simpl_SL_0_360, axes=TRUE, xlim=c(110, 330), ylim=c(-80, 80))

which keeps the correct coordinates but at the cost of losing the polygons 
and a couple of artefact lines at 180?.

Roger

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Hans-Jörg Bibiko
    

    

      
      #
    

  


  

      
On 08.01.2009, at 15:05, Roger Bivand wrote:

        

            

      
      
      
    

        
Thanks! I'm just playing with it.
Actually I do need the polygons as well.

A VERY quick&dirty&crude solution would be to pack the wrld_simpl  
polygons three times in a new SpatialPolygons object:

-the original one from -180 to 180
-a transformed one from 180 to 540
-a transformed one from -540 to -180

then by specifying xlim only one should be able to center the map for  
any longitude (to get only that clipping).
But I know it's a crude approach and then if one wants to put  
additional points/lines/etc. on it one has to treble them as well.

But anyway, thanks so much for the hints. If I found a way I'll post it.

--Hans

**********************************************************
Hans-Joerg Bibiko
Max Planck Institute for Evolutionary Anthropology
Department of Linguistics
Deutscher Platz 6     phone:   +49 (0) 341 3550 341
D-04103 Leipzig       fax:     +49 (0) 341 3550 333
Germany               e-mail:  bibiko at eva.mpg.de

  
  


  


      

    

    
22 days later

    

      
      

        


  

        

      Hans-Jörg Bibiko
    

    

      
      #
    

  


  

      
Hi,

one simple possibility to draw seamless maps based on e.g.  
TM_WORLD_BORDERS_SIMPL-0.2 is to draw each polygon in question twice  
meaning first the 'normal' polygon and afterwards the transformed  
polygon by simply adding 360 to each x-value (for maps with lon from  
-180 to 540).

The following frugal example illustrates it:
< 
------------------------------------------------------------------------>
library("maptools")
load(url("http://spatial.nhh.no/R/etc/ 
TM_WORLD_BORDERS_SIMPL-0.2.RData"))

wrld_simpl <- wrld_simpl[wrld_simpl$NAME != "Antarctica",]

# Aruba is only a dummy to init the plot plus axes
plot(wrld_simpl[wrld_simpl$NAME=="Aruba",],
   xlim=c(130, 240), # which viewport
   ylim=c(-20, 90),  # should be displayed
   asp=1,
   axes=T, pbg='white')


nirvana <- sapply(wrld_simpl at polygons,
   function(x) sapply(x at Polygons,
     function(y) {
       # draw the normal polygon
       polygon(y at coords,
         col=colors()[which(wrld_simpl$ISO3==x at ID)*2+1],
         border=NA
       );
       # draw the transformed polygon (lon + 360deg)
       polygon(y at coords + rep(c(360, 0), each=dim(y at coords)[1]),
        col=colors()[which(wrld_simpl$ISO3==x at ID)*2+1],
        border=NA
       )
     }
   ))

#polygon's col setting is only for demonstration!
< 
------------------------------------------------------------------------>

or try to set the viewport to:
xlim=c(100, 435)
ylim=c(-90, 90)


etc.

--Hans