Алгоритм поиска ближайших объектов в БД

Теперь обо всем немного поподробнее.

1. Получить и обработать параметры из URL'а.

В моем случае все параметры передавались в строке запроса, поэтому доступ к ним я получал через специальный "суперглобальный" массив $_GET.

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

Реализовать данную задачу можно по-разному, лично я написал класс Params, который помогал мне в этом:

class Params {
    // Параметры по умолчанию и он же список разрешенных параметров
    private $paramsDefault = array(
        'point' => null,
        'distance' => 2,
        'format' => 'ymapsml',
        'jsoncallback' => null
    );

    // Списки доступных значений
    private $enabledParams = array(
        'point' => 'validPoint',
        'distance' =>  'vaildDistance',
        'format' => array('ymapsml', 'json'),
    );
   
    // Параметры
    public $params = null;
   
    function __construct ($params = array()) {
        $this->params = array_merge($this->paramsDefault, $params);         
        $this->params = $this->filter($this->params);
        //$this->printParams();
    }
   
    function __destruct () {}
   
    // Метод оставляет только те параметры, которые есть в списке $paramsDefault и валидны (см. метод isValid)
    private function filter ($params) {}

    // Проверяет валидность значения параметра
    private function isValid ($key, $value) {}
}

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

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

2. По центру и заданной окружности рисуем ограничивающий многоугольник.

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

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

Форму Земли я принял как шар. Если взять такую апроксимацию, то задачи геодезии решаются гораздо проще.

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

Все задачи решает один класс GeodeticTask:

// Класс для решения задач геодезии
class GeodeticTask {

    // Решает прямую задачу геодезии
    function solveDirect ($point, $distance, $angle) {
        $arrPoint = explode(',', $point);
        $angle = self::getAngle($angle);

        // Определяем приращения
        $dx = $distance / (111.2 * 0.8 * cos($arrPoint[1])) * cos($angle);
        $dy = $distance / 111.2 * sin($angle);
       
        return implode(' ', array(
            $arrPoint[0] + $dx,
            $arrPoint[1] + $dy
        ));
    }

    // Из строки получаем угол (временная функция)
    // Нет ничего постоянного, чем временное
    private function getAngle ($str) {
        $arrAngle = explode(',', $str);
        return self::degToRad($arrAngle[0] + $arrAngle[1] / 60 + $arrAngle[2] / 3600);
    }
}

 С помощью метода solveDirect() класса GeodeticTask легко найти координаты ограничивающего многоугольника.

class Polygon {
    private $center; // Центр многоугольника
    private $radius; // Радиус описанной окружности
   
    // Массив координат
    private $points = array();
   
    // Строка с координатами
    public $result = '';

    // На вход передаются параметры урла
    function __construct ($params) {
        if ($point = $params->get('point')) {
            $this->center = $point;
            $this->radius = $params->get('distance');

            $this->create();
        }
    }

    function create () {
        for ($i = 0; $i < 360; $i += 45) {
            array_push($this->points, GeodeticTask::solveDirect($this->center, $this->radius, $i.',0,0.00'));
        }
        // Замыкаем многоугольник
        array_push($this->points, GeodeticTask::solveDirect($this->center, $this->radius, '0,0,0.00'));
       
        $this->result = implode(',', $this->points);
    }
}

Полученные координаты из многоугольника я сохраняю в строку. Это нужно для формирования запроса к БД.

Для угла и хранения точек лучше создать свои структуры данных, а не хранить строками. Я делал временное решение, которое так и осталось временным :)

3. Из БД выбираем все точки, который попали в ограничивающий многоугольник.

Для хранения координат точечных объектов (как, например, станций метро) используется стандартный тип данных POINT (spatial index) в MySQL. Соответственно данные хранятся в виде R-tree.

Можно было хранить отдельно координаты как два разных поля. Такой метод подходит, если вы всегда выбираете данные из прямоугольной области. Если есть желание делать выборку из более сложных областей (как в моем случае), то лучше все же использовать spatial index. Возможно, что у вас есть более красивое решение? С удовольствием бы послушал.

В итоге результат на выборку данных выглядит примерно следующим образом:

SELECT s.title, asText(s.coords) as coords, s.address, l.color
FROM `stations` as s, `lines` as l
WHERE MBRContains(
    GeomFromText('
        Polygon((
            ".$polygon->result."
        ))'
    ) , coords) AND l.id = s.line_id

Переменная $polygon - это объект класса Polygon, который был описан выше. Результат этого запроса поместим в переменную $stations, с которой и будем далее работать.

4. Находим расстояния от центра многоугольника до всех найденных объектов.

Для того, чтобы реализовать эту задачу - необходимо будет решить еще одну задачку геодезии - обратную.

Обратная геодезическая задача (ОГЗ) заключается в определении по геодезическим координатам двух точек на земном эллипсоиде длины и дирекционного угла направления между этими точками.

В класс GeodeticTask добавляем еще один метод solveInverse(), который будет решать обратную задачу геодезии. Алгоритм я взял с сайта gis-lab.info. Вот так он примерно будет выглядеть:

function solveInverse ($point1, $point2) {
    $lat1 = self::degToRad($point1[1]);
    $lat2 = self::degToRad($point2[1]);
    $lng1 = -1 * self::degToRad($point1[0]);
    $lng2 = -1 * self::degToRad($point2[0]);

    $dlon_E = mod($lng2 - $lng1, 2 * pi()) ;
    $dlon_W = mod($lng1 - $lng2, 2 * pi());

    $dphi = log(tan($lat2 / 2 + pi() / 4) / tan($lat1 / 2 + pi() / 4));
    $q = (abs($lat2 - $lat1) > 0.000000001) ? ($lat2 - $lat1) / $dphi : cos($lat1);

    $distance = 0;
   
    if ($dlon_E < $dlon_W) {
        $distance = 6372795 * pow((pow($q, 2) * pow($dlon_E, 2) + pow(($lat2 -  $lat1), 2)), 0.5);
    } else {
        $distance = 6372795 * pow(pow($q, 2) * pow($dlon_W, 2) + pow($lat2 - $lat1, 2), 0.5);
    }
    return $distance;
}

private function degToRad ($deg) {
    return (pi() * $deg) / 180;
}

private function radToDeg ($rad) {
    return ($rad * 180) / pi();
}

 Также понадобится функция самописная функция mod (стандартная не подходит):

function mod ($delimoe, $delitel) {
    return $delimoe - ($delitel * floor($delimoe / $delitel));
}

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

5. Выбираем ближайшие n.

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

6. Создаем XML-документ.

7. Транслируем его с помощью XSLT в любой нужный формат.