Líneas en objetos sp reproyectados con proyección Mollweide

Question

Estoy teniendo un problema de impar de tal manera que cuando reproyecto mapas a una proyección Mollweide en R usando sp, obtengo líneas adicionales que se extienden a través del globo. ¿Qué está pasando aquí, y tiene alguna sugerencia para solucionarlo? Aquí hay un ejemplo reproducible:

library(sp)
library(maps)
library(maptools)

worldmap <- map("world", plot=F)
plot(worldmap, type="l)

Hasta ahora, todo va bien. Tengo

Así que ahora reproyecto. Este es un CRS que estoy usando para ser compatible con algunos otros proyectos. Nótese cómo las líneas que abarcan todo el mundo aparecen entonces.

worldmapLines <- map2SpatialLines(worldmap,
                                  proj4string=CRS("+proj=longlat +datum=WGS84"))

worldM <- spTransform(worldmapLines,  
                      CRS("+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs"))

plot(worldM)

Answer 1

Tus datos superan el límite de 180°E en el este de Rusia y en alguna isla del Pacífico. Si quieres que el mapa se vea bien en grados, tienes que cortar tus datos de origen a 179,9°E/W.

Ver mi respuesta aquí (aunque se trata de una visión centrada en el Pacífico)

QGIS muestra los archivos de forma de los países del mundo centrados en el océano Pacífico utilizando Robinson, Miller Cylindrical u otra proyección

Answer 2

Así que, basándome en la respuesta de AndreJ más arriba, preparé esto.

cleanSpLinesForProjection <- function(w){
  #We'll have to go through this one, line by line...
  slot(w, "lines") <- lapply(slot(worldmapLines, "lines"), function(x) {
    coords <- slot(x, "Lines")[[1]]@coords

    #get the rows with too large longitude values
    rIDX <- which(coords >= 180, arr.ind=T)[,1]

    #if there are some, replace them
    if(length(rIDX>0)) coords <- coords[-rIDX,]

    #replace the slot
    slot(x, "Lines")[[1]]@coords <- coords
    x
  })

  w
}

Esto funciona.

worldM <- spTransform(cleanSpLinesForProjection(worldmapLines),  
                      CRS("+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs"))

plot(worldM)

Me gustaría que hubiera un argumento para spTransform o una función genérica para todos los objetos sp. Ah, bueno. He aprendido algo sobre sp aquí y he creado un código que seguramente utilizaré en el futuro. Y ahora tengo que aplicarlo a un objeto SpatialPolygons... ¿Qué podría salir mal?

Answer 3

Sus datos no son válidos:

> bbox(worldmapLines)
         min       max
x -179.95721 190.29080
y  -85.44308  83.57391

y el gráfico es correcto, dado que lo proporcionaste con datos erróneos.

Con la ayuda de Roger Bivand: una versión corregida de este conjunto de datos está disponible en maptools inténtalo:

library(maptools)
data(wrld_simpl)
plot(spTransform(wrld_simpl, CRS=("+proj=moll +lon_0=0 +x_0=0 +y_0=0")))

Roger mencionó además: No hay ninguna manera fácil de romper las geometrías en los límites del "otro lado del mundo". wrld_simpl se basaba en parte en recenter() pero necesitaba intervenciones manuales (esto era con GPCLib, antes de rgeos).

Answer 4

Puede cortar el área más allá de los 180 grados y cambiar a su representación habitual (no hay nada realmente malo en tener grados más allá de 180, pero no es lo esperado por la mayoría del software).

Lo que tienes

library(sp)
library(maps)
library(maptools)
worldmap <- map("world", plot=F)
worldmapLines <- map2SpatialLines(worldmap, proj4string=CRS("+proj=longlat +datum=WGS84"))

Ahora puedes hacerlo:

library(raster)
w1 <- crop( worldmapLines, extent(-180, 180,-90,90))
w2 <- crop( worldmapLines, extent(180, 200,-90,90))
w2 <- shift(w2, -360)
w <- bind(w1, w2)
x <- spTransform(w, CRS("+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs"))
plot(x)