Используя Visual Python, я думаю, что теперь я решил это:
# Rotation first described for geo purposes: http://www.uwgb.edu/dutchs/mathalgo/sphere0.htm
# /4201267/vraschenie-3d-vektora
# http://vpython.org/
from visual import *
from math import *
import sys
def ll2cart(lon,lat):
# http://rbrundritt.wordpress.com/2008/10/14/conversion-between-spherical-and-cartesian-coordinates-systems/
x = cos(lat) * cos(lon)
y = cos(lat) * sin(lon)
z = sin(lat)
return x,y,z
def cart2ll(x,y,z):
# http://rbrundritt.wordpress.com/2008/10/14/conversion-between-spherical-and-cartesian-coordinates-systems/
r = sqrt((x**2) + (y**2) + (z**2))
lat = asin(z/r)
lon = atan2(y, x)
return lon, lat
def distance(lon1, lat1, lon2, lat2):
# http://code.activestate.com/recipes/576779-calculating-distance-between-two-geographic-points/
# http://en.wikipedia.org/wiki/Haversine_formula
dlat = lat2 - lat1
dlon = lon2 - lon1
q = sin(dlat/2)**2 + (cos(lat1) * cos(lat2) * (sin(dlon/2)**2))
return 2 * atan2(sqrt(q), sqrt(1-q))
if len(sys.argv) == 1:
sys.exit()
else:
csv = sys.argv[1]
# Points A and B defining the rotation:
LonA = radians(float(sys.argv[2]))
LatA = radians(float(sys.argv[3]))
LonB = radians(float(sys.argv[4]))
LatB = radians(float(sys.argv[5]))
# A and B are both vectors
# The crossproduct AxB is the rotation pole vector P:
Ax, Ay, Az = ll2cart(LonA, LatA)
Bx, By, Bz = ll2cart(LonB, LatB)
A = vector(Ax,Ay,Az)
B = vector(Bx,By,Bz)
P = cross(A,B)
Px,Py,Pz = P
LonP, LatP = cart2ll(Px,Py,Pz)
# The Rotation Angle in radians:
# http://code.activestate.com/recipes/576779-calculating-distance-between-two-geographic-points/
# http://en.wikipedia.org/wiki/Haversine_formula
RotAngle = distance(LonA,LatA,LonB,LatB)
f = open(csv,"r")
o = open(csv[:-4] + "_translated.csv","w")
o.write(f.readline())
for line in f:
(lon, lat) = line.strip().split(",")
# Point C which will be translated:
LonC = radians(float(lon))
LatC = radians(float(lat))
# Point C in Cartesian coordinates:
Cx,Cy,Cz = ll2cart(LonC,LatC)
C = vector(Cx,Cy,Cz)
# C rotated to D:
D = rotate(C,RotAngle,P)
Dx,Dy,Dz = D
LonD,LatD = cart2ll(Dx,Dy,Dz)
o.write(str(degrees(LonD)) + "," + str(degrees(LatD)) + "\n")