Как геокодировать миллион точек на Spark по-быстрому?

    В моем предыдущем проекте перед нами встала задача провести обратное геокодирование для множества пар географических координат. Обратное геокодирование — это процедура, которая паре широта-долгота ставит в соответствие адрес или название объекта на карте, к которому принадлежит или близка заданная координатами точка. То есть, берем координаты, скажем такие: @55.7602485,37.6170409, и получаем результат либо «Россия, Центральный федеральный округ, Москва, Театральная площадь, дом такой-то», либо например «Большой театр».

    Если на входе адрес или название, а на выходе координаты, то эта операция — прямое геокодирование, об этом мы, надеюсь, поговорим позже.

    В качестве исходных данных у нас на входе было примерно 100 или 200 тысяч точек, которые лежали в кластере Hadoop в виде таблицы Hive. Это чтобы был понятен масштаб задачи.

    В качестве инструмента обработки в конце концов был выбран Spark, хотя в процессе мы попробовали как MapReduce, так и Apache Crunch. Но это отдельная история, возможно заслуживающая своего поста.

    Прямолинейное решение доступными средствами


    Для начала мы попробовали подойти к проблеме так сказать «в-лоб». В качестве инструмента имелся сервер ArcGIS, который предоставляет REST-сервис обратного геокодирования. Пользоваться им достаточно просто, для этого выполняем http GET запрос со следующим URL:

    http://базовый-url/GeocodeServer/reverseGeocode?<параметры>
    

    Из множества параметров достаточно задать location=x,y (главное не перепутайте, кто из них широта, а кто — долгота ;). И вот мы уже имеем JSON с результатами: страна, регион, город, улица, номер дома. Пример из документации:

    {
     "address": {
      "Match_addr": "Beeman's Redlands Pharmacy",
      "LongLabel": "Beeman's Redlands Pharmacy, 255 Terracina Blvd, Redlands, CA, 92373, USA",
      "ShortLabel": "Beeman's Redlands Pharmacy",
      "Addr_type": "POI",
      "Type": "Pharmacy",
      "PlaceName": "Beeman's Redlands Pharmacy",
      "AddNum": "255",
      "Address": "255 Terracina Blvd",
      "Block": "",
      "Sector": "",
      "Neighborhood": "South Redlands",
      "District": "",
      "City": "Redlands",
      "MetroArea": "Inland Empire",
      "Subregion": "San Bernardino County",
      "Region": "California",
      "Territory": "",
      "Postal": "92373",
      "PostalExt": "",
      "CountryCode": "USA"
     },
     "location": {
      "x": -117.20558993392585,
      "y": 34.037880040538894,
      "spatialReference": {
       "wkid": 4326,
       "latestWkid": 4326
      }
     }
    }
    

    Можно дополнительно указать, какие виды ответов мы хотим — почтовые адреса, или скажем POI (point of interest, это чтобы получить ответы вида «Большой театр»), или нужны ли нам пересечения улиц, например. Также можно указать радиус, внутри которого от указанной точки будет произведен поиск именованного объекта.

    Для быстрой проверки качества ответа можно прикинуть расстояние между исходной точкой в параметрах запроса, и точкой полученной — location в ответе сервиса.

    Казалось бы, сейчас все будет хорошо. Но не тут-то было. Наш экземпляр ArcGIS был достаточно медленный, серверу было выделено кажется 4 ядра, и примерно 8 гигабайт оперативной памяти. В итоге задача на кластере могла прочитать наши 200 тыс точек очень быстро, но упиралась в REST и производительность ArcGIS. И геокодирование всех точек занимало часы. При этом на Hadoop мы выделяли всего каких-то 8 ядер, и немного памяти, но так как загрузка сервера ArcGIS была близка к 100% в течение многих часов, лишние ресурсы в кластере нам ничего не давали.

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

    Второе приближение, внедряем кеш


    Для начала мы выяснили, что многие точки в нашем наборе имеют повторяющиеся координаты. Причина тут простая — очевидно, точность GPS недостаточно хороша, чтобы координаты расположенных двух метрах друг от друга точек отличались на выходе, или же в исходную базу просто вводились координаты, полученные не с GPS, а из другой базы. В общем, неважно, почему это так, главное что это вполне типичная ситуация, так что кеш результатов, полученных от сервиса, позволит выполнить геокодирование каждой пары координат только однократно. А памяти для кеша мы вполне можем себе позволить.

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

    Таким несложным образом мы смогли получить ускорение примерно до 10 раз, что примерно на глаз соответствует числу повторов координат в исходном наборе. Это уже было приемлемо, но все равно очень медленно.

    Хорошо, сказал нам в это время наш заказчик, если мы не можем вычислить адреса быстрее, можем ли мы быстро определить хотя бы город? И мы взялись за…

    Упрощенное решение, внедряем Geomerty API


    Что у нас есть, чтобы определить город? У нас была геометрия регионов России, административно-территориальное деление, примерно с точностью до районов города.

    Взять эти данные можно например вот тут. Что там лежит? Это база данных административных границ РФ, для уровней от 2 (страна) до 9 (городские округа). Формат либо geojson, либо CSV (при этом сама геометрия в формате wkt). Всего в базе примерно 20 тыс записей.

    Новое упрощенное решение задачи выглядело так:

    1. Загружаем данные АТД в Hive.
    2. Для каждой точки с координатами ищем в таблице территориального деления полигоны, куда эта точка входит.
    3. Найденные полигоны сортируем по уровню.

    В итоге на выходе получаем что-то вроде: Россия, Центральный федеральный округ, Москва, такой-то административный округ, район такой-то, то есть, список территорий, в которые входит наша точка.

    Загрузка АТД


    Чтобы загрузить CSV попроще, мы применим kite. Этот инструмент умеет очень хорошо строить схему для Hive на основе заголовков колонок в CSV. Фактически, импорт сводится к трем командам (одна из которых повторяется для каждого файла уровня):

    kite-tools csv-schema admin_level_2.csv --class al --delimiter \; >adminLevel.avrs
    kite-tools create dataset:hive:/default/levelswkt -s adminLevel.avrs
    kite-tools csv-import admin_level_2.csv dataset:hive:/default/levelswkt --delimiter \;
    ...
    kite-tools csv-import admin_level_10.csv dataset:hive:/default/levelswkt --delimiter \;
    

    Что мы тут сделали? Первая команда построила нам Avro-схему для csv, для чего мы указали некоторые параметры схемы (класс, в данном случае), и разделитель полей для CSV. Далее стоит посмотреть на полученную схему глазами, и возможно внести некоторые исправления, потому что Kite при анализе данных смотрит не на все строки нашего файла, а лишь на некоторую выборку, поэтому может делать иногда неверные предположения о типе данных (увидел три числа — решил что колонка числовая, а дальше пошли строки).

    Ну а дальше на основе схемы мы создаем dataset (это общий термин Kite, обобщающий таблицы в Hive, таблицы в HBase, и кое-то еще). В данном случае default это база данных (для Hive тоже что и схема), а levelswkt — название нашей таблицы.

    Ну и последние команды загружают CSV файлы в наш датасет. После успешного завершения загрузки вы уже можете выполнить запрос:

    select * from levelswkt;

    где-нибудь в Hue.

    Работа с геометрией


    Для работы с геометрией в Java мы выбрали Java Geometry API от компании Esri (разработчик ArcGIS). В принципе, можно было взять и другие фреймворки, некоторый выбор Open Source имеется, например широко известный пакет JTS Topology Suite, или Geotools.

    С первой задачей нам позволяет справиться еще один фреймворк от той же компании Esri, именуемый Spatial Framework for Hadoop, и основанный на первом. По сути, из него нам нужен так называемый SerDe, модуль сериализации-десериализации для Hive, который позволяет представить пачку geojson файлов как таблицу в Hive, колонки которой берутся из атрибутов geojson. А сама геометрия становится еще одной колонкой (с бинарными данными). В итоге у нас есть таблица, одна из колонок которой — геометрия некоего региона, а остальные — его атрибуты (название, уровень в АТД, и пр.). Эта таблица доступна Spark приложению.

    Если же мы загружаем в базу данные в формате CSV, то мы имеем колонку, где геометрии лежат в текстовом виде, в формате WKT. В этом случае Spark может создать объект геометрии при выполнении, применяя Geometry API.

    Мы выбрали именно CSV формат (и WKT), по одной простой причине — как все знают, Россия занимает на карте область с координатами Чукотки за 180 меридианом. У формата geojson есть ограничение — все полигоны внутри него должны быть ограничены 180 градусами, а те что пересекают 180 меридиан, нужно порезать на две части по нему. В итоге, при импорте геометрии в Geometry API мы получаем объект, для которого Geometry API неправильно определяет Bounding Box (объемлющий прямоугольник) для границы России. Получается ответ -180,180 по долготе. Что конечно же не правда — в реальности Россия занимает примерно от 20 до -170 градусов. Это проблема Geometry API, на сегодня она возможно уже исправлена, но тогда нам пришлось ее обходить.

    У WKT такой проблемы нет. Вы спросите, зачем нам Bounding Box? Потом расскажу ;)

    Остается решить так называемую задачу PIP, point in poligon. Это снова умеет Java Geometry API, для нас это просто, одна геометрия типа Point, вторая полигон (Multipoligon) для региона, и один метод contains.

    В итоге, второе решение, и тоже в лоб, выглядело так: Spark приложение загружает АТД, включая геометрии. Из них строится нечто типа Map название->геометрия (на самом деле чуть сложнее, так как АТД вложены друг в друга, и нам нужно искать только в тех нижних уровнях, которые входят в уже найденный верхний. В итоге тут некое дерево геометрий, которое по исходным данным еще нужно построить). А дальше мы строим Spark Dataset с нашими точками, и к каждой точке применяем свою UDF, которая проверяет вхождение точки во все геометрии (по дереву).

    Написание новой версии занимает примерно день работы, благо в комплекте Spatial Framework for Hadoop были неплохие примеры прямо для решения PIP задачи (правда, несколько другими средствами).

    Запускаем, и… о ужас, быстрее-то не стало. Снова часы. Пора подумать над оптимизацией.

    Оптимизированное решение, QuadTree


    Причина тормозов достаточно очевидная — скажем, геометрия России, т.е. внешние границы, это мегабайты geojson, здоровенный полигон, и не один. Если вспомнить, как решается задачка PIP, то одним из широко известных методов является построение луча из точки скажем куда-то вверх, до бесконечности, и определение, в скольких точках он пересекает полигон. Если число точек четное, точка вне полигона, если нечетное — внутри.

    Вот скажем описание из вики.



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

    Для этого нам подойдет дерево квадрантов. Благо реализация есть все в том же Geometry API (и много где еще).



    Решение на основе дерева выглядит примерно так:

    1. Загружаем геометрии АТД
    2. Для каждой геометрии определяем объемлющий прямоугольник
    3. Помещаем ее в QuadTree, получаем в ответ индекс
    4. индекс запоминаем

    Далее, при обработке точек:

    1. Спрашиваем у QuadTree, в какие из известных ей геометрий может входить точка
    2. Получаем индексы геометрий
    3. Только для них проверяем вхождение решением задачи PIP

    На это все уходит еще часа четыре разработки. Запускаем еще раз, и видим, что задачка завершилась уж как-то очень быстро. Проверяем — все нормально, решение получено. И все — за пару минут. QuadTree дает нам ускорение на несколько порядков.

    Итоги


    Итак, что мы в итоге получили? Мы получили механизм обратного геокодирования, который эффективно параллелится на Hadoop кластере, и который решает нашу исходную задачу с 200 тыс точек примерно за минуты. Т.е. мы спокойно можем применить это решение к миллионам точек.

    Какие ограничения у данного решения? Первое, очевидное — оно основано на доступных нам данных АТД, которые а) могут быть не актуальны б) ограничены только Россией.

    И второе — мы не умеем решать задачу обратного геокодирования для объектов, отличных от замкнутых полигонов. А это значит — для улиц в том числе.

    Развитие


    Что с этим можно сделать?

    Чтобы иметь актуальные геометрии АТД, самое простое — взять их из OpenStreetMap. Над ними конечно придется немного поработать, но это вполне решаемая задача. Если будет интерес, о задаче загрузки данных OpenStreetMap в кластер Hadoop я расскажу как-нибудь в другой раз.

    Что можно сделать для улиц и домов? В принципе, улицы есть в том же OSM. Но это уже не замкнутые структуры, а линии. Чтобы определить, что точка «близка» к некоторой улице, придется построить для улицы полигон из равноудаленных от нее точек, и проверить попадание в него. В итоге эдакая сосиска получается… выглядит это примерно так:



    Насколько точка близка? Это является параметром, и примерно соответствует тому радиусу, внутри которого ищет объекты ArcGIS, и о котором я упоминал в самом начале.

    Таким образом мы находим улицы, расстояние от точки до которых меньше некоторого предела (скажем, 100 метров). И чем этот предел меньше — тем быстрее работает алгоритм, но тем больше шансов, что вы не найдете ни одного совпадения.

    Очевидная проблема состоит в том, что вычислить эти так называемые buffers заранее невозможно — их размер является параметром сервиса. Их нужно строить на лету, после того, как мы определили искомый район города, и отобрали из базы OSM улицы, которые проходят через этот район. Улицы, впрочем, можно отобрать и заранее.

    Дома, которые находятся в найденном районе, тоже никуда не движутся, поэтому и их список можно построить заранее — но попадание в них все равно придется считать налету.

    То есть, нужно предварительно построить индекс вида «район города»-> список домов со ссылками на геометрии, и аналогичный для списка улиц.

    И как только мы определили район, достаем список домов и улиц, строим по улицам границы, и уже только для них решаем задачу PIP (вероятно, с применением тех же оптимизаций, как и для границ регионов). Дерево квадрантов для домов в этом случае очевидно можно тоже построить заранее, и куда-то сохранить.

    Основная наша цель — максимально сократить объем вычислений, при этом максимально рассчитав и сохранив заранее все, что можно рассчитать и сохранить. В этом случае процесс будет состоять из медленного этапа построения индексов, и второго этапа расчета, который уже будет проходить быстро, близко к онлайн-варианту.
    Share post

    Similar posts

    AdBlock has stolen the banner, but banners are not teeth — they will be back

    More
    Ads

    Comments 23

    • UFO just landed and posted this here
        0
        Про 200 тыс в основном. То что они лежали в Hive, определяло только выбор инструмента в целом.
        0
        А использование одного из платных сервисов не рассматривалось? Вижу вот у гугла 5$ за 1000, возможно есть и подешевле
          0
          Во-первых, за 200000 точек это будет уже $10000, во-вторых, это скорее всего будет медленный REST.

          Но самое главное — у нас интранет, где нет никакого гугля :)
          +1

          ArcGIS судя по всему просто использует модуль PostGIS, т.е. можно было бы просто загнать данные карт в локальный постгрес и потом выбирать напрямую из него по координатам с помощью ST_Intersects, ST_CoveredBy, ST_Within и тд. Эта схема относительно легко масштабируется.

            0
            ArcGIS в нашем случае живет на Oracle, так что использовать PostGIS точно не может :)

            Но вообще вариант взять координаты из базы ArcGIS рассматривался, и был отброшен, потому что наша база не обновляется — денег не заплатили.
              0
              Кстати, Spatial Framework for Hadoop включает в себя все те же функции, ST_Contains например, потому что их набор — стандарт OGC. Хотя да, spatial индексы, как PostGIS, оно само для нас не построит, поэтому если тупо применить ST_Contains — скорее всего будет медленно.
                0

                Не совсем понял, что мешает залить в базу тот же OpenStreetMap, придется конечно немного поразбиратся, куда без этого. Сейчас многие БД предоставляют возможности работать с гео данными, даже SQLite можно поднять в памяти с соотв. расширением. Вопрос сводится к следующему, будет ли реализация реверс геокодинга на хадупе/спарке, быстрее/проще/лучше, чем аналогичная (и уже готовая) реализация в одной из таких БД? Причем сама БД статическая, скопировал на вторую машину запустил и получил гарантированный х2 перформанс.

                  0
                  >чем аналогичная (и уже готовая) реализация в одной из таких БД?
                  Сможете назвать парочку готовых?

                  >Причем сама БД статическая
                  Вообще-то нет. База адресов и POI вовсе не статическая, и процесс ее периодического обновления далеко не такая простая вещь, как может показаться. Особенно если вы берете данные из OSM — в этом случае очень желательно свои данные туда же и возвращать.

                  >быстрее/проще/лучше, чем аналогичная
                  Описывал я это все вовсе не как готовое решение, которое проще или лучше (хотя оно и сильно быстрее многих), а как подход, позволяющий быстро реализовать такое самостоятельно. Цели сделать лучше не ставилось, и далеко не очевидно, что именно тут считать лучшим.

                  Развертывание нескольких машин с Postgres вполне возможно, но мне не кажется решением без недостатков. Основной из которых в том, что если задачи геокодирования это лишь малая часть задач кластера (а у нас именно так), то отдельный Postgres — это отдельная немаленькая морока по администрированию. Он либо работает постоянно, потребляя ресурсы когда не нужен, либо надо как-то организовать его запуск по требованию. А так да, вполне себе решение. У нас были планы развертывания ArcGIS с его базой внутри кластера, для этих целей — и в частности по этим же причинам они не реализованы (там еще и лицензии примешиваются).

                  >Не совсем понял, что мешает залить в базу тот же OpenStreetMap, придется конечно немного поразбиратся, куда без этого.
                  Ничего не мешает. Про это — следующий пост, который уже ждет своей очереди.
              +1
              Вот еще отличный вариант с использованием Lucene — opensourceconnections.com/blog/2014/04/11/indexing-polygons-in-lucene-with-accuracy
                0
                Lucene — это не спортивно :)

                В смысле, я пока рассказал идею, как это сделать вручную самому, в 400 строк своего кода на Spark. Про то, что есть другие инструменты, в том числе пригодные для применения в кластере — возможно расскажу в следующий постах.
                  +2
                  А в чем спорт?

                  С Lucene заранее например индексируем точки городов мира отсюда, а потом загружаем на каждом Executor в Spark и опрашиваем.

                  Поиск в радиусе на индексе из ~3 млн точек и точности в ~1 метр берет окоро 2 мс на запрос. А вот поиск топ ближайших городов работает значительно медленнее 700-1000 мс на запрос, т.к. индекс по всей видимости не используется.
                    0
                    >А в чем спорт?
                    В том, чтобы сделать обратный геокодер в виде примерно 400 строк кода. Я никого не убеждаю, что надо делать именно так, это просто был пример проблем, которые приносит с собой масштаб, и как их можно быстро решить.

                    Что до точек — то в моем понимании с точками все сильно проще. Расчет расстояния между точками — это константное время, ну может чуть-чуть сложнее, чем на плоскости. Уже для улиц это не так. Ваш пример вполне применим для домов, например, если считать их точками (в реальности в OSM далеко не все дома точки, но построить некий центр для дома вполне можно заранее).
                  0
                  Кстати, хорошая статья.

                  Мы реально смотрели в сторону Solr, но скорее для задачи прямого геокодирования, там это нужнее.
                    0
                    Ага, статья грамотная. Я собираюсь с автором кстати пообщаться в ближайшее время. Есть у меня к нему ряд вопросов.
                  +1
                  А геоспарк не пробовали? Вполне неплохо себя показывает. Если индекс а-ля геохэш запилить, то и вообще отлично. Мы вот юзаем, не жалуемся. У нас порядка миллионов точек на одну карту, карты считаются на потоке.
                    0
                    Я тут описал историю в том виде, как она развивалась в жизни. На тот момент — нет, не пробовали. Потом уже дозрели, что нужно внедрять что-то готовое. Смотрели на GeoSpark, и еще на GeoMesa и семейство продуктов от тех же авторов, которые в тоже время и авторы JTS.

                    Собственно, главный вывод из этой истории состоял в том, что нужно попытаться построить spatial индекс в том или ином виде, а иначе будет плохо. Ну и как бы GeoSpark это умеет. В целом же нам GeoMesa чуть больше понравилась, хотя порог входа наверное чуть выше.

                    И да, геохеш на самом деле я зря не упомянул, потому что это по сути одна из альтернатив quad tree.
                      0
                      Ну мы как-то сразу с MR, который дёргает st_чегоНибудь() из PostgreSQL с PostGIS, перешли на геоспарк с индексами по геохешу, без промежуточных шагов. Меньше чем за полгода техпроцесс перевели, а итоговое ускорение расчёта получилось в пару сотен раз. На изначальном стеке расчёт карты по десятку тысяч точек был тем ещё предприятием, теперь бывает до сотен миллионов обсчитываем в приемлемое время на кластере из пяти нод среднего размера.

                      Довольно здорово дела упрощает ещё алгоритмическая сетка — заказчик был из Токио, там у них есть такая Japan Mesh. Она иерархическая, зависит только от координат, и её можно использовать сразу как индекс. В будущем планируем заюзать что-то такое же, только для всего мира.
                        0
                        Дело в том, что собственно геометрии у нас в основном вторичны. А большая часть логики завязана на адреса, а адреса — это уже совсем не ST_..., а совсем другая логика.
                          0
                          Почему другая? Адрес — это разве не метаданные аутлайна?
                            0
                            Ну если коротко — то нет. Ну или вернее, там все сильно сложнее. Во-первых, OSM не хранит в тэгах (по большей части) все компоненты адреса, так что для определения скажем города вам придется решать геометрическую задачу, и только дом и корпус к примеру прописан тэгом (причем далеко не очевидно, на чем именно он будет стоять). И это рекомендованное поведение — почему-то кое-кто считает, что геометрические отношения для этого применять правильнее.

                            А во-вторых, что важнее, там слишком много вариаций записи, и если обратное геокодирование делается достаточно примитивно, то прямое — это уже ужас-ужас.
                              0
                              Хм. Ну, в наших задачах чаще попадается либо сопоставление с postcode/zip/почтовым индексом (это всегда мозаика из аутлайнов конкретных участков, вплоть до 1 дома), либо в сильно ограниченной области, где не только город, но и район города заранее точно известен.

                              Или дата аналитик как-то сам подготавливает портянку аутлайнов в geojson с метаданными по адресам, или ячейками грида, или чем-то другим, что нужно для текущей карты. Откуда он их берёт, кроме OSM, я, если честно, не интересовался.

                              Ваша задача звучит посложнее.
                                0
                                >Откуда он их берёт, кроме OSM, я, если честно, не интересовался
                                Ну, у нас есть похожий проект в другом департаменте, связанный с оценкой стоимости недвижимости, если его совсем просто описать. При этом резонно считается, что стоимость объекта очень сильно зависит от локации, и эта локация — примерно квартал. Ну и что вы думаете — сидят люди, и руками эти самые кварталы рисуют, потому что нигде их нет, ни на одной карте. Теоретически можно попытаться построить их по улицам, например, но эта задачка как минимум уровня построения маршрутов, и я пока слабо себе представляю, как ее решать.

                                А наша задачка — она не то что сложнее, она просто другая. Прямое геокодирование, если его делать по уму, включает приличный объем NLP. И это совсем не геометрия.

                  Only users with full accounts can post comments. Log in, please.