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

查看：221
本文介绍了在绘制世界地图时使用与主要子午线不同的中心的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！

问题描述

我将 maps 包中的世界地图覆盖到 ggplot2 栅格几何上。然而，这个栅格不是以主子午线（0度）为中心，而是在180度（大约白令海和太平洋）。以下代码获取地图并重新映射180度地图： 

  require（地图）
 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代替接近主子午线的点，用中的调用替换：

<$ p （世界地图，{
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 产生 TRUE ）。

  2. 在不删除某些多边形的情况下，应该可以摆脱水平条纹。这可能是解决方法1。也解决了这一问题。 。它的工作原理是：


    1. 将世界地图从 maps 包转换为 SpatialLines 具有地理（lat-long）CRS的对象。

    2. 投影 SpatialLines 映射到以素数子午线为中心的PlateCarée（aka等距圆柱形）投影。 （此投影与地理映射非常相似）。
    3. 剪切两段，否则将被地图的左右边缘剪切。 （这是使用来自 rgeos 包的拓扑函数完成的。）

    4. 重新投影到PlateCarée投影，集中于所需子午线（<从 PROJ_4 程序中使用的术语中，code> lon_0 spTransform（）在 rgdal 包中）。

    5. 识别（并删除）任何剩余的条纹。我通过搜索通过g.e.的线来自动执行此操作。三个分开的经络中的两个。 （这也使用来自 rgeos 包的拓扑函数。）

    这是显然有很多工作要做，但留下一个最小截断的地图，并且可以使用 spTransform（）轻松重新投影。要使用 base 或 lattice 图形将这些覆盖在光栅图像上，我首先重新映射栅格，同时使用 spTransform（）。如果你需要它们，网格线和标签也可以预计为匹配 SpatialLines 地图。





     库（sp）
    库（映射）
    库（maptools）## map2SpatialLines（），pruneMap（）
    库（rgdal）## CRS（），spTransform（）$ b $库库（rgeos）## readWKT（），gIntersects（），gBuffer（），gDifference（）
     
     ## map映射到SpatialLines映射
     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整齐地沿着反对子午线
     sliceAtAntipodes <  - 函数（SLmap，lon_0）{
     ##预备
     long_180<  - （lon_0 %% 360） -  180 
     llCRS<  -  CRS（+ proj = longlat + ellps = 地图对象的## CRS 
     eqcCRS < -  CRS（+ proj = eqc）
     ##将地图重新​​映射到Equidistant圆柱/板块Caree投影
     SLmap< ;  - 线（cbind（long_180，c（-89,89））），ID（线）（cbind（long_180，c（-89,89）））;  - 在子午线上形成一个窄的SpatialPolygon =cutter）
     SL < -  SpatialLines（list（L），proj4string = llCRS）
     SP < -  gBuffer（spTransform（SL，eqcCRS），10，byid = TRUE）
     ##用它来剪切任何它横过
     ii<  -  which（gIntersects（SLmap，SP，byid = TRUE））的SpatialLines段
    ＃用分割版本替换违规行
    ＃ （但是当没有交叉点时跳过（例如，当lon_0 = 0时））
     if（length（ii））{
     SPii < -  g Difference（SLmap [ii]，SP，byid = TRUE）
     SLmap<  -  rbind（SLmap [-ii]，SPii）
    } 
     return（SLmap）
    } 
     
     ## re-中心，并清理剩余的条纹
     recenterAn dClean < - 函数（SLmap，lon_0）{
     llCRS <-CRS（+ proj = longlat + ellps = WGS84）## map包的CRS 
     newCRS < -  CRS（paste（ + proj = eqc + lon_0 =，lon_0，sep =））
     ##接收者
     SLmap<  -  spTransform（SLmap，newCRS）
     ##识别剩余的'scratch-lines '通过搜索
     ##经过3条经线的2条线，间隔120度
     v1 <-spTransform（readWKT（LINESTRING（-62 -89，-62 89）） （读取WKT（LINESTRING（58-89,58 89），p4s = 11CRS），newCRS）
     v3< -spTransform（readWKT LINESTRING（178 -89，178 89），p4s = 11CRS），newCRS）
     ii<  - （（gIntersects（v1，SLmap，byid = TRUE）+ 
     gIntersects（v2，SLmap ，byid = TRUE）+ 
     gIntersects（v3，SLmap，byid = TRUE））> = 2）
     SLmap [-ii] 
    } 
     
     ##把它放在一起：
     Recenter<  -  function（lon_0 = -100，grid = FALSE，...）{ 
     SLmap<  -  makeSLmap（）
     SLmap2<  -  sliceAtAntipodes（SLmap，lon_0）
     recenterAndClean（SLmap2，lon_0）
    } 
     
     ##试试看
     par（mfrow = c（2,2），mar = rep（1,4））
     plot（Recenter（-90），col =grey40）; box（）##以90w 
    图为中心（Recenter（0），col =grey40）; box（）##以素数子午线
    图为中心（Recenter（90），col =grey40）; box（）##以90e 
    图为中心（Recenter（180），col =grey40）; box（）##以国际日期行为中心


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

    which yields the following output:

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

    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屋！

查看全文
相关文章
其他开发最新文章
热门教程
热门工具
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆