Расчет расстояния между географическими объектами по их координатам

В R имеется широкий арсенал инструментов для работы с пространственными данными. Со списком основных из них можно ознакомиться в соответствующем разделе на сайте CRAN. Это сообщение будет посвящено простой, но в то же время часто возникающей задаче - расчету расстояния между географическими объектами по их координатам.

Предположим, что по каким-то причинам нас интересуют четыре местоположения на территории... острова Мадагаскар (я всегда хотел побывать на этом острове - отсюда и выбор примера), показанные ниже в виде точек красного цвета:

Приведенная карта была построена при помощи следующего кода R (подробнее об использовании графического пакета ggplot2 для работы с картами см. здесь):

library(sp)
library(maptools)
library(ggplot2)
 
load(url("http://gadm.org/data/rda/MDG_adm0.RData"))
 
madagascar <- fortify(gadm)
 
m <- ggplot() + geom_map(data = madagascar, 
                    aes(map_id = id), 
                    map = madagascar,
                    fill = "white", color = "black") +
  expand_limits(x = madagascar$long, y = madagascar$lat) +
  coord_map("mercator") + 
  xlab("Долгота") + ylab("Широта") + theme_bw()
 
mp <- m + geom_point(data = data.frame(Lon = c(45.0, 47.2, 45.5, 48.5),
                                       Lat = c(-25.0, -21.0, -17.0, -15.0)),
                     aes(Lon, Lat), color = I("red"), size = 3)
 
print(mp)

Для начала рассчитаем расстояние между двумя наиболее удаленными точками, т.е.:

Возможности для расчета расстояния между географическими объектами по их координатам имеются в нескольких пакетах R. Мы воспользуемся пакетом geosphere, в состав которого входит большой набор функций, реализующих методы пространственной тригонометрии. В частности, в этом пакете имеется несколько функций для расчета кратчайшего расстояния между двумя точками на поверхности Земли разными методами (distCosine(), distVincentySphere(), distVincentyEllipsoid(), distHaversine(), distMeeus()). Применим функцию distHaversine(), которая рассчитывает расстояние по широко используемой формуле гаверсинуса (см. также здесь). Эта функция имеет следующий синтаксис:

distHaversine(p1, p2, r = 6378137),

где p1 - координаты точки (или точек), подаваемые либо в виде вектора из двух значений (долгота, широта), либо в виде матрицы с двумя столбцами (первый столбец - долгота, второй - широта); p2 - то же, что и р1; r = радиус Земли, выраженный в метрах.

Для нашего примера получаем:

install.packages("geosphere")
library(geosphere)

# делим на 1000 для выражения расстояния в км
distHaversine(c(45.0, -25.0), c(48.5, -15.0))/1000 
[1] 1171.7

Следует отметить, что полученный результат основан на предположении, что Земля имеет форму шара. Однако известно, что это на совсем верно - Земля представляет собой эллипсоид. При расчете расстояния между существенно удаленными объектами это обстоятельство может сказаться на точности результата. В таких случаях стоит воспользоваться функцией distVincentyEllipsoid().

Теперь несколько усложним задачу и предположим, что требуется рассчитать расстояния между всеми имеющимися парами точек. Необходимость такого рода расчетов часто возникает, например, в экологии при анализе пространственного распределения особей того или иного вида, когда требуется вычислить среднее расстояние между особями в целом.

Сохраним координаты имеющихся точек в виде матрицы:

coords <- cbind(c(45.0, 47.2, 45.5, 48.5),
                c(-25.0, -21.0, -17.0, -15.0))
 
coords
     [,1] [,2]
[1,] 45.0  -25
[2,] 47.2  -21
[3,] 45.5  -17
[4,] 48.5  -15

Теперь мы можем легко рассчитать все пары расстояний, воспользовавшись базовой R-функцией apply() (подробнее см. здесь):

Dist <- apply(coords, 1, FUN =
              function(eachPoint) distHaversine(eachPoint, coords)/1000)
 
Dist
          [,1]     [,2]     [,3]      [,4]
[1,]    0.0000 499.0594 892.0671 1171.6513
[2,]  499.0594   0.0000 479.8662  681.9332
[3,]  892.0671 479.8662   0.0000  390.6508
[4,] 1171.6513 681.9332 390.6508    0.0000

Как видим, в итоговой квадратной матрице информация о расстояниями между точками дублируется (по обе стороны от главной диагонали). Является ли это проблемой, зависит от стоящей задачи. Например, при выполнении кластерного анализа, многие соответствующие функции примут такую матрицу как есть. Но при необходимости вычислить среднее расстояние (или какие либо другие статистики) потребуются только те значения, которые расположены сверху или снизу от главной диагонали матрицы. Извлечение этих значений не представляет проблемы: для этого необходимо воспользоваться одной из базовых R-функций - lower.tri() или upper.tri(). Для нашего примера получаем:

Dist <- Dist[lower.tri(Dist)]
 
Dist
[1]  499.0594  892.0671 1171.6513  479.8662  681.9332  390.6508

Рассчитать среднее расстояние теперь очень легко:

mean(Dist)
[1] 514.4

--

Другие сообщения по теме:
- Создание пользовательских карт Google
- Создание картограмм: часть 1 и часть 2