Detecting whether a GPS coordinate falls within a polygon on a map

Question

As stated in the title, the goal is to have a way for detecting whether a given GPS coordinate falls inside a polygon or not.

The polygon itself can be either convex or concave. It's defined as a set of edge vectors and a known point within that polygon. Each edge vector is further defined by four coordinates which are the latitudes and longitudes of respective tip points and a bearing relative to the starting point.

There are a couple of questions similar to this one here on StackOverflow but they describe the solution only in general terms and for a 2D plane, whereas I am looking for an existing implementation that supports polygons defined by latitude/longitude pairs in WGS 84.

What API-s or services are out there for doing such collision tests?

Answer 1

Here is a java program which uses a function that will return true if a latitude/longitude is found inside of a polygon defined by a list of lat/longs, with demonstration for the state of florida.

I'm not sure if it deals with the fact that the lat/long GPS system is not an x/y coordinate plane. For my uses I have demonstrated that it works (I think if you specify enough points in the bounding box, it washes away the effect that the earth is a sphere, and that straight lines between two points on the earth is not an arrow straight line.

First specify the points that make up the corner points of the polygon, it can have concave and convex corners. The coordinates I use below traces the perimeter of the state of Florida.

method coordinate_is_inside_polygon utilizes an algorithm I don't quite understand. Here is an official explanation from the source where I got it:

"... solution forwarded by Philippe Reverdy is to compute the sum of the angles made between the test point and each pair of points making up the polygon. If this sum is 2pi then the point is an interior point, if 0 then the point is an exterior point. This also works for polygons with holes given the polygon is defined with a path made up of coincident edges into and out of the hole as is common practice in many CAD packages. "

My unit tests show it does work reliably, even when the bounding box is a 'C' shape or even shaped like a Torus. (My unit tests test many points inside Florida and make sure the function returns true. And I pick a number of coordinates everywhere else in the world and make sure it returns false. I pick places all over the world which might confuse it.

I'm not sure this will work if the polygon bounding box crosses the equator, prime meridian, or any area where the coordinates change from -180 -> 180, -90 -> 90. Or your polygon wraps around the earth around the north/south poles. For me, I only need it to work for the perimeter of Florida. If you have to define a polygon that spans the earth or crosses these lines, you could work around it by making two polygons, one representing the area on one side of the meridian and one representing the area on the other side and testing if your point is in either of those points.

Here is where I found this algorithm: Determining if a point lies on the interior of a polygon - Solution 2

Run it for yourself to double check it.

Put this in a file called Runner.java

import java.util.ArrayList; public class Runner {     public static double PI = 3.14159265;     public static double TWOPI = 2*PI;     public static void main(String[] args) {     ArrayList<Double> lat_array = new ArrayList<Double>();     ArrayList<Double> long_array = new ArrayList<Double>();      //This is the polygon bounding box, if you plot it,      //you'll notice it is a rough tracing of the parameter of      //the state of Florida starting at the upper left, moving      //clockwise, and finishing at the upper left corner of florida.      ArrayList<String> polygon_lat_long_pairs = new ArrayList<String>();     polygon_lat_long_pairs.add("31.000213,-87.584839");       //lat/long of upper left tip of florida.     polygon_lat_long_pairs.add("31.009629,-85.003052");     polygon_lat_long_pairs.add("30.726726,-84.838257");     polygon_lat_long_pairs.add("30.584962,-82.168579");     polygon_lat_long_pairs.add("30.73617,-81.476441");       //lat/long of upper right tip of florida.     polygon_lat_long_pairs.add("29.002375,-80.795288");     polygon_lat_long_pairs.add("26.896598,-79.938355");     polygon_lat_long_pairs.add("25.813738,-80.059204");     polygon_lat_long_pairs.add("24.93028,-80.454712");     polygon_lat_long_pairs.add("24.401135,-81.817017");     polygon_lat_long_pairs.add("24.700927,-81.959839");     polygon_lat_long_pairs.add("24.950203,-81.124878");     polygon_lat_long_pairs.add("26.0015,-82.014771");     polygon_lat_long_pairs.add("27.833247,-83.014527");     polygon_lat_long_pairs.add("28.8389,-82.871704");     polygon_lat_long_pairs.add("29.987293,-84.091187");     polygon_lat_long_pairs.add("29.539053,-85.134888");     polygon_lat_long_pairs.add("30.272352,-86.47522");     polygon_lat_long_pairs.add("30.281839,-87.628784");      //Convert the strings to doubles.            for(String s : polygon_lat_long_pairs){         lat_array.add(Double.parseDouble(s.split(",")[0]));         long_array.add(Double.parseDouble(s.split(",")[1]));     }     //prints TRUE true because the lat/long passed in is     //inside the bounding box.     System.out.println(coordinate_is_inside_polygon(             25.7814014D,-80.186969D,             lat_array, long_array));      //prints FALSE because the lat/long passed in      //is Not inside the bounding box.     System.out.println(coordinate_is_inside_polygon(             25.831538D,-1.069338D,              lat_array, long_array));  } public static boolean coordinate_is_inside_polygon(     double latitude, double longitude,      ArrayList<Double> lat_array, ArrayList<Double> long_array) {               int i;        double angle=0;        double point1_lat;        double point1_long;        double point2_lat;        double point2_long;        int n = lat_array.size();         for (i=0;i<n;i++) {           point1_lat = lat_array.get(i) - latitude;           point1_long = long_array.get(i) - longitude;           point2_lat = lat_array.get((i+1)%n) - latitude;            //you should have paid more attention in high school geometry.           point2_long = long_array.get((i+1)%n) - longitude;           angle += Angle2D(point1_lat,point1_long,point2_lat,point2_long);        }         if (Math.abs(angle) < PI)           return false;        else           return true; }  public static double Angle2D(double y1, double x1, double y2, double x2) {    double dtheta,theta1,theta2;     theta1 = Math.atan2(y1,x1);    theta2 = Math.atan2(y2,x2);    dtheta = theta2 - theta1;    while (dtheta > PI)       dtheta -= TWOPI;    while (dtheta < -PI)       dtheta += TWOPI;     return(dtheta); }  public static boolean is_valid_gps_coordinate(double latitude,      double longitude) {     //This is a bonus function, it's unused, to reject invalid lat/longs.     if (latitude > -90 && latitude < 90 &&              longitude > -180 && longitude < 180)     {         return true;     }     return false; } } 

Demon magic needs to be unit-tested. Put this in a file called MainTest.java to verify it works for you

import java.util.ArrayList; import org.junit.Test; import static org.junit.Assert.*;  public class MainTest { @Test public void test_lat_long_in_bounds(){     Runner r = new Runner();     //These make sure the lat/long passed in is a valid gps      //lat/long coordinate.  These should be valid.      assertTrue(r.is_valid_gps_coordinate(25, -82));     assertTrue(r.is_valid_gps_coordinate(-25, -82));     assertTrue(r.is_valid_gps_coordinate(25, 82));     assertTrue(r.is_valid_gps_coordinate(-25, 82));     assertTrue(r.is_valid_gps_coordinate(0, 0));     assertTrue(r.is_valid_gps_coordinate(89, 179));     assertTrue(r.is_valid_gps_coordinate(-89, -179));     assertTrue(r.is_valid_gps_coordinate(89.999, 179));     //If your bounding box crosses the equator or prime meridian,      then you have to test for those situations still work. } @Test public void realTest_for_points_inside() {     ArrayList<Double> lat_array = new ArrayList<Double>();     ArrayList<Double> long_array = new ArrayList<Double>();     ArrayList<String> polygon_lat_long_pairs = new ArrayList<String>();     //upper left tip of florida.     polygon_lat_long_pairs.add("31.000213,-87.584839");     polygon_lat_long_pairs.add("31.009629,-85.003052");     polygon_lat_long_pairs.add("30.726726,-84.838257");     polygon_lat_long_pairs.add("30.584962,-82.168579");     polygon_lat_long_pairs.add("30.73617,-81.476441");       //upper right tip of florida.     polygon_lat_long_pairs.add("29.002375,-80.795288");     polygon_lat_long_pairs.add("26.896598,-79.938355");     polygon_lat_long_pairs.add("25.813738,-80.059204");     polygon_lat_long_pairs.add("24.93028,-80.454712");     polygon_lat_long_pairs.add("24.401135,-81.817017");     polygon_lat_long_pairs.add("24.700927,-81.959839");     polygon_lat_long_pairs.add("24.950203,-81.124878");     polygon_lat_long_pairs.add("26.0015,-82.014771");     polygon_lat_long_pairs.add("27.833247,-83.014527");     polygon_lat_long_pairs.add("28.8389,-82.871704");     polygon_lat_long_pairs.add("29.987293,-84.091187");     polygon_lat_long_pairs.add("29.539053,-85.134888");     polygon_lat_long_pairs.add("30.272352,-86.47522");     polygon_lat_long_pairs.add("30.281839,-87.628784");      for(String s : polygon_lat_long_pairs){         lat_array.add(Double.parseDouble(s.split(",")[0]));         long_array.add(Double.parseDouble(s.split(",")[1]));     }      Runner r = new Runner();     ArrayList<String> pointsInside = new ArrayList<String>();     pointsInside.add("30.82112,-87.255249");     pointsInside.add("30.499804,-86.8927");     pointsInside.add("29.96826,-85.036011");     pointsInside.add("30.490338,-83.981323");     pointsInside.add("29.825395,-83.344116");     pointsInside.add("30.215406,-81.828003");     pointsInside.add("29.299813,-82.728882");     pointsInside.add("28.540135,-81.212769");     pointsInside.add("27.92065,-82.619019");     pointsInside.add("28.143691,-81.740113");     pointsInside.add("27.473186,-80.718384");     pointsInside.add("26.769154,-81.729126");     pointsInside.add("25.853292,-80.223999");     pointsInside.add("25.278477,-80.707398");     pointsInside.add("24.571105,-81.762085");   //bottom tip of keywest     pointsInside.add("24.900388,-80.663452");     pointsInside.add("24.680963,-81.366577");      for(String s : pointsInside)     {         assertTrue(r.coordinate_is_inside_polygon(             Double.parseDouble(s.split(",")[0]),              Double.parseDouble(s.split(",")[1]),              lat_array, long_array));     } }  @Test public void realTest_for_points_outside() {     ArrayList<Double> lat_array = new ArrayList<Double>();     ArrayList<Double> long_array = new ArrayList<Double>();      ArrayList<String> polygon_lat_long_pairs = new ArrayList<String>();     //upper left tip, florida.     polygon_lat_long_pairs.add("31.000213,-87.584839");     polygon_lat_long_pairs.add("31.009629,-85.003052");     polygon_lat_long_pairs.add("30.726726,-84.838257");     polygon_lat_long_pairs.add("30.584962,-82.168579");     polygon_lat_long_pairs.add("30.73617,-81.476441");     //upper right tip, florida.     polygon_lat_long_pairs.add("29.002375,-80.795288");     polygon_lat_long_pairs.add("26.896598,-79.938355");     polygon_lat_long_pairs.add("25.813738,-80.059204");     polygon_lat_long_pairs.add("24.93028,-80.454712");     polygon_lat_long_pairs.add("24.401135,-81.817017");     polygon_lat_long_pairs.add("24.700927,-81.959839");     polygon_lat_long_pairs.add("24.950203,-81.124878");     polygon_lat_long_pairs.add("26.0015,-82.014771");     polygon_lat_long_pairs.add("27.833247,-83.014527");     polygon_lat_long_pairs.add("28.8389,-82.871704");     polygon_lat_long_pairs.add("29.987293,-84.091187");     polygon_lat_long_pairs.add("29.539053,-85.134888");     polygon_lat_long_pairs.add("30.272352,-86.47522");     polygon_lat_long_pairs.add("30.281839,-87.628784");      for(String s : polygon_lat_long_pairs)     {         lat_array.add(Double.parseDouble(s.split(",")[0]));         long_array.add(Double.parseDouble(s.split(",")[1]));     }      Runner r = new Runner();      ArrayList<String> pointsOutside = new ArrayList<String>();     pointsOutside.add("31.451159,-87.958374");     pointsOutside.add("31.319856,-84.607544");     pointsOutside.add("30.868282,-84.717407");     pointsOutside.add("31.338624,-81.685181");     pointsOutside.add("29.452991,-80.498657");     pointsOutside.add("26.935783,-79.487915");     pointsOutside.add("25.159207,-79.916382");     pointsOutside.add("24.311058,-81.17981");     pointsOutside.add("25.149263,-81.838989");     pointsOutside.add("27.726326,-83.695679");     pointsOutside.add("29.787263,-87.024536");     pointsOutside.add("29.205877,-62.102052");     pointsOutside.add("14.025751,-80.690919");     pointsOutside.add("29.029276,-90.805666");     pointsOutside.add("-12.606032,-70.151369");     pointsOutside.add("-56.520716,-172.822269");     pointsOutside.add("-75.89666,9.082024");     pointsOutside.add("-24.078567,142.675774");     pointsOutside.add("84.940737,177.480462");     pointsOutside.add("47.374545,9.082024");     pointsOutside.add("25.831538,-1.069338");     pointsOutside.add("0,0");      for(String s : pointsOutside){         assertFalse(r.coordinate_is_inside_polygon(             Double.parseDouble(s.split(",")[0]),             Double.parseDouble(s.split(",")[1]), lat_array, long_array));     } } } //The list of lat/long inside florida bounding box all return true. //The list of lat/long outside florida bounding box all return false. 

I used eclipse IDE to get this to run java using java 1.6.0. For me all the unit tests pass. You need to include the junit 4 jar file in your classpath or import it into Eclipse.

Answer 2

I thought similarly as shab first (his proposal is called Ray-Casting Algorithm), but had second thoughts like Spacedman:

...but all the geometry will have to be redone in spherical coordinates...

I implemented and tested the mathematically correct way of doing that, e.i. intersecting great circles and determining whether one of the two intersecting points is on both arcs. (Note: I followed the steps described here, but I found several errors: The sign function is missing at the end of step 6 (just before arcsin), and the final test is numerical garbage (as subtraction is badly conditioned); use rather L_1T >= max(L_1a, L_1b) to test whether S1 is on the first arc etc.)

That also is extremely slow and a numerical nightmare (evaluates ~100 trigonometric functions, among other things); it proved not to be usable in our embedded systems.

There's a trick, though: If the area you are considering is small enough, just do a standard cartographic projection, e.g. spherical Mercator projection, of each point:

// latitude, longitude in radians x = longitude; y = log(tan(pi/4 + latitude/2)); 

Then, you can apply ray-casting, where the intersection of arcs is checked by this function:

public bool ArcsIntersecting(double x1, double y1, double x2, double y2,    double x3, double y3, double x4, double y4)     {      double vx1 = x2 - x1;     double vy1 = y2 - y1;      double vx2 = x4 - x3;     double vy2 = y4 - y3;      double denom = vx1 * vy2 - vx2 * vy1;      if (denom == 0) { return false; } // edges are parallel      double t1 = (vx2 * (y1 - y3) - vy2 * (x1 - x3)) / denom;      double t2;      if (vx2 != 0) { t2 = (x1 - x3 + t1 * vx1) / vx2; }     else if (vy2 != 0) { t2 = (y1 - y3 + t1 * vy1) / vy2; }     else { return false; } // edges are matching      return min(t1, t2) >= 0 && max(t1, t2) <= 1; }