Use un centro diferente al meridiano principal para trazar un mapa mundial

Estoy superponiendo un mapa mundial del paquete de maps a una geometría de ttwig ggplot2 . Sin embargo, este ráster no está centrado en el meridiano principal (0 grados), sino en 180 grados (aproximadamente el Mar de Bering y el Pacífico). El siguiente código obtiene el mapa y vuelve a centrar el mapa en 180 grados:

 require(maps) world_map = data.frame(map(plot=FALSE)[c("x","y")]) names(world_map) = c("lon","lat") world_map = within(world_map, { lon = ifelse(lon < 0, lon + 360, lon) }) ggplot(aes(x = lon, y = lat), data = world_map) + geom_path() 

que produce el siguiente resultado:

enter image description here

Obviamente, hay líneas dibujadas entre polígonos que están en un extremo u otro del meridiano principal. Mi solución actual es reemplazar los puntos cercanos al meridiano principal por NA, reemplazando el llamado within anterior por:

 world_map = within(world_map, { lon = ifelse(lon < 0, lon + 360, lon) lon = ifelse((lon  359), NA, lon) }) ggplot(aes(x = lon, y = lat), data = world_map) + geom_path() 

Lo que lleva a la imagen correcta. Ahora tengo una cantidad de preguntas:

  1. Debe haber una mejor manera de centrar el mapa en otro meridiano. Intenté usar el parámetro de orientation en el map , pero al establecer esto en orientation = c(0,180,0) no orientation = c(0,180,0) el resultado correcto, de hecho no cambió nada al objeto resultante ( all.equal arrojó TRUE ).
  2. Deshacerse de las rayas horizontales debería ser posible sin eliminar algunos de los polígonos. Puede ser que resolver el punto 1. también resuelva este punto.

Aquí hay un enfoque diferente. Funciona por:

  1. Conversión del mapa mundial del paquete de maps en un objeto SpatialLines con un CRS geográfico (lat-long).
  2. Proyectando el mapa SpatialLines en la proyección Plate Carée (también conocida como Equidistant Cylindrical) centrada en el Meridiano de Greenwich. (Esta proyección es muy similar a un mapeo geográfico).
  3. Cortar en dos segmentos que, de lo contrario, quedarían recortados por los bordes izquierdo y derecho del mapa. (Esto se hace usando funciones topológicas del paquete rgeos ).
  4. Reproyectando a una proyección Plate Carée centrada en el meridiano deseado ( lon_0 en la terminología tomada del progtwig PROJ_4 utilizado por spTransform() en el paquete rgdal ).
  5. Identificando (y eliminando) cualquier ‘raya’ restante. Automaticé esto buscando las líneas que cruzan dos de los tres meridianos ampliamente separados. (Esto también usa funciones topológicas del paquete rgeos ).

Obviamente, esto es mucho trabajo, pero deja uno con mapas que están mínimamente truncados, y se pueden reproyectar fácilmente usando spTransform() . Para superponer estos en la parte superior de las imágenes de ttwig con gráficos de base o en lattice , primero reproyecto los rásteres, también utilizando spTransform() . Si los necesita, también se pueden proyectar líneas de cuadrícula y tags para que coincidan con el mapa de SpatialLines .


 library(sp) library(maps) library(maptools) ## map2SpatialLines(), pruneMap() library(rgdal) ## CRS(), spTransform() library(rgeos) ## readWKT(), gIntersects(), gBuffer(), gDifference() ## Convert a "maps" map to a "SpatialLines" map makeSLmap <- function() { llCRS <- CRS("+proj=longlat +ellps=WGS84") wrld <- map("world", interior = FALSE, plot=FALSE, xlim = c(-179, 179), ylim = c(-89, 89)) wrld_p <- pruneMap(wrld, xlim = c(-179, 179)) map2SpatialLines(wrld_p, proj4string = llCRS) } ## Clip SpatialLines neatly along the antipodal meridian sliceAtAntipodes <- function(SLmap, lon_0) { ## Preliminaries long_180 <- (lon_0 %% 360) - 180 llCRS <- CRS("+proj=longlat +ellps=WGS84") ## CRS of 'maps' objects eqcCRS <- CRS("+proj=eqc") ## Reproject the map into Equidistant Cylindrical/Plate Caree projection SLmap <- spTransform(SLmap, eqcCRS) ## Make a narrow SpatialPolygon along the meridian opposite lon_0 L <- Lines(Line(cbind(long_180, c(-89, 89))), ID="cutter") SL <- SpatialLines(list(L), proj4string = llCRS) SP <- gBuffer(spTransform(SL, eqcCRS), 10, byid = TRUE) ## Use it to clip any SpatialLines segments that it crosses ii <- which(gIntersects(SLmap, SP, byid=TRUE)) # Replace offending lines with split versions # (but skip when there are no intersections (as, eg, when lon_0 = 0)) if(length(ii)) { SPii <- gDifference(SLmap[ii], SP, byid=TRUE) SLmap <- rbind(SLmap[-ii], SPii) } return(SLmap) } ## re-center, and clean up remaining streaks recenterAndClean <- function(SLmap, lon_0) { llCRS <- CRS("+proj=longlat +ellps=WGS84") ## map package's CRS newCRS <- CRS(paste("+proj=eqc +lon_0=", lon_0, sep="")) ## Recenter SLmap <- spTransform(SLmap, newCRS) ## identify remaining 'scratch-lines' by searching for lines that ## cross 2 of 3 lines of longitude, spaced 120 degrees apart v1 <-spTransform(readWKT("LINESTRING(-62 -89, -62 89)", p4s=llCRS), newCRS) v2 <-spTransform(readWKT("LINESTRING(58 -89, 58 89)", p4s=llCRS), newCRS) v3 <-spTransform(readWKT("LINESTRING(178 -89, 178 89)", p4s=llCRS), newCRS) ii <- which((gIntersects(v1, SLmap, byid=TRUE) + gIntersects(v2, SLmap, byid=TRUE) + gIntersects(v3, SLmap, byid=TRUE)) >= 2) SLmap[-ii] } ## Put it all together: Recenter <- function(lon_0 = -100, grid=FALSE, ...) { SLmap <- makeSLmap() SLmap2 <- sliceAtAntipodes(SLmap, lon_0) recenterAndClean(SLmap2, lon_0) } ## Try it out par(mfrow=c(2,2), mar=rep(1, 4)) plot(Recenter(-90), col="grey40"); box() ## Centered on 90w plot(Recenter(0), col="grey40"); box() ## Centered on prime meridian plot(Recenter(90), col="grey40"); box() ## Centered on 90e plot(Recenter(180), col="grey40"); box() ## Centered on International Date Line 

enter image description here

Esto puede ser algo complicado, pero puedes hacerlo de la siguiente manera:

 mp1 <- fortify(map(fill=TRUE, plot=FALSE)) mp2 <- mp1 mp2$long <- mp2$long + 360 mp2$group <- mp2$group + max(mp2$group) + 1 mp <- rbind(mp1, mp2) ggplot(aes(x = long, y = lat, group = group), data = mp) + geom_path() + scale_x_continuous(limits = c(0, 360)) 

enter image description here

Con esta configuración, puede establecer fácilmente el centro (es decir, límites):

 ggplot(aes(x = long, y = lat, group = group), data = mp) + geom_path() + scale_x_continuous(limits = c(-100, 260)) 

enter image description here

ACTUALIZADO

Aquí pongo algunas explicaciones:

La información completa se ve así:

 ggplot(aes(x = long, y = lat, group = group), data = mp) + geom_path() 

enter image description here

pero por scale_x_continuous(limits = c(0, 360)) , puede recortar un subconjunto de la región de 0 a 360 de longitud.

Y en geom_path , los datos del mismo grupo están conectados. Entonces, si mp2$group <- mp2$group + max(mp2$group) + 1 está ausente, parece que: enter image description here