【R】geosphere

1. はじめに

geosphereは、geospatial dataの距離や方向を解析したり表示できるパッケージです。

2. インストール

install.packages("geosphere")

3. つかってみる

まずは、ある位置（2点間）の地球表面上のルートを表示します。

library(geosphere)
library(tidyverse)

LA <- c(-118.40, 33.95)
DE <- c(20, 55)
data(wrld)
plot(wrld, type="l")
gc <- greatCircle(LA, DE)
lines(gc, lwd=1, col="yellow")
gci <- gcIntermediate(LA, DE)
lines(gci, lwd=4, col="green")
points(rbind(LA, DE), col="red" , pch=20, cex=2)

３つのルートで囲まれた三角形を表示します。

MS <- c(-93.26, 44.98)
gc1 <- greatCircleBearing(NY, 281)
gc2 <- greatCircleBearing(MS, 195)
gc3 <- greatCircleBearing(LA, 55)
plot(wrld, type= "l" , xlim=c(-125, -70), ylim=c(20, 60))
lines(gc1, col= "green" )
lines(gc2, col= "blue" )
lines(gc3, col= "red" )
int <- gcIntersectBearing(rbind(NY, NY, MS), c(281, 281, 195), rbind(MS, LA, LA), c(195, 55, 55))
int

distm(rbind(int[,1:2], int[,3:4]))
int <- int[,1:2]
points(int)
poly <- rbind(int, int[1,])
centr <- centroid(poly)
poly2 <- makePoly(int)
polygon(poly2, col= "yellow" )
points(centr, pch= "*" , col= "dark red" , cex=2)

> int
            lon      lat      lon       lat
[1,]  -94.40975 41.77229 85.59025 -41.77229
[2,] -102.91692 41.15816 77.08308 -41.15816
[3,]  -93.63298 43.97765 86.36702 -43.97765

> distm(rbind(int[,1:2], int[,3:4]))
           [,1]       [,2]       [,3]       [,4]       [,5]
[1,]        0.0   713665.0   253077.7 20003931.5 19308947.3
[2,]   713665.0        0.0   823532.8 19308947.3 20003931.5
[3,]   253077.7   823532.8        0.0 19751854.8 19195862.0
[4,] 20003931.5 19308947.3 19751854.8        0.0   713665.0
[5,] 19308947.3 20003931.5 19195862.0   713665.0        0.0
[6,] 19751854.8 19195862.0 20003931.5   253077.7   823532.8
           [,6]
[1,] 19751854.8
[2,] 19195862.0
[3,] 20003931.5
[4,]   253077.7
[5,]   823532.8
[6,]        0.0

4. さいごに

geospatialなデータから、距離や方位を計算してくれ、さらに表示してくれるパッケージは貴重ですね。