Accurate determination of the region GPS coordinates

If you are developing one application there is a problem of access for the regions.

/ > The problem of determining the membership of an object to any region of Russia by its GPS coordinates

The first thing we started using is Google APIs
after prescribed aliases to the return lines and payment (to remove the limit on the queries) it worked.
And everything was fine until Google changed the results, such as it was before: Moskovskaya oblast', became Moscow oblast'
Then it was decided not to rely on Google, and to define the region on their own.




The first thing that came to mind was to parse the returned data of Yandex (Yandex-maps version 1.x had a module to work with the regions)
Did so, turned out 2 of the table.
In one were kept the names of the regions, the second coordinates of the polygon describing the region. (83 — regions and the 8,500 — point )

Has been found the algorithm to check whether a point lies in a polygon. The code was rewritten from C to PHP, to disassemble how it works makes no sense, it just works.

$sx,$sy — coordinates of a point that you check the
$coords — coordinates of points of a convex polygon (sorted in order):
$coords = array(
array('x'=> 66.6634, 'y' => '66.4433'),
etc.
)
$x, $y — names of the keys in the array
the

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($version$, 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;
}

Example call:
$coords = array(
array('lang'=> 66.6634, 'lat' => '66.4433'),
array('lang'=> 66.6534, 'lat' => '66.4433'),
array('lang'=> 66.6434, 'lat' => '66.4433'),
array('lang'=> 66.6334, 'lat' => '66.4433'),
);
$in = into_poly(66.4455, 66.2255, &$coords, $x='lng', $y='lat');


Returns true if the point lies inside the polygon.

Now just by iterating over the fields to the point we get the answer.

Everything would be fine, but 8500 points of boundaries of regions in Russia is very small. Scatter get +- 50 km

To gain more precision we need more points of boundaries of regions.
Looking at the open spaces were found a resource that freely gives us these points: gis-lab.info/qa/rusbounds-rosreestr.html (thanks to the enthusiasts for your work)
For some reason I chose KML (probably because knew how to open it — it is a format that supports Google Earth)
Was responsil data and now the table has turned meatier. Table with coordinates of points of boundaries of regions was 620 000 points (34 MB)

Now if you check the old algorithm belonging of a point to the polygon of the region — on an old Core 2 Duo is about 70 seconds.
Not an option. The algorithm should be optimized:

For example, Moscow:
Have 2,064 points of the polygon (in fact, some polygons):


But now it is enough to determine only the possibility that the point will be in this area. Therefore, we use the algorithm of bypass of points and build a convex polygon.
algolist.manual.ru

here implementation on PHP:
the

public static function sort_points($mass, $x = 'x', $y = 'y'){
$current=0;
$next=0;
$p1 = $mass[0];
$mass2 = array();
//define the first point
for ($i=1; $i<count($mass); $i++){
if ( ($p1[$y] > $mass[$i][$y]) || ($p1[$y] == $mass[$i][$y] && $p1[$x] < $mass[$i][$x]) ) {

$current = $i;
}
}
$n = count($mass);

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

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

//start a loop through all the items

do{
$cmax_not_set=1;
for ( $i=0; $i<$n; $i++ ) {
$key = $i;
//continue if element already selected
if ( $mass[$current][$x] == $mass[$key][$x] && $mass[$current][$y] == $mass[$key][$y] ) continue;
//get 1st point
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];

//magic
$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;

}
//writable element found in the array
$mass2[] = $mass[$next];
$current = $next;
}while($first != $next);
return $mass2;
}

At the output we get an array of points of a convex polygon.
For Moscow; 19 points:


In the figure the first layer marked all 2500 points of Moscow + layer, which we received after the calculation of a convex polygon.
Generate for each region a region.

Now we will modify our algorithm will first test the possibility that the point may lie in the region – such regions may be several.
Now after running all of the areas we drive region in which that point is located, the algorithm described above.
Total: 0.177 sec on weak Core 2 Duo

Actually there is another difficulty – the enclaves (Moscow inside Moscow region, as Peter, it was decided to use the priorities, the system first checks whether the point lies in Moscow, and then Moscow region)

Article based on information from habrahabr.ru