在绘制世界地图时使用与本初子午线不同的中心 [英] Use different center than the prime meridian in plotting a world map

问题描述

我正在将 maps 包中的世界地图叠加到 ggplot2 光栅几何图形上.但是，此栅格不是以本初子午线(0 度)为中心，而是以 180 度(大致是白令海和太平洋)为中心.以下代码获取地图并将地图重心 180 度:

I am overlaying a world map from the maps package onto a ggplot2 raster geometry. However, this raster is not centered on the prime meridian (0 deg), but on 180 deg (roughly the Bering Sea and the Pacific). The following code gets the map and recenters the map on 180 degree:

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 调用替换为:

Quite obviously there are the lines draw between polygons that are on one end or the other of the prime meridian. My current solution is to replace points close to the prime meridian by NA, replacing the within call above by:

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

这会导致正确的图像.我现在有很多问题:

Which leads to the correct image. I now have a number of question:

1. 必须有更好的方法将地图居中放置在另一条子午线上.我尝试在 map 中使用 orientation 参数，但是将其设置为 orientation = c(0,180,0) 并没有产生正确的结果，在事实上，它没有对结果对象进行任何更改(all.equal 产生了 TRUE).
2. 应该可以在不删除某些多边形的情况下去除水平条纹.可能是解决第 1 点.也解决了这一点.

1. 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).
2. 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.

推荐答案

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

Here's a different approach. It works by:

1. 将世界地图从 maps 包转换为具有地理(经纬度)CRS 的 SpatialLines 对象.
2. 将 SpatialLines 地图投影到以本初子午线为中心的 Plate Carée(又名等距圆柱)投影中.(此投影与地理地图非常相似).
3. 切割成两段，否则会被地图的左右边缘剪裁.(这是使用 rgeos 包中的拓扑函数完成的.)
4. 重新投影到以所需子午线为中心的 Plate Carée 投影(术语中的 lon_0 取自 spTransform() 中使用的 PROJ_4 程序rgdal 包).
5. 识别(并移除)任何剩余的条纹".我通过搜索跨越 g.e.三个相距很远的经络中的两条.(这也使用了 rgeos 包中的拓扑函数.)

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 地图.

This is obviously a lot of work, but leaves one with maps that are minimally truncated, and can be easily reprojected using spTransform(). To overlay these on top of raster images with base or lattice graphics, I first reproject the rasters, also using spTransform(). If you need them, grid lines and labels can likewise be projected to match the SpatialLines map.

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

这篇关于在绘制世界地图时使用与本初子午线不同的中心的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！