Поверните сферу с координаты1 на координату2, где будет координата3?

У меня есть три координаты (широта, долгота) на сфере. Если бы вы повернули всю сферу из координаты 1 в координату 2, где теперь будет находиться координата 3?

Я пробовал это в Python, используя Great Circle (http://www.koders.com/python/fid0A930D7924AE856342437CA1F5A9A3EC0CAEACE2.aspx?s=coastline), но я получаю странные результаты, поскольку недавно вычисленные точки группируются вместе на экваторе. Я предполагаю, что это как-то связано с расчетом азимута?

Может быть, кто-нибудь знает, как это правильно рассчитать?

Заранее спасибо!

ИЗМЕНИТЬ

Я нашел следующее: http://www.uwgb.edu/dutchs/mathalgo/sphere0.htm

Я предполагаю, что теперь мне нужно вычислить ось вращения и угол поворота из двух точек в декартовых координатах (и 0,0,0)? Я предполагаю, что это должно быть очень просто, что-то связанное с определением плоскости и определением нормальной линии? Может быть кто-то знает, где я могу найти необходимые уравнения?

ИЗМЕНИТЬ 2

Coord1 и coord2 образуют большой круг. Есть ли простой способ найти положение оси нормали большого круга на сфере?

ИЗМЕНИТЬ 3

Похоже, мне удалось ее решить ;) http://articles.adsabs.harvard.edu//full/1953Metic...1...39L/0000039.000.html сделал свое дело.


person joosthoek    schedule 28.03.2012    source источник


Ответы (2)


Хорошо, я не знаю точной формулы, я думаю, что это будет простое умножение матриц, но вот как вы можете понять без этого.

  1. Преобразуйте координаты так, чтобы полюса вращения находились на 90,0 и -90,0 соответственно, и чтобы линия вдоль вашего вращения от coord1 до coord2 находилась на «экваторе» (это должно быть просто дельта-широта, дельта-длина)

  2. тогда вращение - это просто изменение долготы, и вы можете применить одну и ту же дельту к любой координате3, а затем просто вернуться к исходным координатам (через отрицательную дельта-широту и отрицательную дельту-длину)

1 и 2 - это то, что будет делать ваша матрица - если вы можете вычислить матрицы для каждого шага, вы можете просто умножить их и получить окончательную матрицу.

person scibuff    schedule 28.03.2012
comment
Звучит логично, спасибо! Я сейчас пытаюсь сделать это преобразование, но еще не понял. Кто-нибудь знает код, который это делает? - person joosthoek; 28.03.2012

Используя Visual Python, я думаю, что теперь решил это:

# Rotation first described for geo purposes: http://www.uwgb.edu/dutchs/mathalgo/sphere0.htm
# http://stackoverflow.com/questions/6802577/python-rotation-of-3d-vector
# 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")
person joosthoek    schedule 25.01.2013