Skip to content

Instantly share code, notes, and snippets.

@jbum
Created September 28, 2017 17:21
Show Gist options
  • Save jbum/f6a04a11ccb8e598c3257f5a5ae740ba to your computer and use it in GitHub Desktop.
Save jbum/f6a04a11ccb8e598c3257f5a5ae740ba to your computer and use it in GitHub Desktop.
# 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment