01 февраля 2014

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



В 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

2 комментария :

Анонимный комментирует...

Сергей,
спасибо за информацию. Функции lower.tri() и upper.tri() извлечь диагональ со значениями не позволили, но зато натолкнули на правильный ход мыслей. В результате решил задачу самостоятельно: оказывается есть отдельная функция для извлечения диагонали - diag(). Как всегда, всё гениальное просто:))
Странно, отслеживаю все Ваши новые записи, а последнюю пропустил. Хотя именно там и была нужная информация.

У меня появился вопрос общего характера. Как научится находить нужные функции? Понятное дело, что есть справка, однако, она требует конкретного запроса. Есть ли у Вас какие-то рекомендации по этому поводу?

Sergey Mastitsky комментирует...

Ответ на вопрос по поводу поиска функций можно найти на все том же ресурсе, т.е. здесь в блоге :)
См. http://r-analytics.blogspot.de/2012/10/sos-r.html

Отправить комментарий