Raycasting algorithm with GPS coordinates

Question

I'm making a small application with Google maps that determines if an entered address is part of a predefined service region.

The user enters an address, and a PHP script fetches the lat/long from the Geocoding API and applies raycasting with a bunch of coordinates that make up the vertices of the region (taken from a KML file Maps generated).

The problem is this: it works most of the time, but some addresses outside the service region incorrectly report as eligible, while some others inside the region are ineligible. At first I thought this was a precision issue with Google maps, but the coordinates generated from addresses in the Geocoding service are spot-on accurate. It probably has something to do with the formula.

Here it is (it's based off of code I found elsewhere):

// $points is an array full of Point objects (the vertices), which contain lat/long variables
// $ctr is simply a counter that we will test to see if it's even/odd
for ($i = 0, $j = sizeof($points) - 1; $i < sizeof($points); $j = $i++) {
    $p1 = $points[$i];
    $p2 = $points[$j];
    // $test_point is the lat/long pair of the user's address
    if ($p1->lat < $test_point->lat && $p2->lat >= $test_point->lat ||  $p2->lat < $test_point->lat && $p1->lat >= $test_point->lat)  {
        if ($p1->long + ($test_point->lat - $p1->lat)/($p2->lat - $p1->lat)*($p2->long - $p1->long) < $test_point->long)  
            $ctr++;
    }
}

Is there something I'm missing here? I tried deriving a formula on my own and I understand the mathematics behind this to some degree, but is it ok to use GPS coordinates from Google maps with this?

There doesn't seem to be a real pattern as to what is incorrectly reported: I tested for things like addresses close to the boundary or ones in corners of the service area, but no luck there. Also something worth noting is that this service area is just a relatively small region in a city, nothing like state or country-wide areas.

I found this algorithm here quite useful: http://bit.ly/1eJSugg — , Feb 24 '14 at 20:27

JvO · Answer 1 · 2013-01-07T22:49:09.833

2

Well.... your second if() does not compensate for the fact that any of the subtractions may result in a negative number; it would only work if the coordinates are strictly ordered.

Update: On http://rosettacode.org/wiki/Ray-casting_algorithmn there's a whole bunch of algorithms in various languages that describe the process in detail (unfortunately, there is no version for PHP). What seems to be missing from your solution is picking a point that's guaranteed to be outside of the polygon; since you're dealing with longitude/lattitude that should be easy. Second, make sure your polygon is closed (i.e. go from the last point back to the first, if Google Maps doesn't already do so)

edited Jan 07 '13 at 22:49

answered Jan 04 '13 at 03:31

JvO

3,036
2
17
32

lat/lon can have negative values – CSᵠ Jan 04 '13 at 04:00
That's not the point; the formule subtracts various variables without knowing which is the correct order. When dealing with vertices & vectors you have to get the order right or you'll flip what is "inside" or "outside" – JvO Jan 04 '13 at 04:29
Ok... but I can't necessarily take the absolute value of each subtraction result and expect the right answer. How do I know when one is the "right" order? – John V Jan 07 '13 at 03:21

score 0 · Answer 2 · answered Jan 08 '13 at 01:08

Assuming that the $points array contains the corners of a polygon describing the coverage area in clockwise (or counterclockwise order), your code looks correct to me. Basically, it's counting the number of polygon edges that intersect a line drawn due east from the given point to the 180th meridian.

I'd perhaps rewrite it like this, just for clarity:

$p0 = end($points);
foreach ( $points as $p1 ) {
    // ignore edges of constant latitude (yes, this is correct!)
    if ( $p0->lat != $p1->lat ) {
        // scale latitude of $test_point so that $p0 maps to 0 and $p1 to 1:
        $interp = ($test_point->lat - $p0->lat) / ($p1->lat - $p0->lat);
        // does the edge intersect the latitude of $test_point?
        // (note: use >= and < to avoid double-counting exact endpoint hits)
        if ( $interp >= 0 && $interp < 1 ) {
            // longitude of the edge at the latitude of the test point:
            // (could use fancy spherical interpolation here, but for small
            // regions linear interpolation should be fine)
            $long = $interp * $p1->long + (1 - $interp) * $p0->long;
            // is the intersection east of the test point?
            if ( $long < $test_point->long ) {
                // if so, count it:
                $ctr++;
            }
        }
    }
    $p0 = $p1;
}

Note that this code will break in all sorts of interesting ways if the region boundary crosses the 180th meridian, so please don't use it if you have any service regions in the middle of the Pacific Ocean.

If you're still having problems, try plotting the polygon described by the $points array on a map; you may find out that it doesn't look like you thought it would, e.g. if some points are listed in the wrong order.

score 0 · Answer 3 · answered Aug 28 '13 at 06:10

There is a bug with this algorithm, when the ray is tangent to the shape. Just add an epsilon to the test point's latitude when it might happen (line 3 of Ilmari's code):

if ($test_point->lat == $p0->lat)
    $test_point->lat += 0.0000000001;

Also see http://rosettacode.org/wiki/Ray-casting_algorithm (corrected URL).

Thanks.

Raycasting algorithm with GPS coordinates

3 Answers3

Linked