28  Основы ГИС

28.1 Обзор

Пространственные аспекты ваших данных могут дать много полезной информации в ситуации вспышек, а также ответить на подобные вопросы:

  • Где находятся текущие очаги заболевания?
  • Как изменились очаги со временем?
  • Каков доступ к медицинским организациям? Нужны ли какие-то улучшения?

Фокус данной страницы по ГИС - ответить на потребности прикладных эпидемиологов при реагировании на вспышки. Мы рассмотрим основные методы пространственной визуализации данных, используя пакеты tmap и ggplot2. Мы также рассмотрим некоторые основные методы управления данными и запросов с помощью пакета sf. И наконец, мы кратко затронем концепции пространственной статистики, например, пространственные отношения, пространственная автокорреляция и пространственная регрессия, используя пакет spdep.

28.2 Ключевые понятия

Ниже мы представим некоторые ключевые понятия. Для более подробного знакомства с ГИС и пространственным анализом, мы рекомендуем вам изучить один из более подробных самоучителей или курсов, указанных в разделе Ресурсы.

Географическая информационная система (ГИС) - ГИС - это рамка или среда для сбора, управления, анализа и визуализации пространственных данных.

Программное обеспечение ГИС

Некоторые популярные ГИС программы допускают взаимодействие по принципу “навести и кликнуть” для разработки карты и пространственного анализа. Эти инструменты имеют свои преимущества - не нужно учиться писать код и легко выбирать вручную и размещать иконки и параметры на карте. Вот две популярных программы:

ArcGIS - Коммерческая ГИС программа, разработанная компанией ESRI, которая очень популярна, но весьма дорогая

QGIS - Бесплатная открытая ГИС программа, которая может делать почти все то же самое, что и ArcGIS. Вы можете скачать QGIS тут

Использование R как ГИС может сначала показаться сложным, поскольку вместо принципа “навести и кликнуть” используется “интерфейс командной строки” (необходимо написать код, чтобы получить желаемый результат). Однако это является серьезным преимуществом, если вам необходимо регулярно готовить одни и те же карты или создавать воспроизводимый анализ.

Пространственные данные

Две основные формы пространственных данных, используемых в ГИС, включают в себя векторные и растровые данные:

Векторные данные - Наиболее частый формат пространственных данных, используемых в ГИС, векторные данные состоят из геометрических характеристик вершин и путей. Векторные пространственные данные можно еще разделить на три широко используемых типа:

  • Точки - Точка состоит из пары координат (x,y), представляющих определенное место в системе координат. Точки являются наиболее простой формой пространственных данных и могут использоваться для обозначения случая (например, дома пациента) или места (например, больницы) на карте.

  • Линии - Линия состоит из двух соединенных точек. Линии имеют определенную длину и могут использоваться для обозначения таких объектов, как дороги или реки.

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

Растровые данные - Альтернативный формат пространственных данных, растровые данные представляют собой матрицу ячеек (например, пикселей), каждая из которых содержит такую информацию, как высота, температура, уклон, лесной покров и т.д. Часто это аэрофотоснимки, спутниковые снимки и т.д. Растровые данные также могут использоваться в качестве “базовых карт” под векторными данными.

Визуализация пространственных данных

Чтобы визуально представить пространственные данные на карте, программное обеспечение ГИС требует от вас задать достаточную информацию о том, где должны быть характеристики по отношению друг к другу. Если вы используете векторные данные, что действительно так в большинстве примеров применения, эта информация, как правило, хранится в шейп-файле:

Шейп-файлы - Шейп-файл - частый формат данных для хранения “векторных” пространственных данных, состоящих из линий, точек или полигново. Один шейп-файл на самом деле является сочетанием как минимум трех файлов - .shp, .shx, и .dbf. Все эти файлы-подкомпоненты должны присутствовать в определенной директории (папке) для того, чтобы шейп-файл был читаемым. Эти связанные файлы можно сжать в ZIP папку и направить по почте, либо скачать с веб-сайта.

Шейп-файл будет содержать информацию о самих характеристиках, а также о том, где их найти на поверхности Земли. Это важно, потому что хотя Земля - шар, карты, как правило, двухмерные; решения о том, как сделать пространственные данные “плоскими” может иметь большое влияение на вид и интерпретацию полученных в результате карт.

Референтная система координат (РСК) - это это система координат, используемая для определения местоположения географических объектов на поверхности Земли. Она состоит из нескольких основных компонентов:

  • Ситема координат - Существует множество различных систем координат, поэтому необходимо знать, в какой системе находятся ваши координаты. Обычно используются градусы широты/долготы, но можно также использовать координаты универсальной поперечной проекции Меркатора UTM.

  • Единицы - знайте, какие единицы используютсядля вашей системы координат (например, десятичные градусы, метры)

  • Основа системы координат - Определенная смоделированная версия Земли. В течение многих лет они пересматривались, поэтому убедитесь, что слои карты используют одну и ту же основу.

  • Проекция - Ссылка на математическое уравнение, которое было использовано для проецирования истинно круглой Земли на плоскую поверхность (карту).

Помните, что обобщить пространственные данные можно и без использования представленных ниже картографических средств. Иногда достаточно простой таблицы по географическому признаку (например, район, страна и т.д.)!

28.3 Начало работы с ГИС

Для создания карты необходимо иметь и продумать несколько ключевых моментов. К ним относятся:

  • Набор данных – Это может быть формат пространственных данных (например, шейп-файлы, как отмечалось выше), а может быть и не пространственный формат (например, просто csv).

  • Если ваш набор данных не в пространственном формате, вам также потребуется референтный набор данных. Референтные данные состоят из пространственного представления данных и соответствующих атрибутов, которые включают материал, содержащий информацию о месте расположения и адресе конкретных характеристик.

    • Если вы работаете с заранее определенными географическими границами (например, административными районами), то референтные шейп-файлы часто можно свободно загрузить из государственных учреждений или организаций по обмену данными. В случае сомнений лучше всего начать с поиска в Google “шейп-файл [регионов]”

    • Если у вас есть информация об адресе, но нет широты и долготы, вы можете использовать механизм геокодирования, чтобы получить референтные пространственные данные для ваших записей.

  • Идея того, как вы хотите представить информацию в ваших наборах данных вашей целевой аудитории. Существует много разных типов карт и важно подумать о том, какой тип карты лучше подходит для ваших потребностей.

Типы карт для визуализации ваших данных

Картограмма - тип тематической карты, в которой цвета, тени или узоры используются для представления географических регионов в зависимости от значения какого-либо атрибута. Например, большее значение может быть обозначено более темным цветом, чем меньшее. Этот тип карт особенно полезен при визуализации переменной и ее изменений в определенных регионах или геополитических зонах.

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

Карта плотности точек - тип тематической карты, на которой точками обозначаются значения атрибутов в данных. Этот тип карты лучше всего использовать для визуализации разброса данных и визуального поиска кластеров.

Карта пропорциональных символов (карта градуированных символов) - тематическая карта, аналогичная картограмме, но вместо цвета для обозначения значения атрибута используется символ (обычно круг) по отношению к значению. Например, большее значение может быть обозначено более крупным символом, чем меньшее значение. Этот тип карты лучше всего использовать, когда необходимо визуализировать размер или количество данных по географическим регионам.

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

Примечание: Основной фокус данной страницы по ГИС основан на контексте реагирования на вспышки заболеваний на местах. Поэтому содержание страницы охватывает основные операции с пространственными данными, визуализацию и анализ.

28.4 Подготовка

Загрузка пакетов

Этот фрагмент кода показывает загрузку пакетов, необходимых для анализа. В данном руководстве мы фокусируемся на использовании p_load() из пакета pacman, которая устанавливает пакет, если необходимо, и загружает его для использования. Вы можете также загрузить установленные пакеты с помощью library() из базового R. См. страницу [Основы R] для получения дополнительной информации о пакетах R.

pacman::p_load(
  rio,           # для импорта данных
  here,          # для пути к файлу
  tidyverse,     # для вычистки, работы с данными и построения графиков (включает пакет ggplot2)
  sf,            # для управления пространственными данными, используя формат простых свойств (Simple Feature)
  tmap,          # для создания простых карт, работает как для интерактивных, так и статичных карт
  janitor,       # для вычистки имен столбцов
  OpenStreetMap, # для добавления базовой карты OSM на карту ggplot
  spdep          # пространственная статистика
  ) 

Вы можете посмотреть обзор всех пакетов, которые работают с пространственными данными в CRAN “Spatial Task View”.

Пример данных о случаях

Для демонстрационных целей мы будем работать со случайной выборкой из 1000 случаев из имитированной эпидемии Эболы из датафрейма linelist (с точки зрения вычислений работу с меньшим количеством случаев легче отобразить в руководстве). Если вы хотите выполнять действия параллельно, кликните, чтобы скачать “чистый” построчный список (как .rds файл).

Поскольку мы делаем случайную выборку, при самостоятельном выполнении кодов результаты могут несколько отличаться от представленных здесь.

Импортируем данные с помощью функции import() из пакета rio (он работает с разными типами файлов, такими как .xlsx, .csv, .rds - см. детали на странице [Импорт и экспорт]).

# импорт чистого построчного списка случаев
linelist <- import("linelist_cleaned.rds")  

Далее делаем случайную выборку 1000 строк, используя sample() из базового R.

# генерируем 1000 случайных номеров строк из числа строк построчного списка
sample_rows <- sample(nrow(linelist), 1000)

# подмножество построчного списка, чтобы сохранить только строки выборки и все столбцы
linelist <- linelist[sample_rows,]

Теперь нам нужно конвертировать этот построчный список linelist, который относится к классу датафрейм, в объект класса “sf” (пространственные свойства). Учитывая, что в построчном списке есть два столбца “lon” и “lat”, представляющих собой долготу и широту для места проживания каждого случая, это будет легко сделать.

Мы используем пакет sf (пространственные свойствуа) и его функцию st_as_sf() для создания нового объекта, который мы назовем linelist_sf. Этот новый объект, по сути, выглядит так же как и построчный список, но столбцы lon и lat были отмечены как столбцы координат и была присвоена референтная система координат при отображении точек. 4326 определяет наши координаты как основанные на всемирной геодезической системе 1984 (WGS84) - которая является стандартом для GPS координат.

# создаем объект sf
linelist_sf <- linelist %>%
     sf::st_as_sf(coords = c("lon", "lat"), crs = 4326)

Вот так выглядит оригинальный датафрейм linelist. В этой демонстрации мы будем использовать только столбец date_onset и geometry (который был создан из полей широты и долготы, указанных выше, и является последним столбцом в датафрейме).

DT::datatable(head(linelist_sf, 10), rownames = FALSE, options = list(pageLength = 5, scrollX=T), class = 'white-space: nowrap' )

Шейп-файлы административных границ

Сьерра-Леоне: шейп-файлы административных границ

Мы заранее скачали все административные границы для Сьерра-Леоне с веб-сайта Humanitarian Data Exchange (HDX) веб-сайт тут. Альтернативно вы можете скачать их и другие примеры данных для руководства через наш пакет R, который представлен на странице [Скачивание руководствау и данных].

Теперь мы выполним следующие действия, чтобы сохранить шейп-файл Admin Level 3 в R:

  1. Импортируем шейп-файл
  2. Вычищаем имена столбцов
  3. Фильтруем строки, чтобы сохранить только интересующие районы

Чтобы импортировать шейп-файл, мы используем функцию read_sf()из sf. Ей задается путь к файлу с помощью here(). - в нашем случае файл находится внутри нашего проекта R в подпапках “data”, “gis” и “shp”, а имя файла - “sle_adm3.shp” (для дополнительной информации см. страницы [Импорт и экспорт] и [Проекты R]). Вам нужно задать свой собственный путь к файлу.

Далее мы используем clean_names() из пакета janitor, чтобы стандартизировать имена столбцов шейп-файла. Мы также используем filter(), чтобы сохранить только строки с admin2name “Western Area Urban” или “Western Area Rural”.

# чистый уровень ADM3
sle_adm3 <- sle_adm3_raw %>%
  janitor::clean_names() %>% # стандартизируем имена столбцов
  filter(admin2name %in% c("Western Area Urban", "Western Area Rural")) # фильтр, чтобы сохранить некоторые районы

Ниже вы видите как выглядит каждый шейп-файл после импорта и вычистки. Пролистайте вправо, чтобы увидеть столбцы с admin уровень 0 (страна), admin уровень 1, admin уровень 2, и наконец, admin уровень 3. У каждого уровня есть текстовое имя и уникальный идентификационный “pcode”. pcode увеличивается на каждом увеличивающемся уровне admin level, например, SL (Сьерра-Леоне) -> SL04 (Western - западная) -> SL0410 (Western Area Rural - Западный сельский район) -> SL040101 (Koya Rural - Сельский район Койа).

Данные о населении

Сьерра-Леоне: Население по ADM3

Опять же, эти данные можно скачать с (ссылка тут) или через наш пакет R epirhandbook, как объясняется [на этой странице][Скачивание руководства и данных]. Мы используем import() для загрузки .csv файла. Мы также передаем импортированный файл в clean_names(), чтобы стандартизировать синтаксис имен столбцов.

# Население по ADM3
sle_adm3_pop <- import(here("data", "gis", "population", "sle_admpop_adm3_2020.csv")) %>%
  janitor::clean_names()

Вот как выглядит файл по населению. Пролистайте вправо, посмотрите, что у каждой юрисдикции есть столбец с мужским населением (male), женским населением (female), общим населением (total), а также разбивка населения по столбцам по возрастной группе.

Медицинские организации

Сьерра-Леоне: Данные о медицинских организациях из OpenStreetMap

Опять же, мы скачиваем места расположения медицинских организаций с HDX тут или с помощью инструкций на странице [Скачивание руководства и данных].

Мы импортируем шейп-файл с точками медицинских организацией с помощью read_sf(), также вычищаем имена столбцов, затем фильтруем, чтобы сохранить только те точки, которые отмечены как “hospital”, “clinic” или “doctors”.

# Шейп-файл по медицинским организациям из OSM
sle_hf <- sf::read_sf(here("data", "gis", "shp", "sle_hf.shp")) %>% 
  janitor::clean_names() %>%
  filter(amenity %in% c("hospital", "clinic", "doctors"))

Вот получившийся в результате датафрейм - пролистайте вправо, чтобы увидеть название организации и координаты geometry.

28.5 Построение графика координат

Самый простой способ построить координаты X-Y (долгота/широта, точки), в случае этих случаев - нарисовать их как точки напрямую из объекта linelist_sf, который мы создали на этапе подготовки.

Пакет tmap предлагает простые возможности составления как статических (режим “график - plot”), так и интерактивных карт (режим “просмотр - view”) буквально с помощью нескольких строк кода. Синтаксис tmap похож на ggplot2 в том, что команды добавляются друг к другу +. Более подробно читайте в этой виньетке.

  1. Устанавливаем режим tmap. В этом случае используем режим графика “plot”, который создает статические результаты.
tmap_mode("plot") # выберите "view" или "plot"

Ниже мы показываем только точки. В tm_shape() задаются объекты linelist_sf. Затем мы добавляем точки с помощью tm_dots(), уточняяя размер и цвет. Поскольку linelist_sf является объектом sf, мы уже назначили два столбца, которые содержат координаты широты/долготы и референтную систему координат (РСК):

# Только случаи (точки)
tm_shape(linelist_sf) + tm_dots(size=0.08, col='blue')

Отдельно точки не очень информативны. Поэтому нужно также наложить карту административных границ:

Мы снова используем tm_shape() (см. документацию), но вместо того, чтобы задать шейп-файл с точками случаев, мы задаем шейп-файл с административными границами (полигонами).

С помощью аргумента bbox = (bbox означает “ограничительная рамка”) мы можем уточнить границы координат. Сначала мы покажем отображение карты без bbox, а затем с ним.

# Только административные границы (полигоны)
tm_shape(sle_adm3) +               # шейп-файл административных границ
  tm_polygons(col = "#F7F7F7")+    # показываем полигоны светло серым
  tm_borders(col = "#000000",      # показываем границы с помощью цвета и толщины линии
             lwd = 2) +
  tm_text("admin3name")            # текст столбца для отображения для каждого полигона


# Также как до этого, но с ограничительной рамкой
tm_shape(sle_adm3,
         bbox = c(-13.3, 8.43,    # угол
                  -13.2, 8.5)) +  # угол
  tm_polygons(col = "#F7F7F7") +
  tm_borders(col = "#000000", lwd = 2) +
  tm_text("admin3name")

Теперь вместе точки и полигоны:

# All together
tm_shape(sle_adm3, bbox = c(-13.3, 8.43, -13.2, 8.5)) +     #
  tm_polygons(col = "#F7F7F7") +
  tm_borders(col = "#000000", lwd = 2) +
  tm_text("admin3name")+
tm_shape(linelist_sf) +
  tm_dots(size=0.08, col='blue', alpha = 0.5) +
  tm_layout(title = "Distribution of Ebola cases")   # определяем заголовок карты

Хорошее сравнение опций картирования в R можно прочитать в этом посте в блоге.

28.6 Пространственные соединения

Вы можете быть знакомы с соединением данных из одного набора данных с другим. На странице [Соединение данных] руководства обсуждается несколько методов. Пространственное соединение служит аналогичным целям, но использует в большей степени пространственные отношения. Вместо того, чтобы полагаться на общие значения в столбцах, которые должны совпадать с наблюдениями, вы можете использовать их пространственные отношения, например, если одна характеристика находится внутри другой, либо является ближайшим соседом к другому, либо находится внутри буферной зоны определенного радиуса от другого и т.п.

Пакет sf предлагает разные методы пространственного соединения. См. дополнительную документацию по методу st_join и типам пространственного соединения по этой ссылке.

Точки в полигонах

Пространственное присваивание административных единиц для случаев

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

Ниже мы отобразим пространственное пересечение расположения случаев (точки) с границами ADM3 (полигоны):

  1. Начинаем с построчного списка (точки)
  2. Пространственное соединение с границами, установка типа соединения в “st_intersects”
  3. Используем select(), чтобы сохранить только некоторые новые столбцы с административными границами
linelist_adm <- linelist_sf %>%
  
  # соединяем файл с административными границами с построчным списком на основе пространственного пересечения
  sf::st_join(sle_adm3, join = st_intersects)

Вса столбцы из sle_adms были добавлены в построчный список! У каждого случая теперь есть столбцы, указывающие административные уровни, к которым он относится. В данном примере мы хотим сохранить только два новых столбца (уровень admin 3), поэтому вы выбираем старые имена столбцов и всего два дополнительных интересующих с помощью select():

linelist_adm <- linelist_sf %>%
  
  # соединяем файл с административными границами с построчным списком на основе пространственного пересечения
  sf::st_join(sle_adm3, join = st_intersects) %>% 
  
  # сохраняем старые имена столбцов и два новых интересующих столбца admin
  select(names(linelist_sf), admin3name, admin3pcod)

Ниже в целях иллюстрирования вы можете увидеть первые 10 случаев и их юрисдикции на уровне admin 3 (ADM3), которые были присоединены на основе того, где точка пространственно пересекалась с формами полигонов.

# теперь используем имена ADM3, присоединенные к каждому случаю
linelist_adm %>% select(case_id, admin3name, admin3pcod)
Simple feature collection with 1000 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -13.27122 ymin: 8.447753 xmax: -13.20613 ymax: 8.490746
Geodetic CRS:  WGS 84
First 10 features:
     case_id     admin3name admin3pcod                   geometry
4280  241a7d       West III   SL040208 POINT (-13.25955 8.456019)
3260  ceee63         East I   SL040203 POINT (-13.21588 8.487489)
3133  501f29        East II   SL040204 POINT (-13.22259 8.483358)
4173  6fbf2f       West III   SL040208 POINT (-13.26315 8.486153)
117   a11d21 Mountain Rural   SL040102 POINT (-13.21382 8.465117)
2619  18ee39       West III   SL040208  POINT (-13.26451 8.45718)
1041  1655de       West III   SL040208 POINT (-13.26031 8.457142)
4641  821030        East II   SL040204 POINT (-13.21894 8.484974)
2865  af82b2       West III   SL040208 POINT (-13.26751 8.461085)
161   8e1819 Mountain Rural   SL040102 POINT (-13.21726 8.465353)

Теперь мы можем описать наши случаи по административной единице - до пространственного соединения мы не могли этого сделать!

# создаем новый датафрейм, содержащий количество случаев по административным единицам
case_adm3 <- linelist_adm %>%          # начинаем с построчного списка с новыми столбцами admin
  as_tibble() %>%                      # конвертируем в таблицу tibble для более удобного отображения
  group_by(admin3pcod, admin3name) %>% # группируем по административной единице по имени и pcode 
  summarise(cases = n()) %>%           # обобщаем и считаем строки
  arrange(desc(cases))                     # сортируем в порядке уменьшения

case_adm3
# A tibble: 10 × 3
# Groups:   admin3pcod [10]
   admin3pcod admin3name     cases
   <chr>      <chr>          <int>
 1 SL040102   Mountain Rural   274
 2 SL040208   West III         218
 3 SL040207   West II          192
 4 SL040204   East II          116
 5 SL040201   Central I         62
 6 SL040203   East I            61
 7 SL040206   West I            30
 8 SL040205   East III          29
 9 SL040202   Central II        16
10 <NA>       <NA>               2

Мы также можем создать столбчатую диаграмму количества случаев по административным единицам.

В этом примере мы начнем ggplot() с linelist_adm, чтобы мы могли применить функции факторов, такие как fct_infreq(), которая упорядочивает столбцы по частоте (см. советы на странице [Факторы]).

ggplot(
    data = linelist_adm,                       # начинаем с построчного списка, содержащего информацию об административной единице
    mapping = aes(
      x = fct_rev(fct_infreq(admin3name))))+ # ось x - административные единицы, упорядоченные по частоте (в обратном порядке)
  geom_bar()+                                # создаем столбцы, высота - количество строк
  coord_flip()+                              # меняем местами оси X и Y, чтобы легче читать административные единицы
  theme_classic()+                           # упрощаем фон
  labs(                                      # заголовки и подписи
    x = "Admin level 3",
    y = "Number of cases",
    title = "Number of cases, by adminstative unit",
    caption = "As determined by a spatial join, from 1000 randomly sampled cases from linelist"
  )

Ближайшие соседи

Поиск ближайшей медицинской организации/участка прикрепления

Может быть полезным узнать, где расположены медицинские организации относительно очагов заболевания.

Мы можем использовать метод соединения st_nearest_feature из функции st_join() (пакет sf), чтобы визуализировать ближайшую к отдельным случаям медицинскую организацию.

  1. Мы начинаем с построчного списка с шейп-файлом linelist_sf
  2. Мы пространственно соединяем с помощью sle_hf, то есть с расположениями медицинских организаций и клиник (точки)
# Ближайшая медицинская организация к каждому случаю
linelist_sf_hf <- linelist_sf %>%                  # начинаем с шейп-файла с построчным списком  
  st_join(sle_hf, join = st_nearest_feature) %>%   # данные о ближайшей клинике соединяются с данными о случаях 
  select(case_id, osm_id, name, amenity) %>%       # сохраняем интересующие столбцы, включая id, имя, тип, а также геометрию, для медицинских организаций
  rename("nearest_clinic" = "name")                # переименовываем для ясности

Ниже мы можем увидеть (первые 50 строк), что у каждого случая теперь есть данные о ближайшей клинике/больнице

Мы можем увидеть, что клиника “Den Clinic” является ближайшей медицинской организацией примерно для ~30% случаев.

# Подсчет случаев по медицинской организации
hf_catchment <- linelist_sf_hf %>%   # начинаем с построчного списка, включая данные о ближайшей клинике
  as.data.frame() %>%                # конвертируем из шейп-файла в датафрейм
  count(nearest_clinic,              # считаем строки по "имени" (клиники)
        name = "case_n") %>%         # присваиваем новому столбцу с подсчетом имя "case_n"
  arrange(desc(case_n))              # сортируем в порядке уменьшения

hf_catchment                         # выводим на консоль
                         nearest_clinic case_n
1                            Den Clinic    355
2       Shriners Hospitals for Children    316
3         GINER HALL COMMUNITY HOSPITAL    173
4                             panasonic     55
5 Princess Christian Maternity Hospital     38
6                     ARAB EGYPT CLINIC     29
7                                  <NA>     18
8                  MABELL HEALTH CENTER     16

Чтобы визуализировать результаты, мы можем использовать tmap - в этот раз в интерактивном режиме для облегчения просмотра

tmap_mode("view")   # устанавливаем режим tmap на интерактивный  

# строим визуализацию случаев и точек клиник 
tm_shape(linelist_sf_hf) +            # визуализируем случаи
  tm_dots(size=0.08,                  # случаи окрашены по ближайшей клинике
          col='nearest_clinic') +    
tm_shape(sle_hf) +                    # визуализируем медицинские организации большими черными точками
  tm_dots(size=0.3, col='black', alpha = 0.4) +      
  tm_text("name") +                   # накладываем название организации
tm_view(set.view = c(-13.2284, 8.4699, 13), # корректируем масштаб (центр координат, приближение)
        set.zoom.limits = c(13,14))+
tm_layout(title = "Cases, colored by nearest clinic")

Буферные зоны

Мы также можем изучить, сколько случаев находятся на расстоянии 2.5 км (~30 минут) ходьбы от ближайшей медицинской организации.

Примечание: Для более точных расчетов расстояний лучше перепроецировать объект sf в соответствующую местную систему картографических проекций, например UTM (Земля проецируется на плоскую поверхность). В данном примере для простоты мы будем придерживаться географической системы координат Всемирной геодезической системы (WGS84) (Земля представлена в виде сферической/круглой поверхности, поэтому единицы измерения - десятичные градусы). Мы будем использовать общее преобразование: 1 десятичный градус = ~111 км.

См. дополнительную информацию о проекциях карт и системах координат в этой статье esri. В этом блоге рассказывается о разных типах проекций карт и о том, как выбирать подходящую проекцию в зависимости от интересующего района и контекста вашей карты/анализа.

Сначала, создаем круговой буфер с радиусом ~2.5 км вокруг каждой медицинской организации. Это делается с помощью функции st_buffer() из tmap. Поскольку единица карты - десятичные градусы широты/долготы, так интерпретируется “0.02”. Если ваша система координат карты в метрах, надо задавать число в метрах.

sle_hf_2k <- sle_hf %>%
  st_buffer(dist=0.02)       # десятичные градусы, примерно соответствующие 2.5 км 

Ниже мы строим сами буферные зоны следующим образом:

tmap_mode("plot")
# создаем круговые буферы
tm_shape(sle_hf_2k) +
  tm_borders(col = "black", lwd = 2)+
tm_shape(sle_hf) +                    # визуализируем медицинские организации большими красными точками
  tm_dots(size=0.3, col='black')      

**Во-вторых, мы ищем пересечение этих буферов со случаями (точками), используя st_join() с типом соединения st_intersects*. То есть данные о буферах присоединяются к точкам, с которыми они пересекаются.

# пересечение случаев с буферами
linelist_sf_hf_2k <- linelist_sf_hf %>%
  st_join(sle_hf_2k, join = st_intersects, left = TRUE) %>%
  filter(osm_id.x==osm_id.y | is.na(osm_id.y)) %>%
  select(case_id, osm_id.x, nearest_clinic, amenity.x, osm_id.y)

Теперь мы можем посчитать результаты: nrow(linelist_sf_hf_2k[is.na(linelist_sf_hf_2k$osm_id.y),]) из 1000 случаев не пересекаются ни с одним буфером (это значение отсутствует), то есть живут дальше, чем 30 минут ходьбы от ближайшей медицинской организации.

# Случаи, у которых нет пересечения с буферами медицинских организаций
linelist_sf_hf_2k %>% 
  filter(is.na(osm_id.y)) %>%
  nrow()
[1] 1000

Мы можем визуализировать результаты таким образом, чтобы случаи, не пересекающиеся с буфером, отражались красным.

tmap_mode("view")

# Сначала отобразим случаи точками
tm_shape(linelist_sf_hf) +
  tm_dots(size=0.08, col='nearest_clinic') +

# визуализируем медицинские организации большими черными точками
tm_shape(sle_hf) +                    
  tm_dots(size=0.3, col='black')+   

# затем наложим буферы медицинских организаций полилиниями
tm_shape(sle_hf_2k) +
  tm_borders(col = "black", lwd = 2) +

  # выделим случаи, которые не попали в буферы ни одной медицинской организации
# красными точками  
tm_shape(linelist_sf_hf_2k %>%  filter(is.na(osm_id.y))) +
  tm_dots(size=0.1, col='red') +
tm_view(set.view = c(-13.2284,8.4699, 13), set.zoom.limits = c(13,14))+

# добавим заголовок  
tm_layout(title = "Cases by clinic catchment area")