Почему переход от земных координат к галактоцентрическим координатам в астропсии не сохраняет расстояния?

Я конвертирую GCRS объектов в галактоцентрические координаты и обнаруживаю, что расстояния между двумя точками не сохраняются при этом преобразовании.

import astropy.units
import astropy.coordinates
import astropy.time
from numpy.linalg import norm

t = astropy.time.Time('1999-01-01T00:00:00.123456789')

def earth2galaxy(lat):
    '''                                                                         
    Convert an Earth coordinate to a galactocentric coordinate.                 
    '''
    # get GCRS coordinates                                                      
    earth = astropy.coordinates.EarthLocation(lat=lat*astropy.units.deg,
                                            lon=0,
                                            height=0)
    pos, _ = earth.get_gcrs_posvel(obstime=t)
    cartrep = astropy.coordinates.CartesianRepresentation(pos.xyz,
                                                          unit=astropy.units.m)

    gcrs = astropy.coordinates.GCRS(cartrep)

    # convert GCRS to galactocentric                                            
    gc = gcrs.transform_to(astropy.coordinates.Galactocentric)

    return earth, gcrs, gc

earthA, gcrsA, gcA = earth2galaxy(0)
earthB, gcrsB, gcB = earth2galaxy(0.01)

print(norm(earthA-earthB))
print(norm(gcrsA.cartesian.xyz-gcrsB.cartesian.xyz))
print(norm(gcA.cartesian.xyz-gcB.cartesian.xyz))

Этот код дает

1105.74275693
1105.74275232
971.796949054

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

Раньше я обходил это, используя точки A и B, преобразуя точки A и C = A + c*AB, где c — большое число. Затем я восстановил преобразованный B', отменив масштабирование B' = A' + A'C' / c. Однако мне кажется, что я должен обратиться к фактическому корню проблемы, а не к этому обходному пути.


python debugging astropy
person wht    schedule 14.06.2019    source источник

Ответы (2)


arrow_upward
3
arrow_downward

Это может быть просто проблема точности с плавающей запятой. Если я смотрю на декартовы значения, то x, y и z имеют порядок 1e6, 1e6 и 1e2 для системы координат GCRS, но они имеют порядок 1e20, 1e10 и 1e17 соответственно для галактической системы координат.

Учитывая точность 1e-15 для 8-байтовых чисел с плавающей запятой (numpy.finfo('f8').eps), это означает, что x-значение галактической координаты может быть точным только до 1e5 (метров). Тогда взятие нормы (с преобладающей неопределенностью x значения) также приведет к точности порядка 1e5 метров, что намного больше, чем фактическое расстояние.

Тот факт, что рассчитанные значения все еще близки друг к другу, в значительной степени является удачей (хотя у этого будет основная причина, такая как некоторое усреднение отклонений).

Это также согласуется с тем фактом, что вы не видите проблемы (или меньше проблемы) для больших смещений. Хотя сам тестировал, но все равно вижу разницу порядка 1e4~1e5). Чтобы быть точным, используя широту 0 и 10, я получаю:

GCRS:     1104451.74511518
Galactic: 1108541.8206286128

Если мое предположение верно, то мой совет прост: используйте подходящую систему координат для ваших координат и принимайте во внимание соответствующие неопределенности (как точность машины, так и точность используемой системы координат).

person 9769953    schedule 18.06.2019

arrow_upward
1
arrow_downward

Я думаю, что "0 0" в основном попадает в самую точку, но для этого случая есть альтернативное решение: использовать dtype с плавающей запятой, который имеет большую точность - например. float128 dtype (который вы бы передали ключевому слову dtype слова CartesianRepresentation), вероятно, расстояние было бы намного меньше.

К сожалению, есть ошибка, из-за которой типы dtypes возвращаются к float64 (о которой я сообщил: https://github.com/astropy/astropy/issues/8870), но если это исправлено, в зависимости от вашего компьютера и поддерживаемой им точности, это может помочь использование более точного dtype.

person eteq    schedule 19.06.2019
comment
Немного не по теме, но будет это быть актуальная (NumPy) документация по float128? Потому что это также, возможно, немного разочаровывает: несмотря на имена, np.float96 и np.float128 обеспечивают такую ​​же точность, как np.longdouble, то есть 80 бит на большинстве машин x86 и 64 бита в стандартных сборках Windows. - person 9769953; 20.06.2019
comment
Да, это правда. Если у вас есть машина Solaris, очевидно, это правда 128 бит, но на x86_64 (я подозреваю, что это то, что есть у большинства людей), это 80, что составляет всего 1e-15 -> 1e-18. Но для исходного случая OP это имеет большое значение, поскольку я думаю, что они находятся прямо на краю 1e-15? - person eteq; 21.06.2019