testing if a Geodesic polygon contains a point c++

Question

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?

Answer 1

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).

ADDENDUM 1

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)

ADDENDUM 2

It's just possible that this method is exact even for an ellipsoid because geodesics through the center of project are straight lines in the projection.

Answer 2

In python, there's a library called shapely...one call will do it once you've constructed your polygon and point:

point.within(polygon)

See: https://stackoverflow.com/questions/217578/point-in-polygon-aka-hit-test/2922778#2922778

this was the result of a google search "c++ polygon point within"