Checking if a longitude/latitude coordinate resides inside a complex polygon in an embedded device?

Question

I need the user to be able to draw a complex polygon on a map and then have the application check if a given longitude/latitude resides within that polygon.

I was only able to find algorithms that were using a simple x/y cartesian coordinate system that doesn't compensate for the curvature of the earth.

The user draws the polygon on a PC, where the points are transferred over radio to a embedded device, which then needs to check if the given polygon resides within it's current position (taken from GPS).

As this is for an embedded device I am not able to use huge libraries, rather I need the algorithm to perform the check myself or a very small library. But I seem to be unable to find any such algorithm.

Solution 1

Here's an implementation I wrote in C# for a Polygon class that contains a list of vertices. It doesn't consider the curvature of the Earth. Rather, you would pre-process the polygon into smaller segments prior to running this.

The performance of this algorithm is very good. Even for polygons with thousands of edges it completes in around one or two milliseconds on my desktop.

The code has been optimised quite a bit and so isn't that readable as psuedo-code.

public bool Contains(GeoLocation location)
{
    if (!Bounds.Contains(location))
        return false;

    var lastPoint = _vertices[_vertices.Length - 1];
    var isInside = false;
    var x = location.Longitude;
    foreach (var point in _vertices)
    {
        var x1 = lastPoint.Longitude;
        var x2 = point.Longitude;
        var dx = x2 - x1;

        if (Math.Abs(dx) > 180.0)
        {
            // we have, most likely, just jumped the dateline (could do further validation to this effect if needed).  normalise the numbers.
            if (x > 0)
            {
                while (x1 < 0)
                    x1 += 360;
                while (x2 < 0)
                    x2 += 360;
            }
            else
            {
                while (x1 > 0)
                    x1 -= 360;
                while (x2 > 0)
                    x2 -= 360;
            }
            dx = x2 - x1;
        }

        if ((x1 <= x && x2 > x) || (x1 >= x && x2 < x))
        {
            var grad = (point.Latitude - lastPoint.Latitude) / dx;
            var intersectAtLat = lastPoint.Latitude + ((x - x1) * grad);

            if (intersectAtLat > location.Latitude)
                isInside = !isInside;
        }
        lastPoint = point;
    }

    return isInside;
}

The basic idea is to find all edges of the polygon that span the 'x' position of the point you're testing against. Then you find how many of them intersect the vertical line that extends above your point. If an even number cross above the point, then you're outside the polygon. If an odd number cross above, then you're inside.

Solution 2

Good explanations and simple C code that you can convert to your needs

http://alienryderflex.com/polygon/

Combine the polygon check with a RTree for quick culling of the search tree if you have many non overlapping polygons.

Comments

[PeterJ]

What sort of distances are you looking at? I do this sort of work a lot and if you're talking about relatively small polygons like say < 1000 meters those simple algorithms will work fine relative to the accuracy of a standard GPS. Earth curvature doesn't have much effect over smaller distances.

[Drew Noakes]

Note that if you're working with GIS, you might like this SO sister-site: gis.stackexchange.com

[ChewToy]

@PeterJ I'm looking to use this both on a very-wide scale (hundreds of kilometers) but still require very-close precision as it will also be used within segmented areas that are < 1 km

[PeterJ]

Fair enough, but remember unless segments are a long way apart it won't make a difference. For example a 200KM highway with a series of segments a few KM apart won't make much difference. It's only if you say have an edge on a polygon 100KM long that it's really going to make a difference, which is pretty rare in my experience. Most real-life things like say a city border aren't that straight in practice so it's worth some thought about likely errors and if you really need it. At highway speeds you'll probably also need some dead reckoning too if you really need that precision.

[ChewToy]

Thanks for your reply! What kind of precision have you been able to achieve using this algorithm? Considering it doesn't compensate for the curvature of the earth, does it only work with semi-small distances or will it work decently even with larger polygons (say 50km or so)?

[Drew Noakes]

@ChewToy, precision varies with distance. If you need more precision that it gives you, find the midpoint of two adjacent vertices on the sphere, and insert a new vertex. Do this repeatedly until the polygon is within your required tolerances. Alternatively, modify the algorithm to perform the intersection test in a way that accounts for the curvature of the earth.

[Admin]

Many thanks for the algorithm, what does the condition Bounds.Contains(location) do?

[Drew Noakes]

@VladimirGhetau that is an optional optimisation. If you have a rectangular bounding box that contains the extent of all points in the polygon then it might be faster to do a quick test on the box before walking all the vertices. If the point is not within the bounds then it cannot be within the polygon.

[Admin]

Hi Drew, this algorithm will tell you if a random vertex is within a polygon. I am trying to exclude the case when the random vertex is actually a vertex forming the polygon by removing any >= and <= comparison and turning it into > and <, any advice on this?

[Admin]

Just figured, what I need is the Ray Casting algorithm, this algorithm doesn't return True if the point I'm looking for matches a vertex.

[Mohamad Mousheimish]

Can you please explain Bounds.Contains(location) this syntax. I can't find any class in C# called Bounds. And since I didn't find it, I've removed it from the method. It's not working as it should of course.

[Drew Noakes]

@MohamadMousheimish Bounds and GeoLocation types are not included. You can adapt the algorithm to your own types for these kinds of things.