Привет Хабр. История начинается с того что мы решили сделать гео сервис с возможностью размещения меток на карте самими пользователями.
И когда решили залить в базу 1 миллион маркеров то поняли, что даже если запрашивать маркеры только в определенном радиусе то все работает очень медленно и кластеризация на клиенте тоже не вариант :)

А где-то под этим лесом находится манхетен




Кластеризация это очень интересная задача, есть клиентская, есть серверная и все крупные картографические сервисы ее поддерживают. Известно множество видов кластеризации, но я в своих поисках выбрал Grid-кластеризацию. Про нее расказывали сотрудники Яндекса в нескольких своих докладах на конференциях. Она используется в яндекс картах для кластеризации маркеров.
Yandex кластеризация


Немного теории


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

Термины:

Тайл(tile, плитка) — то по сути просто разделение карты на квадраты для определенного zoom-а. На первом уровне у нас карта будет состоять из 1 квадрата(тайла) на 0 уровне, на 1 уровне из 4-х квадратов(тайлов), на втором уровне на 16 на третьем на 64 и т.д. По сути это 4 в степени уровня.


Quadkey — это уникальный идентификатор тайла на любом уровне. На первом уровне нумеруем тайлы 0,1,2,3 справа налево, сверху вниз. Затем на следующем уровне для каждого тайла из верхнего уровня делаем такую же нумера но уже добавляя quadkey родителя т.е. 00,01,02,03,10,11,12,13 и т.д. Более понятно это становится по рисунку:

Посути мы получили число в 4-чной системе исчисления. Теперь у нас каждый тайл имеет свой уникальный идентификатор на любом уровне zoom-а. Из zoom-а, широты и долготы можно получить уже сам quadkey тайла т.е. любую точку можно преобразовать в quadkey соответсвующего уровня zoom-а. А вот так будет выглядеть quadkey для конкретных координат:
longitude -79.377807617187500
latitude 43.653785705566406
quadkey 13940830302567 (base4: 03022313122033033011213)

Описание этого подхода на англиской языке от Microsoft:
Bing Maps Tile System
В этой статье реализации всех этих преобразований написана на C#, я побыстрому переписал ее на PHP:
TileSystem.php
class TileSystem
{
  const EarthRadius = 6378137;
  const MinLatitude = -90;
  const MaxLatitude = 90;
  const MinLongitude = -180;
  const MaxLongitude = 180;
  const MaxZoom = 23;

  private static function clip( $n, $minValue, $maxValue)
  {
    return min( max( $n, $minValue), $maxValue);
  }

  public static function mapSize( $levelOfDetail)
  {
    return 256 << $levelOfDetail;
  }

  public static function GroundResolution( $latitude, $levelOfDetail)
  {
    $latitude = self::clip( $latitude, self::MinLatitude, self::MaxLatitude);
    return cos( $latitude * pi() / 180) * 2 * pi() * self::EarthRadius / self::mapSize( $levelOfDetail);
  }

  public static function mapScale( $latitude, $levelOfDetail, $screenDpi)
  {
    return self::groundResolution( $latitude, $levelOfDetail) * $screenDpi / 0.0254;
  }

  public static function latToPixelY( $latitude, $levelOfDetail)
  {
    $latitude = self::clip( $latitude, self::MinLatitude, self::MaxLatitude);
    $sinLatitude = sin( $latitude * pi() / 180);
    $y = 0.5 - log((1 + $sinLatitude) / (1 - $sinLatitude)) / (4 * pi());
    $mapSize = self::mapSize( $levelOfDetail);
    $pixelY = (int) self::clip( $y * $mapSize + 0.5, 0, $mapSize - 1);
    return $pixelY;
  }

  public static function longToPixelX( $longitude, $levelOfDetail)
  {
    $longitude = self::clip( $longitude, self::MinLongitude, self::MaxLongitude);
    $x = ($longitude + 180) / 360;
    $mapSize = self::mapSize( $levelOfDetail);
    $pixelX = (int) self::clip( $x * $mapSize + 0.5, 0, $mapSize - 1);
    return $pixelX;
  }

  public static function pixelXToLong( $pixelX, $levelOfDetail)
  {
    $mapSize = self::mapSize( $levelOfDetail);
    $x = ( self::clip( $pixelX, 0, $mapSize - 1) / $mapSize) - 0.5;
    $longitude = 360 * $x;
    return $longitude;
  }

  public static function pixelYToLat( $pixelY, $levelOfDetail)
  {
    $mapSize = self::mapSize( $levelOfDetail);
    $y = 0.5 - (self::clip( $pixelY, 0, $mapSize - 1) / $mapSize);
    $latitude = 90 - 360 * atan( exp(-$y * 2 * pi())) / pi();
    return $latitude;
  }

  public static function pixelYToTileY( $pixelY)
  {
    return $pixelY / 256;
  }

  public static function pixelXToTileX( $pixelX)
  {
    return $pixelX / 256;
  }

  public static function tileXYToPixelXY( $tileX, $tileY, &$pixelX, &$pixelY)
  {
    $pixelX = $tileX * 256;
    $pixelY = $tileY * 256;
  }

  public static function tileXYToQuadKey( $tileX, $tileY, $levelOfDetail)
  {
    if ($levelOfDetail == 0)
      return "0";

    $quadKey = "";
    for ($i = $levelOfDetail; $i > 0; $i--) {
      $digit = 0;
      $mask = 1 << ($i - 1);
      if (($tileX & $mask) != 0) {
        $digit++;
      }
      if (($tileY & $mask) != 0) {
        $digit++;
        $digit++;
      }
      $quadKey .= $digit;
    }
    return $quadKey;
  }

  public static function latLongToQuadKeyDec( $lat, $long, $zoom)
  {
    return self::quadKey4ToDec( self::latLongToQuadKey( $lat, $long, $zoom));
  }

  public static function latLongToQuadKey( $lat, $long, $zoom)
  {
    $pixelX = \foci\utils\TileSystem::longToPixelX( $long, $zoom);
    $pixelY = \foci\utils\TileSystem::latToPixelY( $lat, $zoom);
    $tileX = \foci\utils\TileSystem::pixelXToTileX( $pixelX);
    $tileY = \foci\utils\TileSystem::pixelYToTileY( $pixelY);
    return \foci\utils\TileSystem::tileXYToQuadKey( $tileX, $tileY, $zoom);
  }

  public static function pixelXYToQuadKey( $pixelX, $pixelY, $zoom)
  {
    $tileX = \foci\utils\TileSystem::pixelXToTileX( $pixelX);
    $tileY = \foci\utils\TileSystem::pixelYToTileY( $pixelY);
    return \foci\utils\TileSystem::tileXYToQuadKey( $tileX, $tileY, $zoom);
  }

  public static function pixelXYToQuadKeyDec( $pixelX, $pixelY, $zoom)
  {
    return base_convert( self::pixelXYToQuadKey( $pixelX, $pixelY, $zoom), 4, 10);
  }

  public static function quadKeyToTileXY( $quadKey, &$tileX, &$tileY, &$levelOfDetail)
  {
    $tileX = $tileY = 0;
    $levelOfDetail = strlen( $quadKey);
    for( $i = $levelOfDetail; $i > 0; $i--)
    {
      $mask = 1 << ($i - 1);
      switch( $quadKey[$levelOfDetail - $i])
      {
        case '0':
          break;

        case '1':
          $tileX |= $mask;
          break;

        case '2':
          $tileY |= $mask;
          break;

        case '3':
          $tileX |= $mask;
          $tileY |= $mask;
          break;

        default:;
      }
    }
  }

  public static function tileXYToQuadKeyDec( $tileX, $tileY, $levelOfDetail)
  {
    return base_convert( self::tileXYToQuadKey( $tileX, $tileY, $levelOfDetail), 4, 10);
  }

  public static function quadKey4ToDec( $quadkey)
  {
    return base_convert( $quadkey, 4, 10);
  }
}


Переходим к практике


Область видимости


Тут надо учесть важный момент, не нужно запрашивать все кластеры на карте на конкретном zoom, а нужно получать только те, которые входят в область видимости клиента(в нашем случае это android и web).

Варианты получения кластеров в области видимости:
  • Умный клиент. Клиент сам определяет какой zoom ему нужен и какие кластеры войдут в область его видимости и запрашивает их все сразу например при инициализации view либо запрашиват только те которые вошли в область видимости при перемещении по карте
  • Ленивый клиент. Клиент просто отсылает свою область видимости на сервер и отображет то что будет прислано сервером.

Мы выбрали Ленивый клиент. Почему?
  • Логика клиента очень проста. Посути клиент работает как view, что получил то и отображает.
  • Это простота удобна если как у нас есть несколько клиентов на разных платформах. Не надо дублировать хитрую логику определения фокусов в области видимости. Вся логика сосредоточина а одном месте — на сервере.
  • Как показала практика требуется малое время на реализацию.

Какова будет логика сервера?
  • На сервер приходят координаты углов прямоугольной области видимости( например верхний левый угол, нижний правый угол).
  • Сервер зная размеры этой прямоугольной области определяет тайл какого размера в нее входит полностью(и по ширине и по высоте)
  • Дальше сервер зная размер тайла может определить zoom, поэтому его передавать вместе с координатами нет необходимости.
  • Затем надо учесть что этот zoom надо увеличить, для более удобного отображения.
  • И последним сервер опеределит quadkey-и попадающие в область видимости, сделает выборку кластеров и отдаст клиенту.

Хранение


Для хранения quadkey-ев нам надо выбрать уровень zoom, возьмем 23, его будет вполне достаточно. БД будем использовать MySQL. Таблица хранения маркеров будет иметь следующую структуру:
CREATE TABLE `marker` (
  `id` int(11) NOT NULL AUTO_INCREMENT,
  `latitude` FLOAT(19,15) NOT NULL,
  `longitude` FLOAT(19,15) NOT NULL,
  `quadkey` bigint(20) NOT NULL DEFAULT 0,
  PRIMARY KEY (`id`),
  INDEX `quadkey ` (`quadkey `)
) ENGINE=InnoDB DEFAULT CHARSET=utf8;

Индекс по quadkey позволяет нам легко делать выборки количества маркеров в тайле(кластере).

Для хранения quadkey в MySQL можно использовать bigint(20), просто переводим из 4-чной системы в десятичную и в БД и приложении работаем с 10-чным quadkey.

Код получения кластеров


Код получения кластеров будет выглядеть примерно так:
  public function getClusters( $quad_key, $zoom)
  {
    $q1 = $quad_key . str_repeat("0", TileSystem::MaxZoom - $zoom);
    $q2 = $quad_key . str_repeat("3", TileSystem::MaxZoom - $zoom);
    $mask = $q1;
    $mask_plus = 1;
    $mask = $quad_key . str_repeat("3", $mask_plus);
    $mask = $mask . str_repeat("0", TileSystem::MaxZoom - $zoom - $mask_plus);

    $q1 = TileSystem::quadKey4ToDec( $q1);
    $q2 = TileSystem::quadKey4ToDec( $q2);
    $mask = TileSystem::quadKey4ToDec( $mask);

     return $this-> findBetweenQuadKeys( $q1, $q2, $mask);
  }

 public function findBetweenQuadKeys( $quad_key_l, $quad_key_r, $mask)
  {
    $entities = [];

    $sql = "SELECT `quadkey`, MIN(id) AS id, COUNT(id) AS count, AVG(longitude) AS longitude, AVG(latitude) AS latitude ";
    $sql .= "FROM `marker ` ";
    $sql .= "WHERE `quadkey` BETWEEN $quad_key_l AND $quad_key_r ";
    $sql .= "AND `is_deleted`= '0' ";
    $sql .= "GROUP BY (`quadkey` & $mask) ";

    if( $stmt =  $this->_connection->query( $sql))
    {
      foreach( $stmt as $row)
        $entities[] = $this->_createEntityFromRow( $row);
    }

    return $entities;
  }

Переменная $mask позволяет нам группировать quadkey-и c определенным zoom. Переменные $quad_key_l и $quad_key_r позволяют учитывать quadkey-и всех маркеров для данного родителя т.е если $quad_key_l=03022313122033033000000 и $quad_key_r=03022313122033033033333 то мы найдем в базе все маркеры для 030223131220330330.
Так будет выглядеть результирующий sql запрос:
SELECT `quadkey`, MIN(id) AS id, COUNT(id) AS count, AVG(longitude) AS longitude, AVG(latitude) AS latitude
FROM marker 
WHERE quadkey between 15499682906112 and 15499683168255 
GROUP BY (quadkey & 15499683151872);

-- более понятная 4-чня форма
-- WHERE quadkey between "03201203031012000000000" and "03201203031012333333333"
-- GROUP BY (quadkey & "03201203031012330000000")


Что получилось


К сожалению скриншотов с кластерами на базе в 1 миллион маркеров не сохранилось(хотя это было нечто), но есть текущая небольшая база, и вот как это выглядит:
web

Android


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

Результаты


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

Источники


www.mapbox.com/guides/how-web-maps-work
msdn.microsoft.com/en-us/library/bb259689.aspx

Комментарии (29)


  1. kashey
    14.10.2015 08:57
    +3

    1. Конкретно этот quadkey — Morton code.
    2. 16ый уровень разбиения, который можно сохранить в 32 битный uint описывает квадрат со стороной примерно 300 метров. Для предварительной кластеризации этого достаточно.
    3. Ленивый клиент — это не самое умное решение, но за перевод в пиксельные координаты отдельное спасибо.
    4. Маску лучше представлять как маску. Т.е. 03201203031012000000000/12.

    В общем случае эта задача называется «понижением размерности» 2Д данных через использование какой либо spatial filling curve, которая сводит задачу к 1Д варианту. Люблю эту наркоманию.


    1. kashey
      14.10.2015 09:04
      +3

      В статье W for Wikipedia, в разделе Z for Z-code есть чуть более техническое описание работы Z кривых, плюс ссылки на sql+php и mongo+js решения для серверной кластеризации.
      Вообще об этой технологии я пытаюсь рассказывать уже лет 8, жаль что до нужных людей мои презентации не доходят :(


      1. guyfawkes
        14.10.2015 09:46

        А где можно увидеть Ваши презентации?


        1. kashey
          14.10.2015 10:12

          Последняя — events.yandex.ru/lib/talks/2354 — 2014 год, первая нормальная была 2012.codefest.ru/lecture/77 — 2012 год, даже примерчик на гитхаб выложил — github.com/theKashey/quadLoader.
          В принципе «читабельно» в первый раз я описал Z и кластеризацию в статье про gdeetotdom(2010), которая сейчас уехала на geektimes. С тех пор в каждой второй статье на хабре я это дело упоминаю, так и или иначе — уж больно широкие применения у этих чтук.
          Это не просто «число» — одновременно это nested set, B-tree индекс и материализация Linear quad-tree. И это самая простая кривая — их на любой взыскательный вкус напридумывали.


          1. xanm
            14.10.2015 11:31

            Да действительно, ваша презентация была как раз одним из источников. Еще на HighLoad 2014 девушка из яндекс карт рассказывала по этой теме.


    1. xanm
      14.10.2015 09:51

      Отлично, мого полезной информации. Спасибо. Когда искал что-то подобное найти было сложновато и наткнулся только на обрывки информации.
      А на счет ленивого клиента — это был выбран просто самый быстрый способ реализации. Естественно в дальнейшем от него надо уходить.


      1. kashey
        14.10.2015 10:15

        На самом деле, лично использую quad коды в основном для выборки обьектов, а не для кластеризации. Это так — приятный бонус, да и не все кластеризуется. Да и на клиенте можно поработать немного :)
        В частности Z коды используются не только для загрузки данных из базы — именно они (лично у меня) используются для передачи команды на сервер.
        В первый раз такой подход был замечен на wikimapia, я его активно перенял. С тех пор особо популярен он не стал — если и грузят по тайловой сетке, то в обычных x,y,z координатах.


        1. xanm
          14.10.2015 10:24

          Нам удалось еще использовать их для посика кластеров с определенными тегами маркеров.


    1. Fesor
      15.10.2015 01:32

      ммм а чем квадкоды отличаются от геохэшей?


  1. ESQUELETO
    14.10.2015 10:56

    А почему не используют K-Means кластеризацию? Или это дорого на сервере?

    Это ведь суперлениво для клиента — посылаешь на сервер область видимости и максимальное желаемое количество результирующих кластеров и получаешь ответ.


    1. xanm
      14.10.2015 11:34

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


    1. kashey
      14.10.2015 11:52

      K-means достаточно долгий, в CPU, вариант, но для многих бэкендеров нет ничего проще, чем нагрузить кластер спарков его расчетом…
      Но применимость k-means на картах просто низка. У него есть проблема — он не «стабильный».
      Например вы смотрите на Москву. Запрашиваете данные только для Москвы. Но… вы должны произвести кластеризацию всех обьектов по всему миру на 10ом зуме. Чтобы при сдвиге у вас не изменился набор исходных данных и все «взрывообразно» не перекластеризовалось.
      Потом приходит необходимость делать иеархический k-means, и ты понимаешь что это того не стоит.


      1. Fesor
        15.10.2015 01:27

        Потом приходит необходимость делать иеархический k-means, и ты понимаешь что это того не стоит.

        сделали, жалеем)


  1. Gard
    14.10.2015 11:04

    Не понятно какая задача решается:
    если нужно правильно хранить, то лучше взять хранилище с поддержкой геометрических индексов. Они есть уже и в SQL и в noSQL решениях. C MySQL не работал, поэтому не могу сказать, что есть в нём, но в среде открытых есть PostGIS.
    если нужна именно кластеризация, то грид быстрый, и самый не красивый из всех алгоритмов. Когда данных много, пины кластеров образуют практически регулярную сетку. Это хорошо заметно на скрине Яндекса в статье. Раз кластеры считаются на сервере, то можно выбрать более ресурсоёмкий алгоритм с красивым результатом.
    Если посмотреть на скрины того «что получилось»: алгоритм выбора репрезентативной точки выбран не удачно. Например, маркер который, видимо, должен принадлежать Кубе просто висит над морем.


    1. xanm
      14.10.2015 11:28

      Задача решалась реализовать выбранный способ кластеризации.
      Да, в других БД есть поддержка хранения гео данных, но вот на счет кластеризации я не уверен. И здесь просто вырезки из большого проекта где используется MySQL.
      Естественно что алгоритм не идеален, понятно что мы находим средную координату в квадрате и она может посать и на вершину горы и в другю страну и в озеро и в море. Все эти недостатки можно устранить в будущем.


      1. kashey
        14.10.2015 11:54

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


        1. xanm
          14.10.2015 12:01

          Да, это очень простой вариант, так мы всегда попадем в реальную точку.


  1. smartnlg
    14.10.2015 11:21

    Если не секрет… Где вы использовали данную технологию? Какие задачи решали?


    1. xanm
      14.10.2015 11:38

      Использовали в одном гео-сервисе, название сказать не могу(посчитают пиаром).
      Задача была размещать на карте большое количество маркеров( миллионы), и отображать это пользователям без проблем с производительностью.


      1. lifestar
        14.10.2015 11:42

        Не посчитают, не боись.
        Так где же?


        1. xanm
          14.10.2015 11:57
          +1

          Demo версия карты fociapp.com/map. Там же есть ссылка на приложение. А вообще планируетс статья по этому проекту в разделе «Я пиаруюсь».


          1. kashey
            14.10.2015 12:33

            Немного хромает реализация — тут и метки «мигают» при сдвиге карты и кластеризация немного подтупливает.
            Я бы с сервера отдавал бы больше точек — еще бы на один уровень маску группировки сдвинул бы и «докластеризовывал» бы на клиенте.
            Ну и лимиты при выборке из базы увеличил бы (или совсем убрал) — количество данных регулируется кластеризацией.
            Сейчас реально бывают «дырки» — в европе данные есть, в Австралии — нет. Убираем европу, и внезапно точки появляются. От этого еще хорошо разбиение запроса на тайлы лечит.


            1. xanm
              14.10.2015 12:44

              Согласен, уже обдумывал эти правки, но пока немного не хватает времени. Сама клиентская часть еще достаточно срая.


  1. siberiano
    14.10.2015 12:33

    Очень хочу увидеть сравнение этого приложения с Postgresql/PostGIS и идекса типа R-tree. Что быстрее?

    Просто в MySQL вся поддержка геоданных сделана по остаточному принципу и to be implemented. В PostGIS она полноценная и работает очень хорошо.


    1. xanm
      14.10.2015 12:47

      Сделать сравнение было бы доастаочно интересно, попробую найти на это время :)
      В MySQL рассматривал Spatial Data, но там насколько помню они еще не готовы к полноценному использованию.


      1. siberiano
        14.10.2015 14:35

        Вот я в 2009 их пробовал, тоже не были готовы. :)


    1. Fesor
      15.10.2015 01:02

      Из опробованных решений еще стоит упоминуть ElasticSearch и плагин geocluster-facet. Ну и эластика умеет агрегацию по сетке делать из коробки (использую геохэши).


  1. duke_pro
    15.10.2015 16:43

    Кстати на практике оказалось что использование 50% медианы куда лучше выглядит чем среднее арифметическое.


    1. xanm
      15.10.2015 17:57

      Возможно, надо будет попробовать.