Colorless Green Ideas

Rの geosphere パッケージ(1)：2点間の距離と方位角を計算する

概要
Rの geosphere パッケージの導入方法、地球上の2点間の距離・方位角の求め方の説明。

はじめに

これから何回かに分けて、R言語の geosphere パッケージについて説明していきたいと思う。このパッケージは、球面三角法 (spherical trigonometry) に関する計算を行うパッケージであり、地球上の2点間の距離を計算したり、大圏航路を求めたりすることができる。

今回はまずパッケージの導入と、2点間の距離と方位角の計算方法について説明していきたいと思う。

geosphere パッケージの導入

geosphere パッケージはCRANに入っているので、以下のようにして簡単にインストールすることができる。

install.packages("geosphere")

実際に使うときには、library()を使って読みこんでおこう。

library("geosphere")

地点の表現

geosphere パッケージを使った計算に当たっては、地球上の地点を経度 (longitude) と緯度 (latitude) を使って表現する必要がある。具体的には、経度を第1の要素とし、緯度を第2の要素とする長さ2の数値ベクトルで表現する。経度が先に来るので要注意。

例えば、東京は東経139.766905度、北緯35.681242度にあるから、c(139.766905, 35.681242) と表すことができる。

また、西経・南緯については負数を用いて表す。例えば、ペルーの首都リマは西経77.029796度、南緯12.044866にあるから、c(-77.029796, -12.044866) と表現する。

とりあえず、東京 (tokyo)、リマ(lima)、ロンドン (london) の3地点を定義しておこう。

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

2点間の距離の計算

geosphere パッケージには、地球上の2点間の距離を求める函数がいくつか用意されている。パッケージのドキュメンテーションを見る限り、正確な値を出したければ distGeo() という函数を使うのが良さそうである。例えば、東京とロンドンの間の距離を出すには以下のようにする。

distGeo(tokyo, london)

すると、9585300 という値が得られる。計算結果の単位はメートルであり、東京とロンドンの間の距離は9,585,300メートル、すなわち約9,585キロメートルということが分かる。

距離の計算のための函数

geosphere パッケージに含まれている2点間の距離を求める函数としては以下のものがある。いずれも結果はメートルで出力される。

  • 地球を真球と見なして距離を求めるもの
    • distCosine()
    • distVincentySphere()
    • distHaversine()
  • 地球を楕円体と見なして距離を求めるもの
    • distMeeus()
    • distGeo()
    • distVincentyEllipsoid()

地球は真球というよりも楕円体に近いので、楕円体と見なしたものの方が正確な値が出る。また、楕円体と見なすものの中でも、distGeo()の方が distMeeus() より正確であるという話がこのパッケージのドキュメンテーションに載っている。

geosphere パッケージの函数を用いた東京とロンドンの間の距離（単位：メートル）
函数名 距離
distCosine 9572412
distVincentySphere 9572412
distHaversine 9572412
distMeeus 9585294
distGeo 9585300
distVincentyEllipsoid 9585300

距離行列

距離行列を作りたければ、geosphere パッケージに含まれているdistm() を使えば良い。以下の例では、東京・ロンドン・リマの3地点についてそれぞれの間の距離を distGeo() 函数を用いて求めている ((実は、fun オプションに何も指定しなければ、デフォルトで distGeo() が用いられる。)) 。

distm(rbind(tokyo, lima, london), fun = distGeo)

これで、以下のような距離行列が得られる。

         [,1]     [,2]     [,3]
[1,]        0 15494082  9585300
[2,] 15494082        0 10160441
[3,]  9585300 10160441        0

結果をもう少し分かりやすくしたければ、行名や列名を適宜つけよう。

cities <- rbind(tokyo, lima, london)
city_names <- stringr::str_to_title(rownames(cities))
distances <- distm(cities, fun = distGeo)
colnames(distances) <- city_names
rownames(distances) <- city_names

こうやってできた distances という距離行列の中身は以下のようになっている。

          Tokyo     Lima   London
Tokyo         0 15494082  9585300
Lima   15494082        0 10160441
London  9585300 10160441        0

方位の計算

ある地点から別の地点への方位角を求める際には、bearing()という関数を用いる。結果は角度で返される。真東が90度、真北が0度、真西が-90度、真南が-180度になる。

例えば、東京からロンドンへの方位角を求めるには、次のようにする。

bearing(tokyo, london)

すると、-23.69978 という結果が得られ、東京から見たときにロンドンはおおよそ北北西方向にあることが分かる。

このパッケージの紹介記事の続きは「Rの geosphere パッケージ(2)：大圏航路を描く」を参照のこと。

注釈