Finding lat,long location at which a given vector intersects the Earth

[bismurphy]

Hi! Let me preface this by saying that I'm writing as much detail as possible to make sure I don't leave out anything useful - hope this post isn't too incredibly long :)

I'm hoping to write a script that will allow me to detect times when the ISS passes in front of the moon, so I can take a neat picture of its silhouette against the glow of the moon. I chose to solve this graphically: I made a script where I punch in my location and time, and it finds times the ISS is near the moon, then I'm able to use sliders to relocate my latitude and longitude. Whenever I change the sliders, it re-draws the ISS and moon positions in the sky, so by adjusting them I can get the two to line up, and then I just read off where the sliders turned out to land, and then I know where to go.

In order to validate the results of that script, I write a second script. The second script takes the time from the first script. At that time, it finds the positions of the moon and the ISS. Then, it draws a vector from the moon to the ISS, and then extends that vector until it intersects the Earth. Then it finds the latitude and longitude where that intersection occurred. We hope to find that the intersection is at the exact same location as our observer was in the first script.

Here is a proof of concept of my first script. It shows the ISS right in the middle of the moon in the sky. (please note: This is a trimmed down version of the script I originally wrote; this is intended as a somewhat-minimal example, for ease of understanding. Therefore some programming choices here may seem odd, but they're just due to leftovers from the fully fleshed out script):

from skyfield.api import load, wgs84, EarthSatellite, Topos
import matplotlib.pyplot as plt
import numpy as np

INIT_LAT,INIT_LON, ELEVATION = 40.33825356986421,-72.86106736111111,0

#Starting time to look for transits
START_TIME = [2021, 9, 20, 17, 0,0] #Remember to use UTC!

DURATION = 3 #days to search through

def plot_sat(TLE, timestamp):
    global LATITUDE,LONGITUDE
    sat = EarthSatellite(*TLE)
    ground = Topos(LATITUDE, LONGITUDE, elevation_m = ELEVATION)
    alt,az,dist = (sat - ground).at(timestamp).altaz()

    #Also get where it was a bit ago so we can draw a little line
    a_bit_ago = ts.tt_jd(timestamp.tt -10/86400)

    prev_alt,prev_az,prev_dist = (sat - ground).at(a_bit_ago).altaz()
    trail = axes.plot([prev_az.degrees,az.degrees],[prev_alt.degrees,alt.degrees],color='w')
    dot = plt.Circle((az.degrees,alt.degrees),0.1,color='r')
    axes.add_artist(dot)
def plot_moon(timestamp):
    global LATITUDE,LONGITUDE
    ground = earth + Topos(LATITUDE, LONGITUDE, elevation_m = ELEVATION)
    alt,az,dist = (moon - ground).at(timestamp).altaz()

    moon_angular_radius = np.arcsin(1737.4/dist.km)*180/np.pi #moon radius
    moon_disk = plt.Circle((az.degrees,alt.degrees),moon_angular_radius,color='w')
    axes.add_artist(moon_disk)
def find_closest_approach(tle):
    global LATITUDE,LONGITUDE
    sat = EarthSatellite(*tle)
    observer = Topos(LATITUDE,LONGITUDE, elevation_m = ELEVATION)
    start = ts.utc(*START_TIME)
    end = ts.tt_jd(start.tt + DURATION)
    times_and_events = sat.find_events(observer, start, end)
    passes = []
    last_start = 0
    for i in zip(*times_and_events):
        event_time = i[0]
        event_type = i[1]
        if event_type == 0:
            last_start = event_time
        if event_type == 2:
            passes.append([last_start,event_time])
    closest_moon_dist = 9999
    closest_moon_time = 0
    for p in passes:
        alts = []
        azes = []
        drawtime = p[0]
        while(drawtime.tt < p[1].tt):
            drawtime = ts.tt_jd(drawtime.tt + 1/86400)
            alt,az,_ = (sat - observer).at(drawtime).altaz()
            moon_alt,moon_az,_ = (moon - (earth+observer)).at(drawtime).altaz()
            if moon_alt.degrees < 5:
                break #break out if the moon is crazy low for this pass
            #this is not proper math but good enough heuristic for now
            moon_dist = angular_separation(alt,moon_alt,az,moon_az)
            if moon_dist < closest_moon_dist and alt.degrees > 10 and moon_alt.degrees > 10:
                closest_moon_dist = moon_dist
                closest_moon_time = drawtime
    return closest_moon_time
def angular_separation(alt1,alt2,az1,az2):
    alt1 = alt1.radians
    alt2 = alt2.radians
    az1 = az1.radians
    az2 = az2.radians
    #http://spiff.rit.edu/classes/phys373/lectures/radec/radec.html
    cosprod = np.cos(np.pi/2 - alt1) * np.cos(np.pi/2 - alt2)
    sinprod = np.sin(np.pi/2 - alt1) * np.sin(np.pi/2 - alt2)
    sinprod = sinprod * np.cos(az1 - az2)
    return np.arccos(cosprod + sinprod)
#Take a skyfield time object with floating seconds and round to nearest int second
def round_seconds(skyfield_time):
    time_utc = list(skyfield_time.utc)
    time_utc[5] = round(time_utc[5])
    return ts.utc(*time_utc)
fig = plt.figure()
axes = fig.add_subplot(1,1,1)
axes.set_aspect(1)

ts = load.timescale()
global PLOT_TIME, LATITUDE, LONGITUDE
LATITUDE, LONGITUDE = INIT_LAT, INIT_LON
planets = load('de421.bsp')
earth = planets['earth']
moon = planets['moon']
sun = planets['sun']

sat_tle = [
"1 25544U 98067A   21273.48726090  .00001212  00000-0  30298-4 0  9993",
"2 25544  51.6456 182.0705 0003968  45.0544 125.1915 15.48872009304915"]
raw_closest = find_closest_approach(sat_tle)
PLOT_TIME = round_seconds(raw_closest)
axes.set_facecolor("black")
axes.set_xlim(0,360)
axes.set_ylim(0,90)
plot_moon(PLOT_TIME)
plot_sat(sat_tle,PLOT_TIME)

plt.show()

This generates a plot image, which you can zoom into, and see the red dot of the ISS landing nearly centered on the moon. Therefore, we've determined that at the precise time of [2021,9,23,4,31,59.0], the ISS passes in front of the moon (Well, passed in front of the moon - I've been working on debugging this for a while :) ).

Now, to validate that, I coded up the mathematical version I mentioned. If we punch in that same time, we hope to see the same location that generated our first image: 40.33825356986421,-72.86106736111111.
Here's that code:

from skyfield.api import load,EarthSatellite,wgs84, Topos
from skyfield.positionlib import Geocentric
from numpy.linalg import norm
def normalize(vector):
    return vector / norm(vector)
def find_earth_intersect_point(origin,vector,earth_radius):
    #Step 1: Turn that vector into a unit vector. Still pointing same direction.
    unit_vector_toward_ground = normalize(vec_from_moon_to_sat)
    #Step 2: Do a binary search to find required vector length.
    #Step 2a: First, extend the vector over and over until you punch through the earth.
    #Note: This step is risky if the vector is near-tangent with the earth.
    #Only use when you know the vector is coming in roughly vertical to the surface.
    length_guess = 1
    while(norm(origin + unit_vector_toward_ground * length_guess) > RE):
        print(norm(origin + unit_vector_toward_ground * length_guess))
        length_guess *= 2
    #Once we find that number, it's under the earth. That's an upper bound on length.
    length_upper_bound = length_guess
    #Meanwhile, the lower bound is the previous guess, which is just half as long.
    length_lower_bound = length_guess/2
    #Step 2b: Now that we have upper and lower bounds, binary search to find the
    #exact required length.
    error = 9999
    #Tolerance of 0.1 km should get us quite close.
    while(abs(error) > 0.1):
        #Standard binary search from here, go higher or lower based on where the guess winds up.
        new_guess = (length_lower_bound + length_upper_bound) / 2
        new_vector_end = norm(origin + unit_vector_toward_ground * new_guess)
        if new_vector_end > RE:
            length_lower_bound = new_guess
        else:
            length_upper_bound = new_guess
        error = new_vector_end - RE
    intersection_point = origin + unit_vector_toward_ground * new_guess
    subpoint = wgs84.subpoint(Geocentric(intersection_point,t=t))
    return subpoint.latitude.degrees,subpoint.longitude.degrees

ts = load.timescale()

t = ts.utc(2021,9,23,4,31,59)
RE = norm(wgs84.latlon(40.33825356986421,-72.86106736111111).at(t).position.km)
print(RE)
planets = load('de421.bsp')
earth = planets['earth']
moon = planets['moon']

sat_tle = [
"1 25544U 98067A   21273.48726090  .00001212  00000-0  30298-4 0  9993",
"2 25544  51.6456 182.0705 0003968  45.0544 125.1915 15.48872009304915"]
sat = EarthSatellite(*sat_tle)
moon = moon - earth

#Get vector from the moon to the satellite; extending this vector will make
#it run into the earth.
vec_from_moon_to_sat = (sat - moon).at(t).position.km
lat,long = find_earth_intersect_point(sat.at(t).position.km,vec_from_moon_to_sat,RE)
print(lat,long)

The output of this code is:
40.15023086852985, -72.86124647265113.
Again, compare to the observer location from the first script:
40.33825356986421,-72.86106736111111
Interestingly, the longitude seems to match (to the precision I would expect), while the latitude has a wide offset. If I punch that latitude into the first script, the ISS totally misses the moon. An offset of 0.18 degrees of latitude corresponds to almost exactly 20 kilometers, which is a pretty wide margin to miss.

Now, the best guesses I can make regarding why the answer would be different would be if my earth radius is being treated differently by the two scripts. But the first script takes an elevation of 0, and the second script makes sure to evaluate the earth radius at the known intersection location, so that doesn't seem like it would end up generating errors.

In the end, I'm wondering if I'm going about this all wrong. Maybe I'm taking the wrong approach with my second script? Is there a better way to validate my first script? What's the best way to do what I'm looking to do, to find the vector intersect point? Maybe I'm running into a floating point error? I could see a possibility that my unit vector ends up slightly misaligned, and that that would make a sub-degree error that, multiplied over the hundreds of kilometers up to the ISS, makes for a 20 kilometer error, but any advice on where to hunt would be hugely appreciated.

Thanks for taking the time to read!

[brandon-rhodes]

A first observation: I would expect the elevation of the resulting position (above or below ground) to be fairly small — say, the distance to the Moon in AU. But if I add:

print('Elevation:', subpoint.elevation)

— right before your return, then I get:

Elevation: 6369.14 au

Is this also a larger number than you were expecting, or have I misunderstood what the loop is trying to accomplish?

[bismurphy]

Hmm, that's peculiar. Well, the wgs84.subpoint call is just intended to convert the ECI intersection_point into a latitude and longitude. The intersection_point is a vector of distances in kilometers, and ultimately the magnitude of that vector is 6369.14 kilometers. So my call to Geocentric must be getting interpreted as being in AU rather than kilometers? That's unusual. That might be another source of floating point error then, if we're trying to get the subpoint of a location that is 6369.14 AU away, rather than 6369.14 km away (from the center of the earth).

[bismurphy]

Oh! And if I import AU_KM from Skyfield's built-in Constants, and change my Geocentric call to Geocentric(intersection_point / AU_KM,t=t), then I get 40.34001300403127, -72.86124647265113 which is MUCH closer! Now just need to determine whether I can avoid doing that scaling and de-scaling on my coordinates, to preserve even more precision.

[psbaltar]

I was also trying to work out ISS transits and ran into the problem of find where the line of sight intersects with Earth. However, trying to find where a line intersects an ellipsoid seemed to be really complicated based on what I was finding in math forums. Thanks to a hint from @brandon-rhodes in #801, I was able to code something based on the SPICE surfpt routine. They actually had a pretty straightforward way of handling it. The basic idea was to apply a transformation to treat the problem as finding the intersection of a sphere, then reverse the transformation.

Here's my script which looks for locations at which the ISS will appear to transit Mars. Though, it should be trivial to modify to look for any satellite transiting any body. I get the positions of the target body and the satellite, then "draw" a vector between them. The location where this vector intersects Earth is where the target body and satellite will be in the same apparent position. At the end, I also plot the locations on a map. Though, I still need to work out how to handle the longitude wrapping.

from skyfield.api import load, EarthSatellite, wgs84
from skyfield.framelib import itrs
from skyfield.positionlib import ICRF
from skyfield.units import Distance, Velocity
from skyfield.functions import length_of, dots

from datetime import datetime, timedelta, timezone
import numpy as np
from numpy import array, where, sqrt

import matplotlib.pyplot as plt
import cartopy.crs as ccrs

def vXYZxv3(vXYZ, v3):
  """
  Given one or more vectors in `vXYZ` and single vector `v3`, return their
  component scalar products.

  Equivalent to 'for vec in vXYZ: v = vec*v3'

  This works whether `vXYZ` has the shape ``(3,)``, or whether it is a whole
  array of corresponding x, y, and z coordinates and has shape ``(3, N)``.
  `v3` must be shape ``(3,)``.
  """
  vX = vXYZ[0]*v3[0]
  vY = vXYZ[1]*v3[1]
  vZ = vXYZ[2]*v3[2]
  v = array([vX, vY, vZ])
  return v

def unitize(v):
  """
  Given a vector, return its unit vector.
  
  This works whether `v` has the shape ``(3,)``, or whether it is a whole
  array of corresponding x, y, and z coordinates and has shape ``(3, N)``.
  """
  return v/length_of(v)

def vproj(a, b):
  """
  Given 2 vectors, return the projection of the first onto the second.

  This works whether `a` and `b` each have the shape ``(3,)``, or
  whether they are each whole arrays of corresponding x, y, and z
  coordinates and have shape ``(3, N)``.
  """
  bhat = unitize(b)
  proj = dots(a, bhat)*bhat
  return proj

def vperp(a, b):
  """
  Given 2 vectors, return the component of the first that is perpendicular
  to the second.

  This works whether `a` and `b` each have the shape ``(3,)``, or
  whether they are each whole arrays of corresponding x, y, and z
  coordinates and have shape ``(3, N)``.
  """
  proj = vproj(a, b)
  perp = a - proj
  return perp

def ray_intersect_ellipsoid(positn, u, a, b, c):
  """
  Based on SPICE surfpt
  https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/surfpt_c.html

  NOTE: Unlike surfpt, `positn` is assumed to be outside of the ellipsoid.
  For cases where `positn` is on or in the ellipsoid, results are invalid.

  Parameters:
    `positn` is the position of an observer with respect to the center of an
    ellipsoid. The vector is expressed in a body-fixed reference frame.
    The semi-axes of the ellipsoid are aligned with the x, y, and z-axes
    of the body-fixed frame.

    `u` is the pointing vector emanating from the observer. `u` does not need
    to be a unit vector.

    `a`, `b`, and `c` are the ellipsoid semi-axes along the X, Y, and Z axes,
    respectively.

  Returns:
    `point` is the point on the ellipsoid pointed to by `u`

    `found` is a boolean indicating if `u` is pointing at the ellipsoid

  This function works whether `positn` and `u` each have the shape ``(3,)``,
  or whether they are each whole arrays of corresponding x, y, and z
  coordinates and have shape ``(3, N)``. `a`, `b`, and `c` are scalars.
  """
  elshape = array([a,b,c])
  invelshape = 1/elshape

  '''
  General procedure:
  1. Apply a linear transformation to the observer position,
     direction vector, and ellipsoid so that the problem is reduced to
     finding the intersection of the vector and the unit sphere.
     The ellipsoid itself doesn't need to be transformed since it
     will be a unit sphere.
  2. Find the intersection of the direction vector and the unit sphere.
  3. Reverse the transformation.
  '''

  x = vXYZxv3(u, invelshape) #u/elshape
  y = vXYZxv3(positn, invelshape) #positn/elshape

  yproj = vproj(y,x)
  p = y - yproj #same as vperp(y,x)

  ymag = length_of(y)
  pmag = length_of(p)

  ux = unitize(x)

  # TODO handle ymag<=1 cases (positn is on/inside ellipsoid)
  '''
  SPICE surfpt checks for two no interesect conditions to return immediately.
  But we will continue with the calculations anyway so this function can
  handle an array of vectors.
  '''
  # no intersection if p is outside of sphere
  found = where(pmag>1, False, True)
  # no intersection if direction vector is pointing away from sphere
  found = where(dots(yproj, x)>0, False, found)

  # before sqrt, clip negative values of d to 0 so that the calculations
  # can be done even for invalid values
  d = 1-pmag**2
  dclip = array([d, d*0]).max(axis=0) # clip negative values
  scale = sqrt(dclip)
  point = p - scale*ux

  # reverse the transformation
  point = vXYZxv3(point, elshape) # point*elshape

  # set point to zero if no intersection to match behavior of SPICE surfpt
  mask = where(found, 1, 0)
  point *= mask

  return point, found


# ---- initialize everything ----

sat_name = 'ISS (ZARYA)'
tle = [
'1 25544U 98067A   22290.76828946  .00015865  00000+0  28397-3 0  9991',
'2 25544  51.6422  89.2789 0003316 312.7838 181.7649 15.50140930364197' ]
sat = EarthSatellite(*tle, name=sat_name)

# need to load both main and planet-specific ephemerides
# https://github.com/skyfielders/python-skyfield/issues/145
ephemeris = load('de421.bsp')
mars_ephemeris = load('mar097.bsp')
mars_ephemeris.segments.extend(ephemeris.segments)

earth = ephemeris['EARTH']
target = mars_ephemeris['MARS']

earth_shape = [6378.1366, 6378.1366, 6356.7519] # from pck00010.tpc

ts = load.timescale()


# ---- the good stuff starts here ----

# build a `Time` object from a datetime arange array
t_range = np.arange( datetime(2022,10,17),  # start
                     datetime(2022,10,18),  # stop
                     timedelta(minutes=1)   # timestep
                       ).astype(datetime)
# `from_datetimes()` requires tzinfo to be set
for i in range(t_range.shape[0]):
  t_range[i] = t_range[i].replace(tzinfo=timezone.utc)
t = ts.from_datetimes(t_range)

# need to use Earth body-fixed coordinates (ITRS) so that the axes
# match the `earth_shape` model
earth_target = earth.at(t).observe(target).frame_xyz(itrs).km
earth_sat = sat.at(t).frame_xyz(itrs).km
target_sat = earth_sat - earth_target

# draw a ray from `earth_sat` in the direction of `target_sat` and find
# were it intersects the ellipsoid defined by `earth_shape`
point, found = ray_intersect_ellipsoid(earth_sat, target_sat, *earth_shape)

# convert the intersection point back to ICRF
point_pos = Distance(km=point)
point_vel = Velocity(km_per_s=point*0) # dummy velocity to build position
point_icrf = ICRF.from_time_and_frame_vectors(t, itrs, point_pos, point_vel)
point_icrf.center = 399

# get the geographic coordinates of the intersection point
geopos = wgs84.geographic_position_of(point_icrf)
lat = geopos.latitude.degrees
lon = geopos.longitude.degrees
elev = geopos.elevation.km # just a sanity check, should be 0 (on the ground)


# ---- print stuff ----
for i in range(t.shape[0]):
  if found[i]:
    print(round(lat[i],6), round(lon[i],6), round(elev[i],6), t[i].utc_iso())
    
    # angle between target and sat should be 0, exit if >0.36"
    loc = wgs84.latlon(lat[i], lon[i], elevation_m=elev[i]*1000)
    s = (sat - loc).at(t[i])
    pos_target = (earth+loc).at(t[i]).observe(target)
    angle = pos_target.separation_from(s)
    print(angle)
    if angle.degrees > 0.0001: print('separation not 0'); exit()


# ---- plot stuff ----

# replace not found points with `np.nan` so they don't get plotted
# should also figure out how to handle lon wrapping from 360 to 0
lat = np.where(lat==0, np.nan, lat)
lon = np.where(lon==0, np.nan, lon)

fig = plt.figure(figsize=(10, 5))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

ax.stock_img()

ax.plot(lon, lat)

plt.show()