# testing if a Geodesic polygon contains a point c++

I am looking for an algorithm, preferrably in C++, that can take a point (lon,lat) and determine whether that point is inside a geodesic polygon.

Is there such an algorithm freely available?

• Have you tried googling "point in polygon"? Or using the "search" box above? May 18, 2015 at 13:27
• The better approaches to this problem take into account (a) whether the polygon will be fixed or will vary from one call to the next; (b) whether the point will be fixed or can vary between calls; (c) if multiple polygons are involved, whether their interiors can possibly overlap; (d) whether information can be precomputed or has to be computed on the fly; (e) accuracy requirements; and (f) whether a spherical or ellipsoidal earth model will be used. It would be useful to include this information in your question. May 18, 2015 at 15:54

The following is a very close approximation to the right answer which works under the assumption that the vertices of the polygon and the point in question all lie within roughly 2000 km of each other:

Download and install the my C++ library GeographicLib. Using the test point as the center of projection, use the Gnomonic class to transform the vertices of the polygon to the ellipsoidal gnomonic projection. Apply your favorite two-dimensional "point in polygon" test to the resulting polygon.

This works because geodesics transform to very nearly straight lines. To quantify "very nearly", assuming that a 2000 km (resp. 1000 km) bound on the points holds, a straight line in the gnomonic projection deviates from a geodesic by no more than 28 m (resp 1.7 m). For a derivation of the ellipsoidal gnomonic projection and a discussion of the errors, see Section 8 of my paper Algorithms for geodesics.

An obvious extension would be to use the error bound to determine whether to perform a refined test in those cases where the point lies close to an edge. In the spherical limit, the method is exact (but the point and the polygon will need to lie in a single hemisphere for the gnomonic projection to be applied).

Here is code which implements this prescription

``````#include <iostream>
#include <vector>
#include <GeographicLib/Gnomonic.hpp>

struct Point {
double x, y;
Point(double _x, double _y) : x(_x), y(_y) {}
};

// Is p0 inside p?  Polygon
bool inside(const Point p0, const std::vector<Point>& p) {
size_t n = p.size();
bool result = false;
for (size_t i = 0; i < n; ++i) {
size_t j = (i + 1) % n;
if (
// Does p0.y lies in half open y range of edge.
// N.B., horizontal edges never contribute
( (p[j].y <= p0.y && p0.y < p[i].y) ||
(p[i].y <= p0.y && p0.y < p[j].y) ) &&
// is p to the left of edge?
( p0.x < p[j].x + (p[i].x - p[j].x) * (p0.y - p[j].y) /
(p[i].y - p[j].y) )
)
result = !result;
}
return result;
}

int main() {
std::vector<Point> p;
double lat0, lon0;
std::cin >> lat0 >> lon0;
Point p0(0.0, 0.0);
GeographicLib::Gnomonic g;
double lat, lon;
while (std::cin >> lat >> lon) {
double x, y;
g.Forward(lat0, lon0, lat, lon, x, y);
p.push_back(Point(x, y));
}
std::cout << inside(p0, p) << "\n";
}
``````

Suppose points.txt contains the outline of Antarctica:

``````-63.1  -58
-72.9  -74
-71.9 -102
-74.9 -102
-74.3 -131
-77.5 -163
-77.4  163
-71.7  172
-65.9  140
-65.7  113
-66.6   88
-66.9   59
-69.8   25
-70.0   -4
-71.0  -14
-77.3  -33
-77.9  -46
-74.7  -61
``````

Then

``````(echo -74 -50; cat points.txt ) | ./inside
==> 0 (outside)
(echo -74 -20; cat points.txt ) | ./inside
==> 1 (inside)
``````

``````point.within(polygon)