Пока вы можете предположить, что ваши координаты широты и долготы образуют прямоугольную сетку (что должно быть разумным приближением в масштабе города, хотя оно не будет работать для больших площадей или очень точных данных), вы можете предположить, что ваше преобразование между координатами будет проективное преобразование. См. этот пост (или этот, если вы предпочитаете оставаться на SO) для получения подробной информации о том, как вычислить соответствующую матрицу преобразования и использовать ее для преобразования координат.
Основываясь на этом сообщении, я провел небольшие вычисления с помощью мудрец:
def pm1(a, b, c, d):
M = Matrix([a , b, c]).transpose()
f = M.adjoint() * d
return M*diagonal_matrix(f)
def pm(a1, a2, b1, b2, c1, c2, d1, d2):
return pm1(a2, b2, c2, d2)*(pm1(a1, b1, c1, d1).adjoint())
P1 = vector(QQ, [0, 0, 1])
P2 = vector(QQ, [0, 1, 1])
P3 = vector(QQ, [1, 1, 1])
P4 = vector(QQ, [1, 0, 1])
Q1 = vector(QQ, [40698291, -74079051, 1000000]) # Upper left
Q2 = vector(QQ, [40659855, -73979638, 1000000]) # Lower left
Q3 = vector(QQ, [40855232, -73835582, 1000000]) # Lower right
Q4 = vector(QQ, [40882919, -73940263, 1000000]) # Upper right
M = pm(Q1, P1, Q2, P2, Q3, P3, Q4, P4)
M.change_ring(RDF)/1e40
Полученная формула выглядит следующим образом:
z = 559.910562534*lat - 539.510073656*lon - 64123.7576703
x = (-5629.59680416*lat - 2176.56828347*lon + 67876.8560722)/z
y = ( 8466.61769096*lat - 11263.0392472*lon - 1178932.12938)/z
Полученные координаты x и y измеряются от нуля до единицы, где ноль обозначает левую соотв. верхний край и один правый соотв. нижний край.
Если предположения о прямоугольной сетке недостаточно, вам необходимо знать подробности о проекции, используемой для создания карты, чтобы вы могли связать положения карты со сферическими координатами.
17.06.2013
