#!/usr/bin/env python
# -*- coding: utf-8 -*-
from math import radians, degrees, sin, cos, asin, acos, sqrt, atan2, pi
import numpy as np
[docs]
def gc_distance(lat1, lon1, lat2, lon2):
'Calculate great circle distance and azimuth (exact. parsed ml)'
# distance
lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
a = sin((lat2-lat1)/2)**2 + cos(lat1) * cos(lat2) * sin((lon2-lon1)/2)**2;
if a < 0: a = 0
if a > 1: a = 1
r = 1
rng = r * 2 * atan2(sqrt(a), sqrt(1-a))
rng = degrees(rng)
# azimuth
az = atan2(
cos(lat2) * sin(lon2-lon1),
cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(lon2-lon1)
)
if lat1 <= -pi/2: az = 0
if lat2 >= pi/2: az = 0
if lat2 <= -pi/2: az = pi
if lat1 >= pi/2: az = pi
az = az % (2*pi)
az = degrees(az)
return rng, az
[docs]
def shoot(lon, lat, azimuth, maxdist=None):
"""Shooter Function
Original javascript on http://williams.best.vwh.net/gccalc.htm
Translated to python by Thomas Lecocq
"""
glat1 = lat * np.pi / 180.
glon1 = lon * np.pi / 180.
s = maxdist / 1.852
faz = azimuth * np.pi / 180.
EPS= 0.00000000005
if ((np.abs(np.cos(glat1))<EPS) and not (np.abs(np.sin(faz))<EPS)):
print("Only N-S courses are meaningful, starting at a pole!")
a=6378.13/1.852
f=1/298.257223563
r = 1 - f
tu = r * np.tan(glat1)
sf = np.sin(faz)
cf = np.cos(faz)
if (cf==0):
b=0.
else:
b=2. * np.arctan2 (tu, cf)
cu = 1. / np.sqrt(1 + tu * tu)
su = tu * cu
sa = cu * sf
c2a = 1 - sa * sa
x = 1. + np.sqrt(1. + c2a * (1. / (r * r) - 1.))
x = (x - 2.) / x
c = 1. - x
c = (x * x / 4. + 1.) / c
d = (0.375 * x * x - 1.) * x
tu = s / (r * a * c)
y = tu
c = y + 1
while (np.abs (y - c) > EPS):
sy = np.sin(y)
cy = np.cos(y)
cz = np.cos(b + y)
e = 2. * cz * cz - 1.
c = y
x = e * cy
y = e + e - 1.
y = (((sy * sy * 4. - 3.) * y * cz * d / 6. + x) *
d / 4. - cz) * sy * d + tu
b = cu * cy * cf - su * sy
c = r * np.sqrt(sa * sa + b * b)
d = su * cy + cu * sy * cf
glat2 = (np.arctan2(d, c) + np.pi) % (2*np.pi) - np.pi
c = cu * cy - su * sy * cf
x = np.arctan2(sy * sf, c)
c = ((-3. * c2a + 4.) * f + 4.) * c2a * f / 16.
d = ((e * cy * c + cz) * sy * c + y) * sa
glon2 = ((glon1 + x - (1. - c) * d * f + np.pi) % (2*np.pi)) - np.pi
baz = (np.arctan2(sa, b) + np.pi) % (2 * np.pi)
glon2 *= 180./np.pi
glat2 *= 180./np.pi
baz *= 180./np.pi
return (glon2, glat2, baz)