Shape of earth seems wrong in Skyfield - is my python correct?
Asked Answered
C

1

6

Mapping the distance from the center of the earth to various (lat, lon) positions using Skyfield shows variation with latitude but independent of longitude (sub-millimeter). This may be a documented approximation in the package, a bug in my script, or something else altogether. Am I doing something wrong here? (besides, of course, using jet)

uhoh topo

import numpy as np
import matplotlib.pyplot as plt
from skyfield.api import load, now

data  = load('de421.bsp')
earth = data['earth']
jd    = now()

epos = earth.at(jd).position.km

lats = np.linspace( -90,  90, 19)
lons = np.linspace(-180, 180, 37)
LATS, LONS = np.meshgrid(lats, lons)

s = LATS.shape

points = zip(LATS.flatten(), LONS.flatten())

rr = []
for point in points:
    la, lo = point
    pos = earth.topos(la, lo).at(jd).position.km
    r   = np.sqrt( ((pos-epos)**2).sum() )
    rr.append(r)

surf = np.array(rr).reshape(s)

extent = [lons.min(), lons.max(), lats.min(), lats.max()]

plt.figure()
plt.imshow(surf.T, origin='lower', extent=extent)
plt.colorbar()
plt.title('uhoh topo')
plt.savefig('uhoh topo')
plt.show()

As a cross-check, I tried some random pairs of locations with the same latitude:

pe = earth.at(jd).position.km
for i in range(10):

    lon1, lon2 = 360.*np.random.random(2)-180
    lat = float(180.*np.random.random(1)-90.)

    p1  = earth.topos(lat, lon1).at(jd).position.km
    p2  = earth.topos(lat, lon2).at(jd).position.km

    r1   = np.sqrt( ((p1-pe)**2).sum() )
    r2   = np.sqrt( ((p2-pe)**2).sum() )

    print lat, lon1, lon2, r2-r1

and got this (the fourth column shows differences of microns):

45.8481950437 55.9538249618 115.148786114 1.59288902069e-08
-72.0821405192 4.81264755835 172.783338907 2.17096385313e-09
51.6126938075 -54.5670258363 -134.888403816 2.42653186433e-09
2.92691713179 -178.553103457 134.648099589 1.5916157281e-10
-78.7376163827 -55.0684703115 125.714124504 -6.13908923697e-10
48.5852207923 -169.061708765 35.5374862329 7.60337570682e-10
42.3767785876 130.850223447 -111.520896867 -1.62599462783e-08
11.2951212126 -60.0296460731 32.8775784623 6.91579771228e-09
18.9588262131 71.3414406837 127.516370219 -4.84760676045e-09
-31.5768658495 173.741960359 90.3715297869 -6.78483047523e-10
Clyte answered 19/1, 2016 at 7:12 Comment(3)
What is the actual question? The variation with latitude from about 6357 km to 6378 km looks reasonable (see, for example, en.wikipedia.org/wiki/Figure_of_the_Earth), so which numbers in your output are you questioning?Quesenberry
The distance should vary with longitude as well - certainly more than a few tens of microns. The oblate spheroid is just a first order approximation to the real shape. And that doesn't even take in to consideration elevation or the shape of the geoid. I'm trying to understand if this is an implicit approximation in the package or a bug or I'm not using the package correctly.Clyte
The fact that the variation with longitude is small but not zero may just be coming from round-off error (these are astronomical distances), so it's possible it is using the ellipsoid approximation. But I don't know if that is the case or not. People do GPS measurements and satellite positions to millimeter resolution (and sometimes accuracy as well!) - this would be a pretty big approximation if it were the case.Clyte
K
1

The topos() method includes a elevation_m=0.0 that you are not specifying. So in each of your calls, it takes its default value of 0.0 meters above sea level. And the definition of “sea level” does not vary with longitude — at a given latitude, sea level is at a fixed distance from the geocenter the entire way around the world.

I am not sure why you call a micron-level error a “pretty big approximation”, but you are indeed running into the limited precision of your machine’s 64-bit floating point math when you subtract two distances that are each about 1 au in length — you are seeing a rounding error down in the last bit or two.

Khalid answered 22/1, 2016 at 14:55 Comment(2)
I checked earth.topos.__doc__ - should have typed help(earth.topos) which has elevation option. On-line docs show lat and lon args, but not elevation m. I know docs are in progress - maybe just add the elevation option to satellite example and use a city that isn't (practically) at sea level (Denver instead of Boston?). Also I thought that even the elevation = 0 sea level surface would be a geoid - a gravity equipotential surface, and would have 2D structure. Now I see WGS84 bi-axial ellipsoid is "sea level" for Geodetic purposes. Thanks - Skyfield rocks!Clyte
The "pretty big approximation" refers to using only the ellipsoid for calculation of near-earth and astrometric events, and is in response to @WarrenWeckesser 's comment, before I knew that an elevation option was available. Surely you'd agree that not including elevation would be a big approximation! btw if you can take a look at this question that would be great. Thanks again!Clyte

© 2022 - 2024 — McMap. All rights reserved.