Точное определение региона по GPS координатам

При разработке одного приложения возникла проблема разграничения доступа для регионов.

Встала проблема определения принадлежности объекта к какому-либо региону России по его GPS координатам

Первое, что мы начали использовать — это API Google,
после того как прописали алиасы к возвращаемым строкам и оплаты доступа (чтобы убрать лимит на запросы) все заработало.
И все было нормально пока гугл не сменил выдачу, например было раньше: Moskovskaya oblast', стало Moscow oblast'
Тут то и было решено не надеяться на гугл, а определять регион своими силами.




Первым, что пришло на ум было: распарсить возвращаемые данные Яндекса (в яндекс-картах версии 1.x был модуль для работы с регионами)
Так и поступили, получилось 2 таблицы.
В одной хранились названия регионов, во второй координаты многоугольника, описывающие регион. (83 — регионы и 8500 — точки )

Был найден алгоритм, для проверки лежит ли точка в многоугольнике. Код был переписан с Си на PHP, разбирать как он работает нет смысла, он просто работает.

$sx,$sy — координаты точки, которую проверяем
$coords — координаты точек выпуклого многоугольника (отсортированы в порядке обхода):
$coords = array(
array('x'=> 66.6634, 'y' => '66.4433'),
и тп.
)
$x, $y — названия ключей в массиве

public static function into_poly($sx, $sy, &$coords, $x='x', $y='y')
{
    Profiler::start('Detect collision', __FUNCTION__);
    $pj=0;
    $pk=0;
    $wrkx=0;
    $yu = 0;
    $yl = 0;
    $n = count($coords);
    for ($pj=0; $pj<$n; $pj++)
    {
        $yu = $coords[$pj][$y]>$coords[($pj+1)%$n][$y]?$coords[$pj][$y]:$coords[($pj+1)%$n][$y];
        $yl = $coords[$pj][$y]<$coords[($pj+1)%$n][$y]?$coords[$pj][$y]:$coords[($pj+1)%$n][$y];
        if ($coords[($pj+1)%$n][$y] - $coords[$pj][$y])
            $wrkx = $coords[$pj][$x] + ($coords[($pj+1)%$n][$x] - $coords[$pj][$x])*($sy - $coords[$pj][$y])/($coords[($pj+1)%$n][$y] - $coords[$pj][$y]);
        else
            $wrkx = $coords[$pj][$x];
        if ($yu >= $sy)
            if ($yl < $sy)
            {
                if ($sx > $wrkx)
                    $pk++;
                if (abs($sx - $wrkx) < 0.00001) return 1;
            }
        if ((abs($sy - $yl) < 0.00001) && (abs($yu - $yl) < 0.00001) && (abs(abs($wrkx - $coords[$pj][$x]) + abs($wrkx - $coords[($pj+1)%$n][$x]) - abs($coords[$pj][$x] - $coords[($pj+1)%$n][$x])) < 0.0001))
            return 1;
    }
    if ($pk%2)
        return 1;
    else
        return 0;
}

Пример вызова:
$coords = array(
array('lng'=> 66.6634, 'lat' => '66.4433'),
array('lng'=> 66.6534, 'lat' => '66.4433'),
array('lng'=> 66.6434, 'lat' => '66.4433'),
array('lng'=> 66.6334, 'lat' => '66.4433'),
);
$in = into_poly(66.4455, 66.2255, &$coords, $x='lng', $y='lat');


Вернет true, если точка лежит внутри многоугольника.

Теперь просто перебором областей для точки мы получаем ответ.

Всё бы было хорошо, но 8500 точек границ регионов для России — это очень мало. Разброс получится +- 50 км

Чтобы получить точность больше – нам нужно больше точек границ регионов.
Поискав на просторах был найден ресурс, который нам свободно отдает эти точки: gis-lab.info/qa/rusbounds-rosreestr.html (спасибо энтузиастам за вашу работу)
Я почему-то выбрал формат KML (наверное потому, знал чем его открывать — это формат, который поддерживает Google Earth)
Распарсил данные и теперь таблицы получились пожирнее. Таблица с координатами точек границ регионов стала 620 000 точек (34 мб)

Теперь если проверять по старому алгоритму принадлежность точки к многоугольнику региона — на стареньком Core 2 Duo — около 70 секунд.
Не вариант. Алгоритм нужно оптимизировать:

Например, Москва:
Имеем 2,064 точек многоугольника (на самом деле несколько многоугольников):


Но сейчас нам достаточно определить всего-лишь возможность, что точка будет в этой области. Поэтому мы воспользуемся алгоритмом обхода точек и построения выпуклого многоугольника.
algolist.manual.ru

вот его реализация на PHP:

public static function sort_points($mass, $x = 'x', $y = 'y'){
    $current=0;
    $next=0;
    $p1 = $mass[0];
    $mass2 = array();
    //определяем первую точку
    for ($i=1; $i<count($mass); $i++){
        if ( ($p1[$y] > $mass[$i][$y]) || ($p1[$y] == $mass[$i][$y] && $p1[$x] < $mass[$i][$x]) ) {
            $p1 = $mass[$i];
            $current = $i;
        }
    }
    $n = count($mass);

    $p0 = $p1;
    $p0[$x]--;

    $mass2[] = $mass[$current];
    $first = $current;

    //начинаем перебор всех элементов

    do{
        $cmax_not_set=1;
        for ( $i=0; $i<$n; $i++ ) {
            $key = $i;
            //продолжаем если такой элемент уже выбран
            if ( $mass[$current][$x] == $mass[$key][$x] && $mass[$current][$y] == $mass[$key][$y] ) continue;
            //берем 1ю точку
            if ($cmax_not_set) {
                $next = $key;
                $cmax_not_set=0;
                continue;
            }

            $v1_x = $mass[$key][$x] - $mass[$current][$x];
            $v1_y = $mass[$key][$y] - $mass[$current][$y];

            $v2_x = $mass[$next][$x] - $mass[$current][$x];
            $v2_y = $mass[$next][$y] - $mass[$current][$y];

            //магия
            $c = $v1_x * $v2_y - $v1_y * $v2_x;

            if ( $c > 0 ) $next = $key;
           // if ( ($c==0) && ( self::_dist($mass[$current], $mass[$key], $x, $y) > self::_dist($mass[$current], $mass[$next], $x, $y) ) ) $next = $key;

        }
        //записываем найденный элемент в массив
        $mass2[] = $mass[$next];
        $current = $next;
    }while($first != $next);
    return $mass2;
}

На выходе мы получим массив точек выпуклого многоугольника.
Для Москвы; 19 точек:


На рисунке первым слоем отмечены все 2500 точек Москвы + слой, который мы получили после расчета выпуклого многоугольника.
Сгенерируем для каждого региона такую область.

Теперь доработаем наш алгоритм и будем сначала проверять возможность, что точка может лежать в области – таких областей может быть несколько.
Теперь после прогона всех областей, мы прогоним области, в которых возможно, что лежит точка, алгоритмом описанным выше.
Итого: 0.177сек на слабеньком Core 2 Duo

На самом деле есть еще одна сложность – это анклавы (Москва лежит внутри Московской области, как и Питер, было решено использовать приоритеты, сначала проверяется лежит ли точка в Москве, а потом уже в Московской области)