Используйте другой центр, чем основной меридиан, при построении карты мира

Я накладываю карту мира из пакета maps на геометрию растра ggplot2. Однако этот растр не сосредоточен на основном меридиане (0 град), а на 180 град (примерно Беринговом и Тихоокеанском). Следующий код получает карту и недавнюю карту на 180 градусов:

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

который дает следующий результат:

Совершенно очевидно, что между многоугольниками, которые находятся на одном конце или другом простого меридиана, есть линии. Мое текущее решение - заменить точки, близкие к первому меридиану, на NA, заменив вызов within выше:

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

Это приводит к правильному изображению. У меня теперь есть ряд вопросов:

Должен быть лучший способ центрирования карты на другом меридиане. Я попытался использовать параметр orientation в map, но установка этого параметра в orientation = c(0,180,0) не дала правильного результата, фактически ничего не изменила для объекта результата (all.equal дано TRUE).
Избавление от горизонтальных полос должно быть возможным без удаления некоторых полигонов. Возможно, решающий пункт 1. также решает эту проблему.

Ответ 1

Здесь другой подход. Он работает:

Преобразование карты мира из пакета maps в объект SpatialLines с географическим (lat-long) CRS.
Проецирование карты SpatialLines в проекцию Plate Carée (аналогично эквидистантной цилиндрической) с центром на Prime Meridian. (Эта проекция очень похожа на географическое отображение).
Резка в двух сегментах, которые в противном случае были бы обрезаны левым и правым краями карты. (Это делается с использованием топологических функций из пакета rgeos.)
Репроектируя проекцию Plate Carée, центрированную на желаемый меридиан (lon_0 в терминологии, взятой из программы PROJ_4, используемой spTransform() в пакете rgdal).
Идентификация (и удаление) оставшихся "полос". Я автоматизировал это, ища линии, пересекающие g.e. два из трех широко разделенных меридианов. (Это также использует топологические функции из пакета rgeos.)

Это, очевидно, много работы, но оставляет ее с минимально усеченными картами и может быть легко перепрограммирована с помощью spTransform(). Чтобы наложить их поверх растровых изображений с графикой base или lattice, я сначала перепрограммирую растры, также используя spTransform(). Если они вам понадобятся, линии сетки и метки могут быть спроектированы так, чтобы они соответствовали карте 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 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

Ответ 2

Это может быть несколько сложно, но вы можете сделать это:

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

С помощью этой настройки вы можете легко установить центр (т.е. пределы):

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

ОБНОВЛЕНО

Вот несколько пояснений:

Все данные выглядят так:

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

но scale_x_continuous(limits = c(0, 360)) вы можете обрезать подмножество области от 0 до 360 долготы.

И в geom_path подключены данные той же группы. Поэтому, если mp2$group <- mp2$group + max(mp2$group) + 1 отсутствует, это выглядит так:

Ответ 3

Это должно работать:

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

Map.wrap разрезает линии, проходящие через карту. Его можно использовать с помощью опции для сопоставления (wrap = TRUE), но это работает только при plot = TRUE.

Одно из оставшихся досаду заключается в том, что в этот момент lat/lon находятся в раде, а не в градусах:

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