itrs_xyz method of GeographicPosition WGS84 object returns bad values

[vholod]

Hi, I am new in Skyfield. I have a lot of fun using this great package. I have an issue that I don't understand and I have stuck because of that issue.

I compute a subpoint ( on the Earth's surface) directly below a satellite and I try to calculate a distance from that point to a reference point. My goal is to find a minimal distance. The reference point is also on the Earth's surface. The problem is that I get wrong distance. I will demonstrate this issue with a simple example where the reference point and the subpoint are pretty the same.

Here is the example: test_lat_lon is <GeographicPosition WGS84 latitude -10.4880 N longitude -27.7536 E elevation 0.0 m> subsat_on_earth is <GeographicPosition WGS84 latitude -10.4880 N longitude -27.7536 E elevation -0.0 m>

Note that test_lat_lon is my reference point. The subsat_on_earth and the test_lat_lon points seem to be the same when comparing the lat lon and elevation. But the itrs_xyz method of both gives different results.

test_lat_lon.itrs_xyz.km is array([ 5550.70173713, -2920.80937139, -1153.36430579]) subsat_on_earth.itrs_xyz.km is array([ 5917.50303437, -3113.82220461, -1230.09455252])

How to solve this problem and why do I get such differences? Thanks in advance.

[brandon-rhodes]

So that I can investigate, can you write a small Python script that builds these two GeographicPosition objects? Without knowing which Skyfield code produced each one, I won't be able to track down where those two code paths might be in error. Thanks!

[vholod]

Thank you very much for your response, Here is a python script with this issue (small script).

import numpy as np
from datetime import datetime, timedelta

import skyfield
from skyfield.api import load, wgs84, EarthSatellite, Topos
from skyfield.framelib import itrs

ts = load.timescale(builtin=True)

# use spesific TLE file.       
TLE = """ISS (ZARYA)
1 25544U 98067A   21340.39932169  .00004577  00000+0  91351-4 0  9997
2 25544  51.6435 211.2046 0003988 277.6071 254.8235 15.48946544315288"""

name, L1, L2 = TLE.splitlines()
L1 = L1.strip()
L2 = L2.strip() 

satellite = EarthSatellite(L1, L2, name,ts)

# Interduce spesific point in the Earth to lookat with all satellites:
test_lat_lon = [-10.48799047, -27.75358867]
test_lat_lon = wgs84.latlon(*test_lat_lon)

# Define times range to seek this spesific occurrence.
minutes = np.arange(0, 60*10, 5) # about two orbits
fdate = satellite.epoch.utc_strftime('%Y,%m,%d,%H,%M,%S')
fdate = fdate.split(',')
times = ts.utc(int(fdate[0]),\
                int(fdate[1]), int(fdate[2]), int(fdate[3]), int(fdate[4]) + minutes)

# test at time 0:
another_test_lat_lon = [-10.487990468777502, -27.753588673056026]
another_test_lat_lon = wgs84.latlon(*another_test_lat_lon)
position1 = test_lat_lon.at(times[0])
position2 = another_test_lat_lon.at(times[0])    
print('start test')
print((position1 - position2).distance().km)
print(np.linalg.norm(test_lat_lon.itrs_xyz.km - another_test_lat_lon.itrs_xyz.km))
print('end test')   
# At time 0 the distance should hbe very small and it is small

st_lons  = np.zeros([len(times),1])
st_lats  = np.zeros([len(times),1])
st_dists = np.zeros([len(times),1])

for t_indx,t in enumerate(times):
    geocentric = satellite.at(t) # Geocentric ICRS position at t,
    # Now Compute the subpoint directly below the satellite  
    # the point on the Earth with the same latitude and longitude
    # as the satellite,
    # but with a height above the WGS84 ellipsoid of zero:
    subsat_on_earth = wgs84.subpoint_of(geocentric) # GeographicPosition WGS84     

    #position1 = subsat_on_earth.at(t).frame_xyz(itrs).km
    #position2 = test_lat_lon.at(t).frame_xyz(itrs).km

    position1 = subsat_on_earth.itrs_xyz.km
    position2 = test_lat_lon.itrs_xyz.km

    st_dists[t_indx]  = np.linalg.norm(position1 - position2)
    st_lons[t_indx]   = subsat_on_earth.longitude.degrees
    st_lats[t_indx]   = subsat_on_earth.latitude.degrees

k = np.argmin(st_dists)
print('closest sat lat {}, lon {}'.format(st_lats[k], st_lons[k]))
print('distance {} km'.format(st_dists[k]))
print('The problem')    
print('The minimal distance should be close to zero, but it is {} km'.format(st_dists[k]))
print('done')

Note that the lat_lon = [-10.487990468777502, -27.753588673056026] is taken from the subpoint of the satellite at iteration index 9. So the minimal distance should be very small = close to zero

[brandon-rhodes]

I've just pushed what could be a fix for the problem you're seeing! Could you try installing the development version of Skyfield and testing your script again? Thanks!

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

[vholod]

Thank you very much, It is working now without the issue a had.

[brandon-rhodes]

Excellent! I'll release a new version right now to get this fix out to the world.