在绘制世界地图时,使用与本初子午线不同的中心

Question

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

这导致了正确的图像.我现在有一些问题:

1. 必须有一种更好的方法将地图定位在另一个子午线上.我尝试使用orientation参数in map,但设置为orientation = c(0,180,0)没有产生正确的结果,实际上它没有改变任何结果对象(all.equalyielding TRUE).
2. 在不删除某些多边形的情况下,应该可以摆脱水平条纹.可能是解决点1也解决了这一点.

Answer 1

这可能有点棘手但你可以做到:

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缺席,它看起来像:

Answer 2

这是一种不同的方法.它的工作原理是:

1. 将世界地图从maps包转换为SpatialLines具有地理(lat-long)CRS的对象.
2. 将SpatialLines地图投影到以Prime Meridian为中心的PlateCatée(aka Equidistant Cylindrical)投影中.(此投影与地理映射非常相似).
3. 切割成两个部分,否则将被地图的左右边缘剪切.(这是使用rgeos包中的拓扑函数完成的.)
4. 重新投影到平板CAREE投影中心的期望子午线(lon_0在术语从所拍摄的PROJ_4使用程序spTransform()在rgdal程序包).
5. 识别(并移除)任何剩余的"条纹".我通过搜索跨越三个分开的经线中的两个的线来自动化.(这也使用了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'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