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

30

Я накладываю карту мира из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()

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

enter image description here

Совершенно очевидно, что между многоугольниками, которые находятся на одном конце или на другом меридиане, есть линии. Мое текущее решение состоит в том, чтобы заменить точки, близкие к первому меридиану, на 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()

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

There must be a better way of centering the map on another meridian. I tried using the orientation parameter in map, but setting this to orientation = c(0,180,0) did not yield the correct result, in fact it did not change anything to the result object (all.equal yielded TRUE). Getting rid of the horizontal stripes should be possible without deleting some of the polygons. It might be that solving point 1. also solves this point.
Sp классы к чему-то, что может использовать ggplot2, довольно легко Так что я открыт для ответов на основе sp. Paul Hiemstra
Если вы просто заинтересованы в центрировании карты вокруг 180 градусов, карта "world2" вmaps пакет или версия с высоким разрешением «world2Hires» вmapdata пакет уже центрирован на 180 град. Остальная часть вашего кода работает нормально. Jim M.
Вас интересует только решение ggplot? Я бы сделал большую часть этого, используяsp и связанные пакеты, но ничего не знаю о конвертацииsp& APOS; sSpatial* объекты (особенно представляющие растры) для ggplot ... Josh O'Brien
Спасибо за комментарий! Я все еще заинтересован в гибком решении, хотя. Paul Hiemstra
@PaulHiemstra - Это (из месяца или двух предыдущих) также решение?stackoverflow.com/questions/5353184/… Josh O'Brien

Ваш Ответ

3   ответа
1

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

 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.

Еще одно неудобство заключается в том, что в этой точке широта / долгота выражена в рад, а не в градусах:

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()
Error: User Rate Limit Exceededggplot2. Paul Hiemstra
Error: User Rate Limit Exceeded
Error: User Rate Limit Exceeded
27

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

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

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

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

enter image description here

UPDATED

Здесь я положил некоторые объяснения:

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

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

enter image description here

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

И вgeom_path, данные той же группы связаны. Так что еслиmp2$group <- mp2$group + max(mp2$group) + 1 отсутствует, выглядит так: enter image description here

Error: User Rate Limit Exceeded
Error: User Rate Limit ExceededmapError: User Rate Limit Exceededmap. Paul Hiemstra
Error: User Rate Limit Exceededstackoverflow.com/questions/47021117/…
Error: User Rate Limit Exceededmap(orientation = c(90, 0, 180), projection="mercator")Error: User Rate Limit Exceeded
Error: User Rate Limit Exceeded
17

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

  1. Converting the world map from the maps package into a SpatialLines object with a geographical (lat-long) CRS.
  2. Projecting the SpatialLines map into the Plate Carée (aka Equidistant Cylindrical) projection centered on the Prime Meridian. (This projection is very similar to a geographical mapping).
  3. Cutting in two segments that would otherwise be clipped by left and right edges of the map. (This is done using topological functions from the rgeos package.)
  4. Reprojecting to a Plate Carée projection centered on the desired meridian (lon_0 in terminology taken from the PROJ_4 program used by spTransform() in the rgdal package).
  5. Identifying (and removing) any remaining 'streaks'. I automated this by searching for lines that cross g.e. two of three widely separated meridians. (This also uses topological functions from the rgeos package.)

Это, очевидно, большая работа, но она оставляет карты с минимально усеченными картами, которые можно легко перепроецировать с помощью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'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

Error: User Rate Limit ExceededrgeosError: User Rate Limit ExceededgBuffer()Error: User Rate Limit ExceededgDifference()Error: User Rate Limit ExceededrgeosError: User Rate Limit ExceededmapError: User Rate Limit Exceeded
Error: User Rate Limit Exceeded
Error: User Rate Limit Exceeded

Похожие вопросы