resumen - paralelos y meridianos

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

Question

Estoy superponiendo un mapa mundial del paquete de maps a una geometría de trama 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:

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 < 1) | (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.

Answer 1

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 programa 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 trama 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 etiquetas 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, e.g., 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

Answer 2

Esto debería funcionar:

 wm <- map.wrap(map(projection="rectangular", parameter=0, orientation=c(90,0,180), plot=FALSE))
world_map <- data.frame(wm[c("x","y")])
names(world_map) <- c("lon","lat")

El map.wrap corta las líneas que cruzan el mapa. Puede usarse a través de una opción para mapear (wrap = TRUE), pero eso solo funciona cuando plot = TRUE.

Una molestia que queda es que en este punto, lat / lon están en rad, no en grados:

world_map$lon <- world_map$lon * 180/pi + 180
world_map$lat <- world_map$lat * 180/pi
ggplot(aes(x = lon, y = lat), data = world_map) + geom_path()

Answer 3

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

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

ACTUALIZADO

Aquí pongo algunas explicaciones:

La información completa se ve así:

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

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: