В 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 = радиус Земли, выраженный в метрах.
Для нашего примера получаем:
Следует отметить, что полученный результат основан на предположении, что Земля имеет форму шара. Однако известно, что это на совсем верно - Земля представляет собой эллипсоид. При расчете расстояния между существенно удаленными объектами это обстоятельство может сказаться на точности результата. В таких случаях стоит воспользоваться функцией distVincentyEllipsoid().
Теперь несколько усложним задачу и предположим, что требуется рассчитать расстояния между всеми имеющимися парами точек. Необходимость такого рода расчетов часто возникает, например, в экологии при анализе пространственного распределения особей того или иного вида, когда требуется вычислить среднее расстояние между особями в целом.
Сохраним координаты имеющихся точек в виде матрицы:
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().
Теперь несколько усложним задачу и предположим, что требуется рассчитать расстояния между всеми имеющимися парами точек. Необходимость такого рода расчетов часто возникает, например, в экологии при анализе пространственного распределения особей того или иного вида, когда требуется вычислить среднее расстояние между особями в целом.
Сохраним координаты имеющихся точек в виде матрицы:
Теперь мы можем легко рассчитать все пары расстояний, воспользовавшись базовой R-функцией apply() (подробнее см. здесь):
Как видим, в итоговой квадратной матрице информация о расстояниями между точками дублируется (по обе стороны от главной диагонали). Является ли это проблемой, зависит от стоящей задачи. Например, при выполнении кластерного анализа, многие соответствующие функции примут такую матрицу как есть. Но при необходимости вычислить среднее расстояние (или какие либо другие статистики) потребуются только те значения, которые расположены сверху или снизу от главной диагонали матрицы. Извлечение этих значений не представляет проблемы: для этого необходимо воспользоваться одной из базовых R-функций - lower.tri() или upper.tri(). Для нашего примера получаем:
Рассчитать среднее расстояние теперь очень легко:
--
Другие сообщения по теме:
- Создание пользовательских карт Google
- Создание картограмм: часть 1 и часть 2
Как видим, в итоговой квадратной матрице информация о расстояниями между точками дублируется (по обе стороны от главной диагонали). Является ли это проблемой, зависит от стоящей задачи. Например, при выполнении кластерного анализа, многие соответствующие функции примут такую матрицу как есть. Но при необходимости вычислить среднее расстояние (или какие либо другие статистики) потребуются только те значения, которые расположены сверху или снизу от главной диагонали матрицы. Извлечение этих значений не представляет проблемы: для этого необходимо воспользоваться одной из базовых 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
спасибо за информацию. Функции lower.tri() и upper.tri() извлечь диагональ со значениями не позволили, но зато натолкнули на правильный ход мыслей. В результате решил задачу самостоятельно: оказывается есть отдельная функция для извлечения диагонали - diag(). Как всегда, всё гениальное просто:))
Странно, отслеживаю все Ваши новые записи, а последнюю пропустил. Хотя именно там и была нужная информация.
У меня появился вопрос общего характера. Как научится находить нужные функции? Понятное дело, что есть справка, однако, она требует конкретного запроса. Есть ли у Вас какие-то рекомендации по этому поводу?
См. http://r-analytics.blogspot.de/2012/10/sos-r.html
Отправить комментарий