在绘制世界地图时使用与主要子午线不同的中心

我将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代替接近主子午线的点，用下面的withinreplace上面的内容：

 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. 在另一条子午线上定位地图肯定有更好的方法。 我尝试在map使用orientation参数，但是将其设置为orientation = c(0,180,0)并没有产生正确的结果，事实上它没有改变任何东西到结果对象（ all.equal yielded TRUE ）。
2. 摆脱横条纹应该是可能的，而不删除一些多边形。 可能是解决第一点也解决了这一点。

Answer 1

这是一个不同的方法。 它的工作原理是：

1. 将maps包中的世界地图转换为具有地理（纬度经度）CRS的SpatialLines对象。
2. 将SpatialLines地图投影到以素数子午线为中心的PlateCarée（aka等距圆柱形）投影中。 （这个预测与地理映射非常相似）。
3. 切割两段，否则将被地图的左侧和右侧边缘剪切。 （这是使用rgeos包中的拓扑函数完成的。）
4. 重新投影到PlateCarée投影，集中在所需的子午线（从rgdal包中的rgdal spTransform()使用的PROJ_4程序lon_0的术语中， rgdal ）。
5. 识别（并删除）任何剩余的“条纹”。 我通过search三条分开的经线中的两条线进行自动化。 （这也使用rgeos包的拓扑函数。）

这显然是很多工作，但留下一个最小截断的地图，并可以使用spTransform()轻松重新spTransform() 。 为了将这些覆盖在具有base或latticegraphics的光栅图像之上，我首先重新指定光栅，同样使用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, eg, 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

这可能有点棘手，但你可以通过：

 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 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时才起作用。

剩下的一个烦恼是，在这一点上，纬度/经度是弧度，而不是度数：

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