Создание границы кластера с использованием подхода K-ближайших соседей
Пару месяцев назад я написал здесь на Medium статью о картировании горячих точек дорожно-транспортных происшествий в Великобритании. Больше всего меня интересовала иллюстрация использования алгоритма кластеризации DBSCAN на географических данных. В статье я использовал географическую информацию о зарегистрированных дорожно-транспортных происшествиях, опубликованную правительством Великобритании. Моей целью было запустить процесс кластеризации на основе плотности, чтобы найти районы, в которых чаще всего регистрируются дорожно-транспортные происшествия. Конечным результатом стало создание набора геозон, представляющих эти горячие точки аварий.
Собирая все точки в данном кластере, вы можете получить представление о том, как кластер будет выглядеть на карте, но вам будет не хватать важной информации: внешней формы кластера. В данном случае речь идет о замкнутом многоугольнике, который можно представить на карте как геозону. Можно предположить, что любая точка внутри геозоны принадлежит кластеру, что делает эту форму интересной информацией: вы можете использовать ее как функцию дискриминатора. Все новые точки выборки, которые попадают внутрь многоугольника, можно считать принадлежащими соответствующему кластеру. Как я намекал в статье, вы можете использовать такие многоугольники, чтобы подтвердить свой риск при вождении, используя их для классификации вашего собственного местоположения по GPS.
От облака к многоугольнику
Теперь вопрос в том, как создать значимый многоугольник из облака точек, составляющих определенный кластер. Мой подход в первой статье был несколько наивным и отражал решение, которое я уже использовал в производственном коде. Это решение повлекло за собой размещение круга с центром в каждой из точек кластера, а затем объединение всех кругов вместе, чтобы сформировать многоугольник в форме облака. Результат не очень красивый и не реалистичный. Кроме того, при использовании кругов в качестве базовой формы для построения конечного многоугольника у них будет больше точек, чем у более обтекаемой формы, что приведет к увеличению затрат на хранение и замедлению работы алгоритмов обнаружения включения.
С другой стороны, этот подход имеет преимущество простоты вычислений (по крайней мере, с точки зрения разработчика), поскольку он использует функцию cascaded_union Shapely для объединения всех кругов вместе. Еще одно преимущество состоит в том, что форма многоугольника неявно определяется с использованием всех точек в кластере.
Для более сложных подходов нам нужно каким-то образом идентифицировать граничные точки кластера, те, которые, кажется, определяют форму облака точек. Интересно, что с некоторыми реализациями DBSCAN [1] вы можете фактически восстановить пограничные точки как побочный продукт процесса кластеризации. К сожалению, эта информация (по всей видимости) недоступна для реализации SciKit Learn [2], поэтому мы должны действовать.
Первый подход, который пришел в голову, заключался в вычислении выпуклой оболочки множества точек. Это хорошо понятный алгоритм, но он страдает от проблемы с обработкой вогнутых форм, таких как этот:
Эта форма неправильно отражает суть лежащих в основе точек. Если бы он использовался в качестве дискриминатора, некоторые точки были бы неправильно классифицированы как находящиеся внутри кластера, когда это не так. Нам нужен другой подход.
Альтернатива вогнутого корпуса
К счастью, этому положению дел есть альтернативы: мы можем рассчитать вогнутую оболочку. Вот как выглядит вогнутый корпус при применении к тому же набору точек, что и на предыдущем изображении:
Или, может быть, вот этот:
Как видите, в отличие от выпуклой оболочки, нет единого определения того, что такое вогнутая оболочка набора точек. В алгоритме, который я представляю здесь, выбор того, насколько вогнутой должны быть ваши корпуса, осуществляется с помощью одного параметра: k - количество ближайших соседей, учитываемых при расчете корпуса. Посмотрим, как это работает.
Алгоритм
Алгоритм, который я представляю здесь, был описан более десяти лет назад Адриано Морейра и Марибель Ясмина Сантос из Университета Минью, Португалия [3]. Из аннотации:
В этой статье описывается алгоритм вычисления огибающей набора точек на плоскости, который генерирует выпуклые на невыпуклых оболочках, которые представляют площадь, занимаемую данными точками. Предлагаемый алгоритм основан на подходе k-ближайших соседей, где значение k, единственного параметра алгоритма, используется для управления «гладкостью» окончательного решения. […]
Поскольку я буду применять этот алгоритм к географической информации, пришлось внести некоторые изменения, а именно при вычислении углов и расстояний [4]. Но это никоим образом не меняет сути алгоритма, который можно в общих чертах описать следующими шагами:
- Найдите точку с наименьшей координатой y (широта) и сделайте ее текущей.
- Найдите k -ближайшие точки к текущей точке.
- Из k -ближайших точек выберите ту, которая соответствует наибольшему правому повороту от предыдущего угла. Здесь мы воспользуемся понятием пеленга и начнем с угла 270 градусов (из-за West).
- Проверьте, не пересекается ли, добавляя новую точку в растущую строку линии, саму себя. Если это так, выберите другую точку из k -ближайшей или перезапустите с большим значением k.
- Сделайте новую точку текущей и удалите ее из списка.
- После k итераций добавьте первую точку обратно в список.
- Петля к номеру 2.
Алгоритм кажется довольно простым, но есть ряд деталей, на которые необходимо обратить внимание, особенно потому, что мы имеем дело с географическими координатами. Расстояния и углы измеряются по-другому.
Код
Здесь я публикую адаптированную версию кода предыдущей статьи. Вы по-прежнему найдете тот же код кластеризации и тот же генератор кластеров в форме облака. Обновленная версия теперь содержит пакет с именем geomath.hulls, в котором вы можете найти класс ConcaveHull. Чтобы создать вогнутый корпус, сделайте следующее:
В приведенном выше коде points - это массив измерений (N, 2), где строки содержат наблюдаемые точки, а столбцы - географические координаты (долгота, широта). . Результирующий массив имеет точно такую же структуру, но содержит только точки, принадлежащие многоугольной форме кластера. Своего рода фильтр.
Поскольку мы будем обрабатывать массивы, вполне естественно ввести NumPy в бой. Все вычисления были должным образом векторизованы, когда это было возможно, и были предприняты усилия для повышения производительности при добавлении и удалении элементов из массивов (спойлер: они вообще не перемещаются). Одно из недостающих улучшений - распараллеливание кода. Но это может подождать.
Я организовал код вокруг алгоритма, как показано в статье, хотя во время перевода были сделаны некоторые оптимизации. Алгоритм построен на нескольких подпрограммах, которые четко обозначены в документе, так что давайте прямо сейчас избавимся от них. Для вашего удобства я буду использовать те же имена, что и в статье.
CleanList [listOfPoints] - очистка списка точек выполняется в конструкторе класса:
Как видите, список точек реализован в виде массива NumPy из соображений производительности. Очистка списка выполняется в строке 10, где сохраняются только уникальные точки. Массив набора данных организован с наблюдениями в строках и географическими координатами в двух столбцах. Обратите внимание, что я также создаю логический массив в строке 13, который будет использоваться для индексации в массиве основного набора данных, облегчая рутинную работу по удалению элементов и, время от времени, добавлению их обратно. Я видел эту технику, называемую «маской», в документации NumPy, и она очень эффективна. Что касается простых чисел, я расскажу о них позже.
FindMinYPoint [listOfPoints]. Для этого требуется небольшая функция:
Эта функция вызывается с массивом набора данных в качестве аргумента и возвращает индекс точки с наименьшей широтой. Обратите внимание, что строки кодируются долготой в первом столбце и широтой во втором.
RemovePoint [вектор, e]
AddPoint [вектор, e] - это несложно благодаря использованию массива indices. Этот массив используется для хранения активных индексов в основном массиве набора данных, поэтому удаление элементов из набора данных выполняется очень просто.
Хотя алгоритм, описанный в документе, требует добавления точки к массиву, составляющему корпус, на самом деле это реализовано как:
Позже переменная test_hull будет снова присвоена hull, когда строка будет считаться непересекающейся. Но здесь я забегаю вперед. Удалить точку из массива набора данных так же просто, как:
self.indices[current_point] = False
Чтобы добавить его обратно, нужно просто вернуть значение этого массива с тем же индексом к истине. Но за все это удобство приходится расплачиваться за то, что нам приходится следить за индексами. Подробнее об этом позже.
NearestPoints [listOfPoints, point, k] - Здесь начинает становиться интереснее, потому что мы не имеем дело с плоскими координатами, поэтому перестаньте использовать Pythagoras и введите Haversine:
Обратите внимание, что второй и третий параметры представляют собой массивы в формате набора данных: долгота в первом столбце и широта во втором. Как видите, эта функция возвращает массив расстояний в метрах между точкой во втором аргументе и точками в третьем аргументе. Когда они у нас есть, мы можем легко получить k -ближайших соседей. Но для этого есть специальная функция, и она заслуживает пояснений:
Функция начинается с создания массива с базовыми индексами. Это индексы точек, которые не были удалены из массива набора данных. Например, если в кластере из десяти точек мы начали с удаления первой точки, массив базовых индексов будет [1, 2, 3, 4, 5, 6, 7, 8, 9]. Далее мы вычисляем расстояния и сортируем полученные индексы массива. Первые k извлекаются и затем используются в качестве маски для получения базовых индексов. Это немного искажено, но работает. Как видите, функция возвращает не массив координат, а массив индексов в массив набора данных.
SortByAngle [listOfPoints, point, angle] - здесь больше проблем, потому что мы не вычисляем простые углы, мы вычисляем пеленги. Они измеряются как ноль градусов по северу с углами, увеличивающимися по часовой стрелке. Вот ядро кода, который вычисляет подшипники:
Функция возвращает массив пеленгов, измеренных от точки, индекс которой находится в первом аргументе, до точек, индексы которых находятся в третьем аргументе. Сортировка проста:
На этом этапе массив кандидатов содержит индексы ближайших k точек, отсортированные по убыванию пеленга.
IntersectQ [lineSegment1, lineSegment2] - вместо того, чтобы использовать собственные функции пересечения линий, я обратился за помощью к Shapely. Фактически, при построении многоугольника мы, по сути, обрабатываем линейную строку, добавляя сегменты, которые не пересекаются с предыдущими. Проверить это просто: мы берем строящийся массив корпуса, конвертируем его в строковый объект Shapely line и проверяем, является ли он простым (не самопересекающимся) или нет.
Короче говоря, строка Shapely line становится сложной, если она самопересекается, поэтому предикат is_simple становится ложным. Легкий.
PointInPolygon [point, listOfPoints]. Это оказалось самым сложным в реализации. Позвольте мне объяснить, посмотрев на код, который выполняет окончательную проверку многоугольника корпуса (проверяет, все ли точки кластера включены в многоугольник):
Функций Shapely для проверки пересечения и включения должно было быть достаточно, чтобы проверить, перекрывает ли окончательный многоугольник корпуса все точки кластера, но это не так. Почему? Shapely не зависит от координат, поэтому он будет обрабатывать географические координаты, выраженные в широте и долготе, точно так же, как координаты на декартовой плоскости. Но мир ведет себя по-другому, когда вы живете на сфере, и углы (или направления) на геодезической не постоянны. Пример в ссылке [4] геодезической линии, соединяющей Багдад и Осака, прекрасно иллюстрирует это. Бывает так, что при некоторых обстоятельствах алгоритм может включать точку, основанную на критерии пеленга, но позже, используя планарные алгоритмы Shapely, будет считаться, что она немного выходит за пределы многоугольника. Вот что здесь делает коррекция на малое расстояние.
Мне потребовалось время, чтобы понять это. Моим помощником по отладке была QGIS, отличное бесплатное программное обеспечение. На каждом этапе подозрительных вычислений я выводил данные в формате WKT в файл CSV для считывания в виде слоя. Настоящая палочка-выручалочка!
Наконец, если многоугольник не покрывает все точки кластера, единственный вариант - увеличить k и повторить попытку. Здесь я добавил немного собственной интуиции.
Prime k
В статье предлагается увеличить значение k на единицу и выполнить алгоритм заново с нуля. Мои ранние тесты с этой опцией были не очень удовлетворительными: время выполнения на проблемных кластерах было бы медленным. Это произошло из-за медленного увеличения k, поэтому я решил использовать другой график увеличения: таблицу простых чисел. Алгоритм уже начинается с k = 3, поэтому его можно было легко расширить, чтобы он развился на основе списка простых чисел. Вот что вы видите в рекурсивном вызове:
Знаешь, я люблю простые числа ...
Взрывать
Вогнутые многоугольники корпуса, сгенерированные этим алгоритмом, все еще нуждаются в некоторой дальнейшей обработке, потому что они будут различать только точки внутри корпуса, но не близко к нему. Решение состоит в том, чтобы добавить к этим узким кластерам немного отступов. Здесь я использую ту же технику, что и раньше, и вот как это выглядит:
Здесь я использовал buffer функцию Shapely, чтобы добиться цели.
Функция принимает Shapely Polygon и возвращает его расширенную версию. Второй параметр - это радиус в метрах добавляемой прокладки.
Запуск кода
Начните с загрузки кода из репозитория GitHub на свой локальный компьютер. Файл, который вы хотите запустить, - это ShowHotSpots.py в основном каталоге. При первом запуске код считывает данные о дорожно-транспортных происшествиях в Великобритании с 2013 по 2016 год и группирует их. Затем результаты кэшируются в виде файла CSV для последующих запусков.
Затем вам будут представлены две карты: первая генерируется с использованием облачных кластеров, а вторая использует алгоритм вогнутой кластеризации, обсуждаемый здесь. Во время выполнения кода генерации многоугольника вы можете увидеть, что сообщается о нескольких сбоях. Чтобы понять, почему алгоритм не может создать вогнутую оболочку, код записывает кластеры в файлы CSV в каталог data/out/failed/. Как обычно, вы можете использовать QGIS для импорта этих файлов в виде слоев.
По сути, этот алгоритм не работает, когда он не находит достаточно точек, чтобы «обойти» форму без самопересечения. Это означает, что вы должны быть готовы либо отбросить эти кластеры, либо применить к ним другую обработку (выпуклая оболочка или сросшиеся пузыри).
Concaveness
Это обертка. В этой статье я представил метод постобработки сгенерированных DBSCAN географических кластеров в вогнутые формы. Этот метод может обеспечить лучшее соответствие внешнего многоугольника кластерам по сравнению с другими альтернативами.
Спасибо за чтение и получайте удовольствие от работы с кодом!
использованная литература
[1] Крышкевич М., Ласек П. (2010) TI-DBSCAN: Кластеризация с помощью DBSCAN с помощью неравенства треугольника. В: Szczuka M., Kryszkiewicz M., Ramanna S., Jensen R., Hu Q. (eds) Rough Sets and Current Trends in Computing. RSCTC 2010. Lecture Notes in Computer Science, vol 6086. Springer, Berlin, Heidelberg [Springer Link]
[2] Scikit-learn: машинное обучение в Python, Педрегоса и др., JMLR 12, стр. 2825–2830, 2011 г.
[3] Морейра А. и Сантос М.Ю., 2007, Вогнутая оболочка: подход с использованием K-ближайших соседей для вычисления области, занимаемой набором точек [PDF]
[4] Вычислить расстояние, азимут и многое другое между точками широты и долготы
