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