Aвтор: Владимир Шитиков
Введение
Целью
многих исследований является анализ конфигурации пространственных объектов и отображение ее структуры на картосхемах. В среде R для решения этой задачи часто используют методы и функции пакетов
sp, maps, RgoogleMaps и связанных с ними ресурсов. При этом для создания основного
слоя карты часто применяются свободно распространяемые shape-файлы, содержащие точечные и контурные графические примитивы,
соответствующие отдельным географическим пунктам или регионам (см. пример здесь и здесь). Однако, по сравнению с картами, полученными на основе
специализированных геоинформационных систем, таких как ArcGIS ESRI, подобная
визуализация может показаться не столь симпатичной.
Ниже рассматриваются некоторые "продвинутые" методы для быстрой визуализации
пространственных данных в среде R, основанные на двух идеях: (1) формирование "на лету" статических карт необходимого качества и масштаба с использованием
актуальной информации серверов GoogleMap, OpenStreetMap, Stamen Maps или
CloudMade и (2) широкое использование грамматики создания графических слоев для
отображения необходимой информации на основе функций пакета ggplot2 (Wickham 2009, 2016; Мастицкий 2016). В результате развития этих концепций был разработан
удобный пакет R, названный ggmap (Kahle & Wickham 2013). Продемонстрируем
некоторые его возможности с использованием собственных данных.
Создание основного картографического слоя
Формирование основого слоя карты в ggmap состоит из двух этапов. На первом этапе с использованием
функции get_map() осуществляется загрузка из Интернета необходимых
картографических данных, их форматирование для создания наиболее подходящего изображения и формирование объекта ggmap, включающего полученный растр и связанные с ним атрибуты. На
втором этапе с использованием функции qmap() осуществляется вывод на
задаваемое графическое устройство основы карты, описанной объектом ggmap.
Источник данных, из которого осуществляется загрузка картоосновы, определяется при вызове функции get_map() параметром source = c("google", "osm", "stamen", "cloudmade"), значения которого соответствуют вышеперечисленным серверам. Кроме того, для загрузки из каждого ресурса можно использовать также специальные функции get_googlemap(), get_openstreetmap(), get_stamenmap() или get_cloudmademap(), для которых список используемых параметров становится более точно определенным. Все они возвращают объект ggmap стандартизованной структуры, однако далее мы будем рассматривать только карты, загруженные с сервисов GoogleMap и Stamen Maps, поскольку использование OpenStreetMap и CloudMade требуют установки дополнительных модулей или предварительной регистрации на сервере и получения доступа к соответствующим API.
Географическая локализация карты выполняется с
использованием двух параметров: location и zoom. Для get_googlemap() стандартом location является пара координат "долгота/широта", определяющих центр карты. Величина радиуса пространственного охвата вокруг
этого центра определяется аргументом масштабирования zoom, который является целым числом от 3 до 20. При zoom = 3 получаем примерно уровень континента, а при zoom = 20 имеем масштаб отдельно стоящих зданий.
Третьим важным аргументом является maptype, который определяет набор загружаемых с сервера графических компонентов и их эстетику. Для GoogleMap возможны следующие
аргументы maptype = c("terrain",
"satellite", "roadmap", "hybrid"). Рассмотрим пример создания основного слоя карты,
выводящей фрагмент территории г. Саратова, на улицах которого, как уверяют, "огней так много золотых". Координаты населенного пункта определим с использованием
функции geocode():
library(ggmap) (reg_center = geocode("Saratov")) lon lat 1 45.9608 51.59237 ggmap(get_map(location = as.numeric(reg_center), source = "google", zoom = 14, maptype = "roadmap"))
После подготовки изображения основного картографического слоя
можно приступить к нанесению дополнительных слоев данных, а также различных статистик
или модельных объектов, полученных с использованием функций ggplot2.
Отображение точечных данных
В качестве примера будем использовать базу
гидробиологических данных Института экологии Волжского бассейна РАН. В ходе
экспедиционных исследований было взято 1158 проб донных сообществ
беспозвоночных организмов в 100 малых и средних реках Среднего и Нижнего
Поволжья, после чего результаты гидробиологической съемки были загружены в таблицы базы
Microsoft Access. Покажем, как можно работать в среде R с данными,
представленными в подобном формате. К сожалению, предоставление этой базы данных в публичный доступ невозможно, что не позволит читателям воспроизвести приведенные ниже примеры. Тем не менее, будем надеятся, что эти примеры окажутся полезными, и читатель легко сможет "спроецировать" их на собственные задачи.
Соединение с базой данных выполним через стандартный API-интерфейс Open DatabaseConnectivity, для чего первым делом создадим новый объект ODBC, используя панель администрирования операционной
системы Windows. Присвоим объекту ODBC имя "База", свяжем его с актуальным
файлом MDB и определим необходимый драйвер. Соединение с источником данных
выполним при помощи функции odbcConnect() из пакета RODBC:
library(RODBC) con <- odbcConnect("База")
Таблица "Координаты" нашей базы данных содержит
географические координаты рек и их отдельных участков (для средних рек выделены
отдельно истоки, среднее и нижнее течения и устье). Выведем на картосхему
названия обследованных рек. Для создания основного слоя карты используем опцию maptype
= "satellite":
SQL = "SELECT Наименование, Avg(Широта) AS lat, Avg(Долгота) AS lon FROM Координаты GROUP BY Наименование" koord <- sqlQuery(con, SQL) myMap_gs <- ggmap(get_googlemap(center = c(50, 51.4), zoom = 6, maptype = "satellite", extent = "device")) myMap_gs + geom_text(aes(x = lon, y = lat, label = Наименование), size = 2.7, colour = "green", data = koord) + ggtitle("Обследованные реки")
Теперь выведем точками на карте координаты мест отбора проб
из таблицы "Координаты", взяв за основу карту из сервиса Stamen
Maps. Отметим, что при использовании карт этого типа необходимо установить пакет ggmap не из репозитория CRAN, а из Github-репозитория его
разработчиков (подробности см. здесь
):library(devtools) install_github("dkahle/ggmap") # -------------------------------- # Участки отбора проб SQL = "SELECT * FROM [Координаты]" koord <- sqlQuery(con, SQL) bbox <- c(left = 45, bottom = 48.5, right = 55, top = 55) myMap_st <- ggmap(get_stamenmap(bbox, maptype = "terrain")) myMap_st + geom_point(aes(x = Долгота, y = Широта), data = koord) + ggtitle("Точки отбора проб")
Обратим внимание, что локализация карты в этом случае может
выполняться иным, часто более удобным способом: с указанием точных
географических координат углов выделяемого прямоугольника поверхности Земли.
Кроме того, для карт Stamen Maps имеется возможность использовать различные
опции параметра maptype = c("terrain",
"terrain-background", "terrain-labels",
"terrain-lines", "toner", "toner-2010",
"toner-2011", "toner-background", "toner-hybrid",
"toner-labels", "toner-lines", "toner-lite",
"watercolor"). Использование стиля toner приводит к резко контрастным черно-белым картам "газетного" вида, а карты watercolor, предназначенные для
визуализации озерных систем, изображаются двумя группами оттенков цвета
(желто-коричневого для суши и голубого для водных объектов). Мы предпочитаем
использовать terrain или terrain-background, где
подавляются все надписи, часто только мешающие анализу карт.
В ходе гидробиологической съемки было обнаружено более 700
видов и более крупных таксонов водных беспозвоночных. Выгрузим из базы показатель Count_Sample, определяющий, сколько раз тот или иной вид встретился в
процессе исследований на каждом из участков рек. Например, рассмотрим, как
распространены в реках Поволжья личинки кровососущих
мокрецов рода Culicoides. Используем
для отображения на карте кружки, диаметр которых пропорционален Count_Sample:
Species = 'CeCul.sp' SQL = paste("SELECT * FROM [_Численность с координатами] WHERE Code = '", Species, "'", sep = "") sumdat <- sqlQuery(con, SQL) name_sp <- as.character(sumdat$NAME.1[1]) head(sumdat[, c(1, 10, 11, 6, 8, 9)]) NAME Широта Долгота NAME.1 Count_Sample Sum_CHISL 1 Анлы 53.93778 52.33430 Culicoides sp. 1 50 2 Б. Саморода 49.11437 46.84341 Culicoides sp. 19 22140 3 Б.Черемшан 54.30345 52.02079 Culicoides sp. 1 60 4 Б.Черемшан 54.50150 50.39124 Culicoides sp. 3 100 5 Байтуган 54.16788 52.28791 Culicoides sp. 11 2886 6 Домашка 52.93876 50.60469 Culicoides sp. 1 10
myMap_st + geom_point(aes(x = Долгота, y = Широта, size = Count_Sample), data = sumdat) + ggtitle(name_sp)
Распределение плотности вероятностей встречаемости видов
Для построения 2D-диаграммы плотности вероятностей выгрузим
из базы данных полный список и координаты проб, в которых встретились виды Culicoides:
SQL = paste("SELECT * FROM [_Проба с координатами] WHERE Code = '", Species, "'", sep = "") pundat <- sqlQuery(con, SQL) name_sp <- as.character(pundat$NAME[1]) nrow(pundat) [1] 222
Графически плотность двухмерного распределения этих точек
на карте можно отобразить, используя контурные линии, которые ограничивают
участки с одинаковой плотностью вероятности. В пакете ggplot2 для нахождения границ контуров плотности служит функция stat_density2d(), которая в свою очередь основана на функции kde2d() из базового R-пакета MASS. Функция scale_fill_gradient() определяет диапазон цветов заливки
каждого из bins = 5 контуров, а
функция scale_alpha() управляет уровнем прозрачности в
зависимости от градации встречаемости вида:
myMap_st + ggtitle(name_sp) + stat_density2d(data = pundat, aes(x = Долгота, y = Широта, fill = ..level.., alpha = ..level..), bins = 5, geom = "polygon") + scale_fill_gradient(low = "green", high = "red") + scale_alpha(range = c(0, 0.3), guide = FALSE)
Предположим теперь, что мы хотим внимательнее рассмотреть
ситуацию с распространением Culicoides в районе оз. Эльтон (сгущение в
нижней части карты). В дополнение к предыдущему скрипту в отдельном слое
расставим местоположения точек отбора проб и добавим легенду с градациями встречаемости личинок мокреца.
ggmap(get_googlemap(center = c(46.64, 49.21), zoom = 10, maptype = "hybrid", extent = "device")) + ggtitle(name_sp) + geom_point(data = pundat, aes (x = Долгота, y = Широта), colour = "red") + stat_density2d(data = pundat, aes(x = Долгота, y = Широта, fill = ..level.., alpha = ..level..), bins = 5, geom = "polygon") + scale_fill_gradient("Плотность\nвстречаемости", low = "green", high = "red") + guides(fill = guide_colorbar(barwidth = 1.5, barheight = 10)) + scale_alpha(range = c(0, 0.3), guide = FALSE)
Литературные источники
- Wickham H. (2009, 2016) ggplot2: Elegant Graphics for Data Analysis. Springer, New York. URL: http://ggplot2.org/book/
- Мастицкий С. Э. (2016) Визуализация данных с помощью ggplot2. М.: ДМК Пресс, 2016
- Kahle D., Wickham H. (2013) ggmap: Spatial Visualization with ggplot2. The R Journal 5: 144–162. URL: http://vita.had.co.nz/papers/ggmap.html
Отправить комментарий