I spent a while looking for a robust algorithm to do this on the sphere for time zone look up and didn't even find good pseudocode, let alone c / c++ I was entirely happy with. I'm going to run through what I found. I'll finish with resources that could be used to put this together fairly easily.
It's Easy …On the Plane
The problem is known as "Point in Polygon".
A frequently used and simple POP algorithm is "Ray Casting". POP on the 2D plane relies on a point at infinity. On the plane this is very easy. There are an infinite number of points at infinity. Pick any one! But there is no such point on the sphere.
You can get away with this if you have a known point either inside or outside of any given query polygon. Given your use case, this is not an onerous requirement: you can easily pick any single point in the sea, and this will be outside of all of the land polygons.
The "Winding Number" POP algorithm also fails (as far as I can see) because on the sphere you can approach any edge in either of two directions.
An Algorithm
I wanted an approach that would work without an auxiliary point and without heuristics (used to generate an auxiliary point from the edge data). If I'm honest, I wanted this because I was convinced it should be possible, no because I genuinely needed it.
For your use case, you can get away with the usual ray casting algorithm and a single point known to be in the ocean, so you don't need to rely on heuristics, although they probably work fairly well anyway.
The approach that I came up with works like this…
- You'll need to loop through your polygons yourself. For each polygon…
- Find a great circle through your query point and through at least one edge of your polygon (a mid point between two corners will do).
- Intersect the great circle with your polygon's edges.
- The inside of the polygon is to your right as you walk its edges in sequence (or to the left, if you prefer). This gives each intersection point sufficient information to know on which side is the inside or outside.
- From the nearest intersection point, you can determine if your query point is inside or outside.
Implementation Tips
If you plan on implementing this (or any of the other POP algorithms), don't try to work in sine or cosine.
Represent your points (polygon corners & query point) as unit vectors. Represent your great circles (polygon edges and the circle the query point is on) as unit vectors normal to the plane that the great is on. Use the dot product and cross product. Don't think in angles. Think in vectors.
It shouldn't be too hard – I didn't need it quite enough to implement it. Get in contact if you'd like it written freelance!
Links From Which You Could Build a Solution
The c++ boost library has a POP implementation that I also didn't like but this is largely because I'm a perfectionist – I imagine it will serve the purpose in nearly all cases.
The tz_world database contains polygons for land masses, and there is a GeoJSON variant of it. You could parse this nicely with the built in NSJSONSerialization class.
Here are some algorithms by NASA for points and spheres (I didn't like their POP though).
