32

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のパラメータを使用してみましmapたが、これをに設定してorientation = c(0,180,0)も正しい結果が得られませんでした。実際、結果オブジェクトには何も変更されませんでした(yield all.equalTRUE
  2. 一部のポリゴンを削除しなくても、横縞を取り除くことができるはずです。ポイント1を解くこともこのポイントを解くということかもしれません。
4

3 に答える 3

31

これはやや難しいかもしれませんが、次の方法で実行できます。

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

ここに画像の説明を入力してください

ただし、によって、経度0〜360の領域のサブセットをトリミングscale_x_continuous(limits = c(0, 360))できます。

そして、geom_pathでは、同じグループのデータが接続されています。したがって、mp2$group <- mp2$group + max(mp2$group) + 1が存在しない場合は、次のようになります。 ここに画像の説明を入力してください

于 2012-05-17T02:52:10.367 に答える
17

これは別のアプローチです。それはによって動作します:

  1. mapsパッケージからの世界地図をSpatialLines地理的(緯度経度)CRSを持つオブジェクトに変換します。
  2. SpatialLines本初子午線を中心としたプレートカレ(別名正距円筒図法)投影にマップを投影します。(この投影法は地理的マッピングと非常によく似ています)。
  3. マップの左端と右端でクリップされる2つのセグメントにカットします。rgeos(これは、パッケージのトポロジ関数を使用して行われます。)
  4. 目的の子午線を中心とするプレートカレ投影に再投影します(パッケージで使用されているプログラムlon_0から取得した用語で)。PROJ_4spTransform()rgdal
  5. 残っている「ストリーク」を特定(および削除)します。大きく離れた3つの子午線のうち2つを横切る線を検索することで、これを自動化しました。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

ここに画像の説明を入力してください

于 2012-05-25T06:58:13.043 に答える
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は、マップを横切る線をカットします。map(wrap = TRUE)のオプションを介して使用できますが、plot=TRUEの場合にのみ機能します。

残っている厄介な点の1つは、この時点で、緯度/経度が度ではなくラジアンであるということです。

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()
于 2016-01-06T15:14:02.740 に答える