R言語 + Leafletで大圏航路

 過去2回(Python + BasemapおよびR言語 + Leaflet)、秋田-ハワイ、阿武-グアムを直線で結んでみた。

 メルカトル図法の地図では距離、面積が正しくないので、2地点を結んだ直線は最短距離にならないが、放射型の正距図法の地図であれば最短距離になる。いわゆる大圏航路。

 正距図法の地図上で2地点を直線で結んで、緯線・経線との交点(緯度・経度)をメルカトル図法の地図上にプロットしていくと大圏航路が描かれる。

 Python + Basemapであれば、BasemapライブラリのHPにあるサンプル(Great Circle from NY to London)で使われているdrawgreatcircle()が、その関数。
 ただし、秋田-ハワイのように経度180°をまたいでいる場合、うまく線が引けない。
 Basemapをバージョンアップすれば線が引ける……。
 ところが、ANACONDAトラブルへと発展。

 R言語にシフトしつつある今日この頃。

 R言語の場合、geosphereパッケージのgcIntermediate()で大圏航路を描ける。

 ※ Colorless Green Ideasのページ(id.fnshr.info/2018/05/20/geosphere2/)を参照。

library(geosphere)
guam <- c(144.75,13.47)
abu <- c(131.47,34.50)
hawaii <- c(-157.86,21.31)
akita <- c(140.10,39.72)

の後、

gcl_ga <- gcIntermediate(guam,abu,addStartEnd = TRUE)

で、阿武-グアムの大圏航路(の行列)が作られる。
 addStartEnd = TRUEをつけないと始点・終点が結果に含まれない。

 ただし、経度180°をまたいでいる秋田-ハワイの間は、うまく線が引けない。

 見覚えのある東西横断線。

 上記Colorless Green IdeasのHPでは、

gcIntermediate(hawaii,akita,addStartEnd = TRUE,breakAtDateLine = TRUE)

といった具合にbreakAtDateLine = TRUEで行列を経度180°までと経度-180°< <0°の2つに分けて、lapply()を使って順次描画させている。

 が、…… plotではなくleafletパッケージを使って描画するところで、つまずいた。

 データ数が多くないので、
 『Leafletで世界地図 Natural Earth
でやったように経度-180°< <0°の地域を360°プラスして経度180°< <360°扱いに変換してみる。

gcl_ha <- gcIntermediate(hawaii,akita,addStartEnd = TRUE)
for(i in 1:nrow(gcl_ha)){if (gcl_ha[i,1] < 0) {gcl_ha[i,1] <- gcl_ha[i,1] + 360}}

 これで一応OK。
 nrow(gcl_ha)はgcl_ha(行列)の行数。

 leafletパッケージでの描画は、

library(leaflet)
map <- leaflet()
map <- addTiles(map)
map <- setView(map,150.00,34.50,3)
map <- addPolylines(map,data=gcl_ga,color="orange",weight=3)
map <- addPolylines(map,data=gcl_ha,color="red",weight=3)
map

 ついでに直線もグレーで加えて描画。
 ちょっとカクカクしている。

R + Leafletで大圏航路   

 

MAP R
スポンサーリンク
ふシゼン
タイトルとURLをコピーしました