2

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? – Vince May 18 '15 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. – whuber May 18 '15 at 15:54
2

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.

-1

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"

  • The question specified a C++ implementation of an algorithm. You could salvage this answer with a link to the source which implements within. – Vince May 19 '15 at 5:03
  • 1
    @Vince fair comment. – user1269942 May 19 '15 at 16:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.