4

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
  • 2
    Could you update your question with some code samples or sample values? I'd like to see how you're going from WGS84 to ECEF and weighting the XYZ average. – Mintx Aug 28 '15 at 14:14
  • @Mintx I've added the algorithms I'm using to convert from WGS84 to ECEF and reverse. These work accurately when I test them by converting a lat/long pair to and from (i.e. returns the same lat long pair after converting back). It seems to be the averaging of the Cartesian coordinates that causes the problem. – J.P. Aug 28 '15 at 14:38
  • 2
    Is the XYZ weighted average falling below the surface of the ellipsoid? That could cause the latitude to be lower than you expect. – mkennedy Aug 28 '15 at 16:00
  • @mkennedy I'm sorry, I'm not entirely sure I understand what you mean by that, could you elaborate? – J.P. Aug 28 '15 at 16:10
  • @mkennedy I understand you now, Mintx has explained further in their answer, thank you – J.P. Aug 28 '15 at 16:39
2

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.

  • Is there any way of compensating for this, or should I find another way to find the average? – J.P. Aug 28 '15 at 17:48
  • 1
    If you have data on either side of the date line, you could add 360 to the negative longitude values. You can't do this for all data, as this would shift the same problem to data on either side of the prime meridian. You could also convert the lat/longs into a projected XY coordinate system. (UTM for example) – Mintx Aug 28 '15 at 17:59
  • I'll mark this as the answer – J.P. Aug 30 '15 at 12:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.