Longitude at which a "line" of constant latitude cuts the "line" between two coordinates

[RobMeades]

Hi there, and thanks for this most excellent code, all nicely in C, an increasing rarity. Not an issue, just a question:

I am implementing a geofence for addition to an [Apache licensed] C library, https://github.com/u-blox/ubxlib. I intend to use ray casting to resolve the "is my point inside a polygon?" question, which requires me to work out the point at which a line of constant latitude cuts a line between two points (that are adjacent vertices of the polygon).

For small geofences, less than 1 kmish, I will do this with plain old linear maths, but I would like to offer a solution that takes the curvature of the earth into account so that geofences of any size may be used.

I can see that the C++ version of this library offers such an intersect function: is there a way to do the same calculation using the C version?

[cffk]

At this point, I'm not planning on porting the intersection capability to the C library. Your choices then are:

• Do the port yourself.
• Write a C wrapper to the C++ class.
• Implement just Intersect::Basic and have the segment intersection use this. This would be right most of the time.

The last might be a reasonable starting point.

HOWEVER, a circle of latitude isn't a geodesic! So you would have to reformulate how you solve this problem. Perhaps you could switch to using meridians?

In either case, I see you running into problems at the poles. I realize that this may not be your typical use case, but still...

If I wanted to tackle this problem and was willing to limit myself to polygons with extent less than 1000 km, say, I would just use my gnomonic projection with the point in question being the center of the projection. (And the gnomonic projection is easy to implement in terms of geodesics.)

[RobMeades]

Hi, and thanks very much for the swift response. Thinking about it, there is no reason why we should not require the use of g++ for this case: all of our code compiles happily enough when consumed by g++ so there should be no issue with making a wrapper.

On the "latitude isn't a geodesic" that is not a subtlety that I'd appreciated until now: not a problem to switch to using meridians, it doesn't matter which direction our ray is heading. Concerning behaviour near the poles, I don't believe that is an issue for the use-cases of our customers; I can add a note about it as a warning.

I must admit that I'm hanging on by my fingernails with this stuff; as you can see I'm barely up to speed with geodesics :-). However I will spend a little longer parsing your suggestion concerning the use of gnomonic projection as a way forward.

[cffk]

Yes, thinking some more about your problem, using the gnomonic projection is the way to go:

• It transforms with ellipsoidal problem into a planar problem.
• The implementation is very simple (see src/Gnomonic.cpp in the C++ distribution).
• Its properties are given in Sec. 8 of Algorithms for geodesics.
• Finally, as long as all the points of the polygon are within 10000 km of the test point, the inside/outside test is exact.

[RobMeades]

The thing I'm having trouble with here is relating the real world to the abstract concept of the Point In Polygon solution, ray casting version.

For the ray casting solution all that matters [I think] is that a ray emitted from the point ends up outside the polygon: the ray can be any shape, as squirly as you like, an infinite spiral or scribble, and can leave/re-enter the polygon any number of times, it just has to end up outside the polygon for the counting to work.

Given that, shouldn't any circle around the earth work, provided it is indeed a circle? Bare in mind that for a geofence I don't care about the distance to the polygon (though that might be a nice piece of additional information if it were not computationally intensive to get).

In fact, do I need to even worry about the earth? Provided the ray I am projecting has no end, does that [maybe] meet the requirement?

[Feel free to close this issue, I don't mean to take up your time unnecessarily with this].

[cffk]

With the gnomonic projection, the points at infinity correspond roughly to a circle of radius 10000 km from the center point. So the point-in-polygon algorithm works provided the polygon is within this circle. The key property of the gnomonic projection is that geodesics project to (nearly) straight lines. An arbitrary mapping of points on the earth to a plane does not have this property (and may, for example, project a circle of latitude to a straight line).

[RobMeades]

Understood, and thanks again for your swift response, that makes clear the unique properties of gnomonic projections. I need to get my head around this, let me sleep on it: I need to resolve in my mind the instances where an intersection with any arbitrary infinite line over the surface would fail the in/out test while the gnomonic projection would succeed.

I think my confusion is that I'm concentrating on the test of inside/outside-ness while that is the trivial part; what actually matter is is the point inside or outside in the first place, and that is why the calculations need to be WGS84 aware, which is [part of] what the gnomonic projection does for me.

[RobMeades]

To test my understanding (he said, failing to sleep), looking at Gnomonic: Forward(), I would call it with the lat/lon of my point as lat0/lon0 and the lat/lon of one end of a side of the polygon, let's call it (a):

Forward(pointLat, pointLon, aVertexLat, aVertexLon, &aVertexX, &aVertexY, NULL, NULL);

I would repeat for the lat/lon of the other end of the side, (b):

Forward(pointLat, pointLon, bVertexLat, bVertexLon, &bVertexX, &bVertexY, NULL, NULL);

I now have everything in an X/Y space, with pointLat/pointLon at the origin, and so calculating the intersection with the side is good 'ole geometry. Is that correct?

[cffk]

Yes, that's it. Except that your C++ invocation is incorrect (don't use the address operator & for a reference argument) and you should use the overloaded function

Forward(pointLat, pointLon, aVertexLat, aVertexLon, aVertexX, aVertexY);
Forward(pointLat, pointLon, bVertexLat, bVertexLon, bVertexX, bVertexY);

[RobMeades]

Excellent: I can resolve the C++ mechanics (haven't used C++ in years, spend most of my time in resource-constrained real-time embedded stuff) it was the principle I hadn't quite got my head around. Now I can sleep :-).

Thank you so much for your help with this, my eyes have been opened; we can sub-module your library and all is, hopefully, straight (for some value of straight) with the world.