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
uhoh
  • 3,713
  • 6
  • 42
  • 95
  • What is the actual question? The variation with latitude from about 6357 km to 6378 km looks reasonable (see, for example, https://en.wikipedia.org/wiki/Figure_of_the_Earth), so which numbers in your output are you questioning? – Warren Weckesser Jan 19 '16 at 08:03
  • 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. – uhoh Jan 19 '16 at 08:12
  • 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. – uhoh Jan 19 '16 at 08:26

1 Answers1

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.

Brandon Rhodes
  • 83,755
  • 16
  • 106
  • 147
  • 1
    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! – uhoh Jan 23 '16 at 02:35
  • 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](http://space.stackexchange.com/q/13523/12102) that would be great. Thanks again! – uhoh Jan 23 '16 at 02:44