recompute_points.py

# these projections come from http://mathworld.wolfram.com/GnomonicProjection.html
def gnomonic(lat, lon, clat, clon):  # point, centroid both given in radians
    # print "g",lat,lon,clat,clon
    cosc = sin(clat) * sin(lat) + cos(clat) * cos(lat) * \
        cos(lon - clon)  # cosine of angular distance
    x = (cos(lat) * sin(lon - clon)) / cosc
    y = (cos(clat) * sin(lat) - sin(clat) * cos(lat) * cos(lon - clon)) / cosc
    return (x, y)

def ignomonic(gx, gy, clat, clon):
    # print "ig",gx,gy,clat,clon
    p = sqrt(gx ** 2 + gy ** 2)
    if p == 0:
        return (clat, clon)
    c = atan2(p, 1)
    lat = asin(cos(c) * sin(clat) + (gy * sin(c) * cos(clat)) / p)
    # lon = clon + atan( (gx * sin(c)) / (p * cos(clat) * cos(c) - gy*sin(clat)*sin(c))   )
    lon = clon + \
        atan2(
            (gx * sin(c)), (p * cos(clat) * cos(c) - gy * sin(clat) * sin(c)))
    return (lat, lon)

# from https://en.wikipedia.org/wiki/Haversine_formula
def haversine(p1, p2):  # assume we are already in radians
    (lat1, lon1) = p1
    (lat2, lon2) = p2
    dLat = (lat2 - lat1)
    dLon = (lon2 - lon1)
    # lat1 = radians(lat1)
    # lat2 = radians(lat2)
    h = sin(dLat / 2) * sin(dLat / 2) + \
        sin(dLon / 2) * sin(dLon / 2) * cos(lat1) * cos(lat2)
    # c = 2 * atan2(h ** .5, (1-h) ** .5)  # 8.19155
    c = 2 * asin(h ** .5)    # 6.9557
    return c

def recomputePoints():
    global pts  # array of points on a sphere, specified as (lat,lon) pairs
    new_pts = []
    for i in xrange(nbrPoints):
        # center point for gnomonic projection (projects to 0,0)
        (clat, clon) = pts[i]
        fx, fy = (0, 0)
        for j in xrange(nbrPoints):
            if j == i:
                continue
            jlat, jlon = pts[j]
            hd = haversine(pts[i], pts[j])
            if hd < maxInfluence:
                # closer points have more influence..
                force = map(pow(hd / maxInfluence, 2), 0, 1, 1, 0) * iScale
                gx, gy = gnomonic(jlat, jlon, clat, clon)
                d = sqrt(gx ** 2 + gy ** 2)  # dist(0,0,gx,gy)
                nfx, nfy = -(gx / d) * force, -(gy / d) * force
                fx += nfx
                fy += nfy
        new_pts.append(ignomonic(fx, fy, clat, clon))

    pts = new_pts