6

I'm trying to reduce a set of lat/lons (each with it's own weight) into one lat/lon that is the center of mass of all of the points.

however when I map the results the average point falls no where near where I think it should be (I've tested on multiple sets of points).

here's my code, in python, based on these two answers:
How to find a point half way between two other points?
Computing an averaged latitude and longitude coordinates

import math
import nummpy
import math

def toCartesian(t):
        latD,longD = t
        latR = math.radians(latD)
        longR = math.radians(longD)
        return (
            math.cos(latR)*math.cos(longR), 
            math.cos(latR)*math.sin(longR), 
            math.sin(latR)
            )

def toSpherical(t):
    x,y,z = t

    r = math.hypot(x,y)
    if r == 0:
        if z > 0:
            return (90,0)
        elif z< 0:
            return (-90,0)
        else:
            return None
    else:
        return (math.degrees(math.atan2(z, r)), math.degrees(math.atan2(y,x)))


xyz = numpy.asarray([0.0,0.0,0.0])
total = 0
for p in points:
        weight = p["weight"]
        total += weight
        xyz += numpy.asarray(toCartesian((p["lat"],p["long"])))*weight

avgXYZ = xyz/total
avgLat, avgLong = toSpherical(avgXYZ)
print avgLat,avgLong
Community
  • 1
Yanofsky
  • 223
  • 2
  • 7
  • 2
    Do you mean weighted mean center? – gm70560 Jun 14 '13 at 14:35
  • 2
    Just glancing over the code for obvious problems, it strikes me that atan2(r,z) ought to be atan2(z,r). With numerical code like this, it's important to write little testing modules for your functions (like toSpherical and toCartesian) and exercise them well so that you know you can trust them. Why don't you do that and see whether your problems go away? – whuber Jun 14 '13 at 14:35
  • @gm70560 The ArcGIS approach you refer to requires (a) data to be projected and (b) that they not extend very widely over the earth's surface. The approach shown in the code here is intended for global applications where the ArcGIS "Mean Center" will fail or give incorrect results. – whuber Jun 14 '13 at 14:37
  • @whuber I am confused as to the use of "center of gravity" in the question, that is all. – gm70560 Jun 14 '13 at 14:51
  • @whuber Yes, that was a problem, however the issue persists after fixing to confirm toSpherical(toCartesian((40.7142,-74.0064))) == (40.7142,-74.0064) – Yanofsky Jun 14 '13 at 15:12
  • 1
    I ran a whole bunch of tests on this and found that my data parsing up until this point was the problem, not the calculation of the weighted average...the above is correct – Yanofsky Jun 14 '13 at 17:09
  • Is there still a question here, then? It looks like not, so we ought to close this thread. However, it's nice that you posted the code. If you want to keep it up, may I suggest you offer it as an answer to one of the related questions you linked to? That will make it available to other interested readers in the future. – whuber Jun 14 '13 at 17:28

0 Answers0