Preface: I'm unsure if GIS or math is the correct community for this question.
I have a set of 9 WGS84 latitude/longitude (degree) pairs for which I need to find the weighted average of.
To prevent things like the dateline causing massively warped averages, I'm converting from WGS84 into ECEF Cartesian (X, Y, Z), weighted averaging those coordinates, then converting back to a WGS84 latitude/longitude (degree) pair.
The problem is that this produces an incorrect result - the latitude seems to 'drift' the further from the equator you are.
The conversion algorithms used (from WGS84 to ECEF, then ECEF to WGS84) work as intended and produce accurate results. The problem appears to lie with the weighted averaging of the Cartesian coordinates - I'm not sure why this is.
Question: How do I produce a weighted average of a set of Cartesian coordinates that, when converted back to WGS84, will not cause the latitude to drift?
Notes:
Naively weighted averaging unconverted WGS84 latitude/longitude pairs produces the correct result (unless the set crosses the dateline, for example).
As a test, I converted the naively weighted averaged WGS84 lat/long pairs to and from ECEF, and this produced the correct result, which means it's definitely something to do with weighted averaging the ECEF coordinates and not the conversion algorithms used.
Sample code (in Python):
WGS84 -> ECEF:
lat = np.radians(lat)
lon = np.radians(lon)
d = e * np.sin(lat)
N = a / np.sqrt(1. - d * d)
x = (N + ele) * np.cos(lat) * np.cos(lon)
y = (N + ele) * np.cos(lat) * np.sin(lon)
z = ((1 - esq) * N + ele) * np.sin(lat)
return x, y, z
ECEF -> WGS84:
r = np.sqrt(x * x + y * y)
Esq = a * a - b * b
F = 54 * b * b * z * z
G = r * r + (1 - esq) * z * z - esq * Esq
C = (esq * esq * F * r * r) / np.power(G, 3)
S = cbrt(1 + C + np.sqrt(C * C + 2 * C))
P = F / (3 * np.power((S + 1 / S + 1), 2) * G * G)
Q = np.sqrt(1 + 2 * esq * esq * P)
r_0 = -(P * esq * r) / (1 + Q) + np.sqrt(0.5 * a * a * (1 + 1.0 / Q) - \
P * (1 - esq) * z * z / (Q * (1 + Q)) - 0.5 * P * r * r)
U = np.sqrt(np.power((r - esq * r_0), 2) + z * z)
V = np.sqrt(np.power((r - esq * r_0), 2) + (1 - esq) * z * z)
Z_0 = b * b * z / (a * V)
h = U * (1 - b * b / (a * V))
lat = np.arctan((z + e1sq * Z_0) / r)
lon = np.arctan2(y, x)
xi = np.sqrt(1 - esq * np.sin(lat) * np.sin(lat))
N = a / xi
ele = (r / np.cos(np.radians(lat))) - N
return np.rad2deg(lat), np.rad2deg(lon), ele
Weighted average:
x = np.average(cart_x, weights)
y = np.average(cart_y, weights)
z = np.average(cart_z, weights)
(where cart_x, cart_y, cart_z are lists of x, y, and z coordinates in same order as their corresponding weights)
Sample data:
Input (WGS84 latitude, WGS84 longitude, Weighting):
53.708794718 -0.876990985207 0.0
56.7088004625 -0.876990985207 0.125
53.9512592773 -0.876990985207 0.125
53.708794718 2.12301475927 0.125
53.708794718 -0.634526425941 0.125
50.7087889735 -0.876990985207 0.125
53.4663301588 -0.876990985207 0.125
53.708794718 -3.87699672969 0.125
53.708794718 -1.11945554447 0.125
...converted to X, Y, Z:
3782945.06909 -57907.7178905 5117626.09215
3508729.39335 -53710.140685 5308169.95693
3761162.0512 -57574.272696 5133553.1923
3780791.31608 140156.069292 5117626.09215
3783156.24975 -41898.5617628 5117626.09215
4046681.01116 -61944.87578 4913093.31984
3804659.61535 -58240.1149506 5101607.59494
3774729.99291 -255812.783436 5117626.09215
3782666.14289 -73915.8369984 5117626.09215
Weighted average of Cartesian coordinates:
3780321.97159 -57867.5646271 5115866.05408
Which produces the incorrect result of:
53.7184377538 -0.876990985207
Expected result of (approximately):
53.708794718 -0.876990985207
Difference: 0.00964303580069 0.0
