Is there a method to tell which hemisphere a point belongs to, in WGS84 geoid, where the equator is an arbitrary ellipsis?

Question

I use the WGS84 geoid (because I'm working with a GPS chip).

I define an ellipsis that splits the surface in two parts, lets name these 'A' and 'B'.

This ellipsis is defined by two points 'a' and 'b' , using their coordinates in WGS84. 'A' is the surface left of the vector 'ab', and therefore 'B' is right of vector 'ab'.

How do I tell if point 'c' belongs to 'A' or not? Any language will do, even pseudo code. I'll figure out how to translate to ATMega assembly later...

Many thanks.

PS: What I currently think of is switching origins so that my ellipsis is the new equator, then applying the transform on my point and finding out if it's North or South. But whereas I could find out how to do it in polar coordinates systems, I'm totally lost in geoids... maybe its not even possible?

PPS: Here is an example:

The points are:

• 50°21'37.9"N 3°05’01.8"E : A
• 50°21'36.3"N 3°05'04.2"E : B
• 50°21'40.2"N 3°05'05.9"E : C1, belongs to A
• 50°21'32.1"N 3°04'52.7"E : C2, does not belong to A

Answer 1

Solution

Using the analysis given at https://gis.stackexchange.com/a/20250 we may find the earth-centered Cartesian coordinates for a point at geodetic latitude f and longitude l are

(x, y, z) = (a*cos(l)*cos(t), a*sin(l)*cos(t), b*sin(t))

where a is the semi-major axis (6 378 137 meters on the WGS 84 ellipsoid), b is the semi-minor axis (approximately 6 356 752.3142 meters on the WGS 84 ellipsoid), and t is the geocentric latitude defined by the equation

tan(t) = (b/a) * tan(f).

Three points on the ellipsoid will determine three triples of geocentric coordinates P, Q, and R. The side of the PQ plane on which R lies is determined by the sign of the determinant

|P;Q;R|

(where P;Q;R is the 3 by 3 matrix created by stacking the three vectors in that given order). Which side? The positive side is the one in which (P, Q, R) form a right-handed basis. On a correctly oriented map, when the distance between P and Q is small, the positive side is to the left of the vector from P to Q.

The result is indeterminate when P is parallel to Q; that is, when those two points are diametrically opposite.

---

Note that this solution does not determine which side of the geodesic between P and Q the point R is on: it determines which side of the 3D plane (determined by P, Q, and the earth's center) the point R is on. I have chosen this particular interpretation because the question refers to "ellipsis" (which I read as meaning ellipse). In general, geodesics on an ellipsoid are not ellipses--they aren't even closed loops!--whereas the intersection of a plane (passing through the center) with an ellipsoid always is an ellipse.

---

Worked example

Consider the points given in the question:

P: 50°21'37.9"N 3°05’01.8"E 
Q: 50°21'36.3"N 3°05'04.2"E 
R1: 50°21'40.2"N 3°05'05.9"E 
R2: 50°21'32.1"N 3°04'52.7"E

In degrees these coordinates (the geodetic latitude and longitude) are

           Lat        Lon Geocentric Lat
P:  50.3605278  3.0838333     50.2659655
Q:  50.3600833  3.0845000     50.2655208 
R1: 50.3611667  3.0849722     50.2666048
R2: 50.3589167  3.0813056     50.2643534

Their corresponding geocentric latitudes are listed next to them. They look about right, because according to the figure at https://gis.stackexchange.com/a/108220, they should be a little bit closer to zero than the geodetic latitudes. With these values (and their longitudes) the geocentric Cartesian coordinates work out to

          X       Y       Z 
P:  4071159  219334 4888470 
Q:  4071194  219383 4888438 
R1: 4071100  219412 4888515
R2: 4071306  219162 4888355

The matrix P;Q;R1 is given by the first three lines of the preceding tableau. Its determinant, 46963622641, is positive: R1 therefore lies on the positive side of PQ (which is to the left of the vector P-->Q on the map). The determinant of P;Q;R2, -110723586947, is negative: R2 therefore lies to the right of P-->Q.

---

R code to carry out these calculations is offered to serve as working pseudocode for any implementation. Data are represented as 2 by 3 arrays: the first row is the geodetic latitude, the second row the longitude; and the columns give the degrees, minutes, and seconds as numbers. The output of side is either +1 (lies to the left), -1 (lies to the right), or 0 (lies on the plane).

point.make <- function(d.lat,m.lat,s.lat, d.lon,m.lon,s.lon) {
  p <- rbind(c(d.lat,m.lat,s.lat), c(d.lon,m.lon,s.lon))
  colnames(p) <- c("Deg", "Min", "Sec")
  rownames(p) <- c("Lat", "Lon")
  return (p) # A 2 by 3 array
}
to.3D <- function(x, a=6378137, b=6356752.3142) {
  to.degrees <- function(x) (((x[3]/60) + x[2])/60 + abs(x[1])) * sign(x[1])
  to.radians <- function(x) x / 180 * pi
  u <- to.radians(apply(x, 1, to.degrees))                 # Geodetic Lat, lon
  t <- atan2(b * sin(u[1]), a * cos(u[1]))                 # Geocentric latitude
  p <- c(a * cos(t) * c(cos(u[2]), sin(u[2])), b * sin(t)) # Geocentric Cartesian coords
  names(p) <- c("X", "Y", "Z")
  return (p) # A 3-vector.
}
side <- function(r, p, q) sign(det(sapply(list(p, q, r), to.3D))) # +1, -1, or 0
#
# Data.
#
p <- point.make(50,21,37.9, 3,05,01.8)
q <- point.make(50,21,36.3, 3,05,04.2)
r1 <- point.make(50,21,40.2, 3,05,05.9)
r2 <- point.make(50,21,32.1, 3,04,52.7)
#
# Results.
#
side(r1, p, q) # Value is 1
side(r2, p, q) # Value is -1