Averaging Cartesian ECEF coordinates then converting to WGS84 without latitude drift

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

Answer

The comment from @mkennedy revealed you may be missing one important part of your weighted average of Cartesian coordinates:

3780321.97159 -57867.5646271 5115866.05408

Converting these ECEF coordinates back to WGS84 actually produces an ellipsoid height parameter:

53.7184377538 -0.876990985207 -2971.299770

Notice the height is negative, and is actually below the WGS84 surface. This makes sense if you think about points on a ball. The "average" of these points will probably be somewhere under the surface. This is why it looks like the latitude value is off.