Rの geosphere パッケージ(2)：大圏航路を描く

はじめに

この記事では、R言語の geosphere パッケージを使って大圏航路を求め、それを地図上に描いてみたいと思う。

準備

まずは、必要となるパッケージを読みこんでおく。

library("geosphere")
library("tidyverse")

また、東京 (tokyo)、リマ(lima)、ロンドン (london) の3地点の経度・緯度を定義しておく。

tokyo <- c(139.766905, 35.681242)
lima <- c(-77.029796, -12.044866)
london <- c(-0.118092, 51.509865)

大圏航路の計算

2点間の大圏航路を求めるには、gcIntermediate()という函数を使う。例えば、東京とロンドンの間の大圏航路ならば、先ほど定義しておいた地点を用いて、次のように計算することができる。

gc1 <- gcIntermediate(tokyo, london, addStartEnd = TRUE)

ここで作成した gc1 は、大圏航路上の地点の経度と緯度の値を並べた行列になっている。なお、オプションの addStartEnd = TRUE は出発点と到着点を結果の中に含めることを指定している。デフォルトではこのオプションは FALSE とされているので注意が必要である。

大圏航路の図示

geosphere パッケージには、簡単な世界地図のデータが入れられており、これを使って大西洋を中心とする世界地図を描くことができる。

data(wrld)
plot(wrld, type="l")

これに先ほど用意した大圏航路を書き加えるには、単に lines() 函数を使えば良い。

data(wrld)
plot(wrld, type="l")
lines(gc1, col="red")

日付変更線をまたぐ場合

東京・ロンドン間の大圏航路を図示したのと同様に、東京・リマ間の大圏航路を図示してみよう。

gc2 <- gcIntermediate(tokyo, lima, addStartEnd = TRUE)
plot(wrld, type="l")
lines(gc2, col="red")

すると、大圏航路らしいものとは別に、左右にのびる水平線が描かれてしまう。これは、東京・リマ間の大圏航路が経度180度線（ほぼ日付変更線に相当）をまたぐために起こる現象である。要するに図の右端と左端は本当は同じところなのだが、データの都合上、右端から左端に行かなくてはならず、それらを結ぶために余計な水平線が引かれてしまったというわけだ。

このことを防ぐためには、gcIntermediate() で大圏航路を計算するときにオプションとして breakAtDateLine = TRUE を指定すれば良い。すなわち、次のようにする。

addStartEnd = TRUE,
breakAtDateLine = TRUE)
plot(wrld, type="l")
lapply(gc2_break, lines, col="red")

breakAtDateLine = TRUE を指定すると、結果が1つの行列ではなく、2つの行列からなるリストになる。つまり、経度180度線を越える前の部分を示す行列と、越えた後の部分を示す行列に分かれるわけである。このため、大圏航路を図示するときは、lapply() 関数を用いて2つの行列を1つずつ図示するようにしている。

応用1：複数の大圏航路の図示

大圏航路の図示については今まで説明してきたところでほぼ説明しつくしているのだが、参考までにもう少し図示の例を挙げてみたい。具体的には、東京から北京・ソウル・札幌といった都市への大圏航路を図示してみたい。

以下に挙げる例では、世界地図のデータとして rworldmap パッケージに含まれているものを用いる。また、図示の際には ggplot2 パッケージを用いる ((以下では、ggplot2 を直接読みこむのではなく、tidyverse を読みこむことで ggplot2 パッケージに含まれる函数を使えるようにしている。)) 。というわけで、まずは必要なパッケージを読みこんでおこう。

library("geosphere")
library("tidyverse")
library("rworldmap")
library("ggrepel")
library("broom")

次に、rworldmap パッケージに含まれている世界地図のデータを取り出す。そのままではggplot2 パッケージでの作図に用いることができないので、broom パッケージの tidy() という函数を使って扱いやすいデータフレームに変換しておく ((ggplot2 パッケージで作図する際のデータは整然データ (tidy data) という形にしておく必要がある。しかし、R で作られるデータは整然データになっていないものも多い。broom パッケージを使うと、R のさまざまなオブジェクトを整然データに変えることができる。)) 。

# rworldmap パッケージの低解像度の世界地図のデータを取得
wmap <- getMap("low")
# ggplot2 パッケージで扱いやすいデータフレームに変換
wmap_df <- tidy(wmap)

大圏航路の計算に移ろう。まず、北京・ソウル・札幌の3都市の経度・緯度の情報を収めた destinations というデータフレームを作成する。次に、それぞれの地点に対して東京との大圏航路を1つずつ計算し、その結果を gcs というデータフレームに入れる。最後に、経度・緯度の情報が入った destinations に大圏航路の計算結果である gcs を合わせる。

destinations <- data.frame(
  City = c("Beijing", "Seoul", "Sapporo"),
  Long = c(116.397726, 126.991339, 141.353668),
  Lati = c(39.904689, 37.579154, 43.062485)
)
gcs <- destinations %>% rowwise() %>%
  do(GC = gcIntermediate(tokyo,
                         c(.$Long, .$Lati),
                         addStartEnd = TRUE))
destinations <- cbind(destinations, gcs)

あとは、このデータをもとに図示するだけだ。

ggplot() +
  # 世界地図の描画
  geom_polygon(data = wmap_df,
               aes(x = long, y = lat, group = group),
               fill = "lightgreen") +
  # 大圏航路を1つずつ図示
  map(destinations$GC, function(gc){
    geom_path(data = as.data.frame(gc),
              aes(x = lon, y = lat))
  }) +
  # 北京・ソウル・札幌の図示
  geom_point(data = destinations,
             aes(x = Long, y = Lati)) +
  geom_text_repel(data = destinations,
                  aes(x = Long, y = Lati, label = City)) +
  # 東京の図示
  geom_point(aes(x = tokyo[1], y = tokyo[2])) +
  geom_text_repel(aes(x = tokyo[1], y = tokyo[2],
                      label = "Tokyo")) +
  # 東経115度から145度、北緯30度から45度の部分だけ切り出し
  coord_cartesian(xlim=c(115, 145), ylim=c(30, 45)) +
  # 軸ラベル
  xlab("東経") + ylab("北\n緯") +
  theme(axis.title.y = element_text(angle = 0, vjust= 0.5))

結果として以下のような図が得られる。

応用2：フライト数の多さを反映した大圏航路の図示

もう1つ図示の例を挙げよう。今度は、2013年にニューヨーク市のジョン・F・ケネディ空港を出発したフライトの大圏航路を図示する。図示に当たっては、フライト数が多い大圏航路ほど目立つようにしたい。

2013年にニューヨーク市の空港から出発したフライトのデータをまとめた nycflights13 というパッケージがある。今回はこのパッケージのデータを用いる。

まずは、2013年のフライトのデータから行先別フライト数を求め、合わせて行先の空港の位置の情報をまとめておく。

library(nycflights13)
# 空港の位置の情報を取り出す
ap_location <- airports %>% select(faa, lon, lat)
# フライトのデータからジョン・F・ケネディ空港を出発したものを
# 抽出し、行先別フライト数を求め、空港の位置の情報と統合
flights %>% filter(origin == "JFK") %>%
  group_by(dest) %>% summarise(num = n()) %>%
  inner_join(ap_location, by = c("dest" = "faa")) -> flt

まとめた情報をもとに、ジョン・F・ケネディ空港からの大圏航路を求める。

# ジョン・F・ケネディ空港の位置
jfk <- c(-73.77893, 40.63975)
# 大圏航路を求める
gcs2 <- flt %>% rowwise() %>%
  do(GC = gcIntermediate(jfk,
                         c(.$lon, .$lat),
                         addStartEnd = TRUE) %>%
       as.tibble())
# 求めた大圏航路を先に作った行先別フライト数等の情報と統合
flt <- bind_cols(flt, gcs2)

ここまでできればデータの準備はほぼ終わりなのだが、図示するときのためにデータフレームを変形しておく。何をしているかというと、大圏航路の情報を収めた個々のデータフレームの中に、フライト数の情報を入れこんでいる。

flt %>% unnest(.id = "temp_id") %>%
  group_by(temp_id, dest, lon, lat) %>% nest() -> flt

後は描画するだけだ。大圏航路そのものの図示は、geom_path() を使っている。ここで size=num, alpha=num と指定することで、フライト数 (num) が多い航路ほど太く濃く図示されるようになっている。

ggplot() +
  # 世界地図の描画
  geom_polygon(data = wmap_df,
               aes(x = long, y = lat, group = group),
               fill = "lightgreen") +
  # 大圏航路を1つずつ図示
  map(flt$data, function(data){
    geom_path(data = data,
              aes(x = lon1, y = lat1,
                  size=num, alpha=num))
  }) +
  # 東経115度から145度、北緯30度から45度の部分だけ切り出し
  coord_cartesian(xlim=c(-170, -50), ylim=c(20, 50)) +
  # グラフの見た目を適宜調整
  scale_size(range = c(0.1, 1.5)) +
  scale_alpha(range = c(0.5, 1)) +
  theme(legend.position="none") +
  theme(aspect.ratio=2/5)

描画した結果は以下の通りである。フライト数が多い行先ほど太く見えるようになっている。ロサンゼルスとサンフランシスコという西海岸の二大都市に行くフライトが多いだけでなく、フロリダの方に行くフライトも結構多いということが分かる。