Transform one Geocentric to multiple earth-locations

[zaltanar]

Hi,

First thanks for your work, your library is really amazing.

But for a specific application, i try to generate astrometric table and graph in alt/az for earth observers. I know the dates in advance, but not the locations.
To gain some time and I/O reads from the Kernel file i try to find a way to pre-compute Geocentric position and use them on demand to transform from Geocentric to a Apparent with the observer coordinate.

I try adding/substrating different positions object with no success. I think i have to Positions class directly, but it's not clear in my head, and the documentation is not talking about this approach.

How can we use an existing Geocentric Position (Apparent or not) and transform it into multiple Earth locations ? I guess it's just a Vector transformation, but don't know how to handle it.

Sorry i don't have a lot of code to show.

from skyfield.api import load, wgs84

def astrometric_to_apparent_with_topo(astrometric, t, location):
    #don't know how to take this function ...
    pass
    return

if __name__ == "__main__":
    planets = load('data/de440.bsp')
    timescale = load.timescale()
    t0 = timescale.now()

    astrometric_mars = planets["EARTH"].at(t0).observe(planets["MARS BARYCENTER"])

    geographic_toulouse = wgs84.latlon(43.6167, 5.7667)

    print(astrometric_to_apparent_with_topo(astrometric_mars, t0, `geographic_toulouse))

[zaltanar]

I manage to make a working (no error) code by substracting two apparent (earth -> body) - (earth + wgs). But i have some differences between results.
I don't know from where the differences come. 30mas maximum (for the moon), i may live with it, but i don't know from where these differences come.

from skyfield.api import load, wgs84

planets = load('data/de440.bsp')
timescale = load.timescale()
t0 = timescale.now()
sun = planets["SUN"]
mercury = planets["MERCURY"]
venus = planets["VENUS"]
earth = planets["EARTH"]
moon = planets["MOON"]
mars = planets["MARS BARYCENTER"]
jupiter = planets["JUPITER BARYCENTER"]
saturn = planets["SATURN BARYCENTER"]
uranus = planets["URANUS BARYCENTER"]
neptune = planets["NEPTUNE BARYCENTER"]
toulouse = earth + wgs84.latlon(43.6167, 5.7667)

bodies = [sun, mercury, venus, moon, mars, jupiter, saturn, uranus, neptune]


# Precompute for optimisation
earth_astrometric = [earth.at(t0).observe(o) for o in bodies]
earth_apparent = [o.apparent() for o in earth_astrometric]

# Compute Vector Earth -> Toulouse and Vector substraction (Geocentric - Eartg/Tls)
tls_earth_astrometric = earth.at(t0).observe(toulouse)
tls_earth_apparent = tls_earth_astrometric.apparent()
tls_astrometric = [o - tls_earth_astrometric for o in earth_astrometric]
tls_apparent = [o - tls_earth_apparent for o in earth_apparent]

#Compute all
tls_all_astrometric = [toulouse.at(t0).observe(o) for o in bodies]
tls_all_apparent = [o.apparent() for o in tls_all_astrometric]

for i in range(0, len(bodies) -1 ):
    print(bodies[i])
    print(tls_astrometric[i].position)
    print(tls_all_astrometric[i].position)

    print(tls_apparent[i].position)
    print(tls_all_apparent[i].position)

[brandon-rhodes]

Thank you for providing a working example! One source of difference could be light travel time. An observer at the geocenter sees a slightly earlier Moon, whose light left earlier, than does an observer up on the planet's surface. But the difference is maybe not big enough to explain 30mas? If the Moon moves through 360° in 29 days, and if the Earth's radius is 0.02 light seconds, then the error from light travel time would be:

360 * 60 * 60 * 1000 / 29 / 24 / 60 / 60 * 0.02 = 10.3 mas

It is also possible that Earth deflection is the difference — but that would appear only in apparent positions. Are you seeing a 30mas difference for both astrometric and apparent positions, or only apparent ones?

[zaltanar]

For the moment, Astrometric seems more precise, but i can't exclude a 10 mas error.

At first i thought about atmospheric, but then i remember that it's only for alt/az.

Thanks, i will run more precise tests tonight to measure this.

[mkbrewer]

Diurnal aberration can add up to 320 mas at the equator.

[brandon-rhodes]

@mkbrewer — Wow, is the effect that large? Thanks, that's probably the difference, then. I’ve never run the numbers to compute it; with the speed of a point on the equator amounting to only around 1/100th of the speed of the Earth in its orbit, I had not realized that in absolute terms it amounted to several hundred mas.

[mkbrewer]

Equation (7.50) on page 267 of the ES.

[zaltanar]

I make some more tests based on your responses.

I use the moon as body, a time of the last lunar node (2021/01/24 21h47 GMT) and plot several position on earth.
I plot mas errors for apparent and astrometric in a notebook : https://gist.github.com/zaltanar/5de40b7df43675f3fe881714015f3422
For apparent it goes up to 600mas errors and 400mas erros for astrometric.

But, a problem i didn't realize when i post my first resolution is that when you substract ICRF childs class you get back an ICRF object, that doesn't fit my use. So i guess i have to do myself the substraction of two ICRF childs. And i can confirm with sub function in ICRF

python-skyfield/skyfield/positionlib.py

Lines 169 to 177 in 2c05c63

	def __sub__(self, body):
	"""Subtract two ICRF vectors to produce a third."""
	# TODO: set center and target of result
	p = self.position.au - body.position.au
	if self.velocity is None or body.velocity is None:
	v = None
	else:
	v = self.velocity.au_per_d - body.velocity.au_per_d
	return ICRF(p, v, self.t)

[mkbrewer]

Yes. In order to do this properly, you need to take the difference in the ICRS vectors including the difference in light time and then transform through GCRS to add the proper aberration and deflection of light due to the gravitational field of the Earth to apparent. Therefore, the initial geocentric calculations are redundant and a waste of time.

[mkbrewer]

Of course, if you don't care about the difference in light time, then the initial geocentric calculations will save ephemeris lookup time. But if you are using one of the JPL ephemerides, those lookups are quite fast.. One possibility for rigorous speed up would be to do the lookup of the Earth once and then apply different observer locations to calculate the different light time delays to the target body. Once it has the proper ephemeris record in memory, this is just an additional Chebyshev polynomial evaluation

[mkbrewer]

You could also just do the target body lookup once and then use the target body velocity to account for the different light time delays. So there are several possibilities here, but they are all at a fairly low level.

[brandon-rhodes]

But, a problem i didn't realize when i post my first resolution is that when you substract ICRF childs class you get back an ICRF object, that doesn't fit my use.

Thanks for pointing that out! It should probably call build_position() instead to allow the class to be chosen based on the center. And, the center: it needs to preserve the center and target! It looks like it was written in 2013 to satisfy a quick short-term need (the commit message is “Add rudimentary subtraction ability to ICRS”, so I apparently recognized it was rudimentary) so it’s now obviously time for an update.

[brandon-rhodes]

As shown just above, I've committed a fix for vector subtraction that will use more specific vector classes like Geocentric if subtraction indeed results in a vector whose origin is the geocenter. To try the new feature out:

pip install -U https://github.com/skyfielders/python-skyfield/archive/master.zip

As for your original question, is this an accurate summary of the discussion so far?

1. Because each location on Earth (besides the poles) has its own velocity around the Earth's axis in addition to the velocity of the geocenter itself, the aberration of light will be different; and unless it happens to be at the exact same distance from the target as the geocenter, the light-travel-time will be different and therefore so will the position of the target; and so vector subtraction alone won't produce a topocentric vector from a geocentric vector.
2. As mentioned above, there is no I/O penalty for computing planet positions once the data has been loaded into RAM; they are one of the fastest things Skyfield does.
3. Which might mean that your application is being slowed down by something else. Have you tried profiling the code to see which Skyfield routine is taking the most time?
4. I assume that a 600 mas error is unacceptable for your requirements, but it is not clear what your requirements are — you might want to review the precision to which coordinates will be given in your output table, so that we can know what the acceptable error is, which would suggest which effects can be ignored while still giving an accurate enough answer.

[zaltanar]

You are too reactive for me.

I had my response since the third post, the discussion after it's just bonuses and some curiosity.

1. Yes, in my initial thought i didn't take into account for diurnal aberration
2. I didn't know at first the caching mechanism
3. No, i'm still in api testing, i plan to tests different methods.
4. After some calculations, as the moon travels 550 mas per second in the sky, i can play with this accuracy error.

I get everything i need now, thanks a lot, i learn everything :

1. Skyfield can substracts Positions
2. The api load data in ram
3. Diurnal aberration exists

I close the issue, i as i have my response. I will share tests results later here. But as it's a side project it may take some days.
Thanks again.