An Earth Occultation calculation for Satellites

[steveehlert]

I am trying to figure out if a satellite's view of a particular RA/DEC target is occulted by the Earth over a large range of times. I have some code which I think should do this, but I am not sure the results make perfect sense. I am also uncertain as to whether or not the coordinates I am choosing for each part of the calculation make sense. Here is the code snippet that calculates the minimum distance between the Earth center and the line of sight between

        target_pos = skyfield.positionlib.position_from_radec(target_ra, target_dec)
        target_ecliptic = target_pos.ecliptic_position().au
        earth_ecliptic = Earth.at(ut).ecliptic_position().au
        sat_ecliptic = satellite.at(ut).ecliptic_position().au

        limit_dist = (EARTH_RADIUS + limiting_altitude) / AU_TO_KM
        #print(limit_dist)
        dist_list = []
        n_times = numpy.shape(ut)[0]
        for j in range(n_times):
            this_sat = sat_ecliptic[:,j]
            this_earth = earth_ecliptic[:,j]
            # this_sat is satellite position relative to Earth, whereas this_earth
            # and target_ecliptic are positions in ecliptic coordinates.
            # This is why we add this_sat and this_earth below.

            diff_sat_target = this_sat + this_earth  - target_ecliptic
            diff_target_earth = target_ecliptic - this_earth

            cross = numpy.cross(diff_sat_target,diff_target_earth)

            lcross = numpy.sqrt(numpy.dot(cross,cross))

            dist_list.append(lcross / numpy.sqrt(numpy.dot(diff_sat_target,diff_sat_target)))
        dist_arr = numpy.array(dist_list)

with the elements of dist_arr being tested against the limit_dist assumed above. The behavior I am finding, however, does not necessarily jive with my expectations given the satellite's nearly equatorial orbit. So three questions: 1) Are all of these positions in a consistent coordinate system? 2) If not, what changes need to be made to put them in one? 3) Would this be something you would be interested in integrating into skyfield moving forward?

[brandon-rhodes]

The three positions you generate with ecliptic_position() are indeed all in the same coordinate system. I cannot quite tell for sure where their origins are, though, because I can’t see how the objects are built (in particular the Earth object).

I am also confused by the comment “this_sat is satellite position relative to Earth, whereas this_earth and target_ecliptic are positions in ecliptic coordinates” — isn't this_sat also in ecliptic coordinates?

But, yes, Skyfield would benefit from built-in math to determine intersections of bodies with lines-of-sight. Some work in that direction is currently underway by @jacoffey3 in #369 for the specific case of whether a satellite is in sunlight or not, and the same math should work for whether any other RA/dec in the sky is occluded by the Earth, once we have that code finished. The related issue is #327.

I wonder what the API should look like for asking whether one position is visible from another, given one or more possible bodies that might be in the way? Feel free to share if you have any ideas!

[steveehlert]

Here are the answers to your questions -

The satellite object is created as skyfield.api.EarthSatellite with a TLE The planets are created through skyfield.api.load_file using a trimmed version of de430 The RA/DEC position is created through skyfield.positionlib.position_from_radec as described above. Making sure the origins are consistent certainly matters!

That actually brings up another question - how is an RA/DEC (a 2-d direction) converted into a 3-d ecliptic position in AU? In other words, what distance is assumed for that conversion? I can imagine that the results may be sensitive to that if I am not careful.

As for the API, I would probably suggest something that takes the time, target position, and occulting body as arguments and only checks one body at a time - in practice there is usually only one body that may be occulting (primarily the body the satellite orbits) and the users will have a pretty good idea which bodies they want to check. So they can simply run it multiple times with different bodies (Earth, Moon, etc) as input. An example below:

EDIT: This is the code I have written for my own fork of skyfield to include this in sgp4lib.py

    def is_target_occulted(self, ut, target_pos, test_body_name, ephemfile = 'de421.bsp', dist_limit = ERAD/1000. + 200.):
        """Returns a True/False as to whether or not the 
        line of sight between the satellite position and a target position of interest
        is occulted by a planetary body. Originally written to account for 
        occultations of particular RA/DEC positions by the Earth. User defines the limiting distance for occultation
        and the ephemeris file that has the target body's ephemeris included.
        Requred inputs are the ut time of observations, the target_pos Skyfield position, and 
        the name of the (possibly) occulting body

        """

        ephemdata = apiload(ephemfile)

        sat_origin_pos = ephemdata['earth'].at(ut)
        test_body = ephemdata[test_body_name]

        sat_ecliptic = sat_origin_pos.ecliptic_position().au +   self.at(ut).ecliptic_position().au
        test_body_ecliptic = test_body.at(ut).ecliptic_position().au 
        target_ecliptic  = target_pos.ecliptic_position().au
        diff_sat_target = sat_ecliptic  - target_ecliptic
        diff_target_body = target_ecliptic - test_body_ecliptic
        #Using cross product formula of
        #https://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html
        #to determine minimum perpendicular distance between line of sight
        #and a third position (the occulting body's center)
        cross = numpy.cross(diff_sat_target,diff_target_body)
        lcross = numpy.sqrt(numpy.dot(cross,cross))
        perp_dist = lcross / numpy.sqrt(numpy.dot(diff_sat_target,diff_sat_target))
        l_sat_target = numpy.sqrt(numpy.dot(diff_sat_target,diff_sat_target))
        l_target_body = numpy.sqrt( numpy.dot(diff_target_body,diff_target_body))
        #Note I am not sure if there is a function for a planetary body's radius like I allude to 
        #below. Note that the return is True of occultation occurs and False if the LOS is clear
        #The alt_limit argument lets you restrict to a region larger than the planetary body's radius itself   
        return ( (perp_dist < (dist_limit / AU_KM ) ) and (l_sat_target > l_target_body) )

[brandon-rhodes]

The satellite object is created as skyfield.api.EarthSatellite with a TLE The planets are created through skyfield.api.load_file using a trimmed version of de430 The RA/DEC position is created through skyfield.positionlib.position_from_radec as described above. Making sure the origins are consistent certainly matters!

Thank you for these details! I think they unravel the mystery.

First: my guess is that you only need two vectors here. One vector will be satellite → target and can be generated from the target’s RA/dec. The other vector is the satellite → Earth vector, which is simply the negation of the satellite’s .at() position, since a satellite vector’s center is the Earth. Happily, this means you do not need an ephemeris at all; the question of whether the Earth occludes a target on the celestial sphere involves only the local question of where the Earth and satellite are relative to each other, regardless of where the two of them are in the Solar System as a whole.

Second: you are quite correct that the problem with your actual approach, that tries to base both vectors at the center of the Solar System, is the length of the default position_from_radec() vector, which is (drum roll) 1 au. Which means the position you are generating is typically going to be very distant from the RA/dec you intended to place it at (assuming the target is way out on the celestial sphere? I am not sure exactly what kind of target you are contemplating, but it sounds like a star or galaxy or radio source).

With the existing position_from_radec() and your existing code, you will want to instead provide a much larger distance.

But, on reflection, I think that position_from_radec() is very badly designed. For one thing, the distance argument includes no units in its name, which is strictly contrary to how I try designing Skyfield routines. For another, the 1 au value is a poor choice on my part. I probably hoped the value of exactly 1.0 looked fake enough that users would immediately recognize it as a placeholder and supply something better. But, really, how often are users going to think to check the length?

So I have just landed 52cd9a4bfbdc7eee19376a7acab275810c1ce969 which provides a better-designed position_of_radec() for you to try out. Could you try:

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

and let me know how the new routine works for you? It not only carefully documents the distance argument, which the older routine failed to do, but has an extremely large distance as the default so that if users make no adjustments themselves, they get a position so far outside the Solar System that its RA/dec should be the same from anywhere inside, eliminating any problems with where they choose to center the vector.

As for the API, I would probably suggest something that takes the time, target position, and occulting body as arguments and only checks one body at a time - in practice there is usually only one body that may be occulting (primarily the body the satellite orbits) and the users will have a pretty good idea which bodies they want to check. So they can simply run it multiple times with different bodies (Earth, Moon, etc) as input. An example below: …

Thank you for the sample code! I will take a look at it when I have a chance. My guess is that we will ultimately want to adapt the JPL code for ellipsoidal bodies instead of approximating all objects as spheres, but I can start from here when working through what the arguments will need to be to cover a case like yours.

[steveehlert]

Hi Brandon!

I did a quick test of my previous code with the new version filled in and it works very well! I noticed independently since I wrote my last comment what the distance was and realized it was causing me trouble. I am indeed trying to observe sources well beyond the Solar System.

I have one suggestion based on the changes and some other testing I have done in the past few days - can you actually make the default distance a little bit smaller? I think a Gigaparsec may actually be too far for the floating point precision of the distance calculations. I was noticing that my occultation fraction across many times suddenly increases to about 100% at this distance. A smaller distance of something like 1e9 AU (of order a kpc) instead of a Gpc gives an occulation fraction of 39%, which is what I expect from independent calculations.

And it is definitely true that elliptical bodies would ultimately be more accurate in an absolute sense, but I have not yet bothered to figure out how large of an effect that is.

EDIT: One thing I would like to do with the occultation is allow the limiting distance of occultation to be set by the user - you'll notice that I wanted 200 km above the surface of the Earth in my calculation above and I can imagine other satellites/observatories having different requirements. For my purposes a spherical Earth is sufficient, so I used the simpler spherical calculation.

Thank you for updating this so quickly!

[brandon-rhodes]

I have one suggestion based on the changes and some other testing I have done in the past few days - can you actually make the default distance a little bit smaller? I think a Gigaparsec may actually be too far for the floating point precision of the distance calculations. I was noticing that my occultation fraction across many times suddenly increases to about 100% at this distance. A smaller distance of something like 1e9 AU (of order a kpc) instead of a Gpc gives an occulation fraction of 39%, which is what I expect from independent calculations.

I chose the gigaparsec because the United States Naval Observatory’s NOVAS library, which Skyfield's internal calculations are based on, uses that value for stars (and other objects) whose measured parallax is zero. It would be nice to remain consistent with them where possible. When I have the chance, I'll try your code above (thanks for providing it!) with one of these gigaparsec vectors to see if I can uncover the problem with the math; I hadn’t thought that cross products or dot products would care.

Edit: (And, the gigaparsec is already used internally in the Skyfield object for stars and star catalogs.)

[brandon-rhodes]

After a few tries, I am not having success in making your code misbehave — so to investigate further I would need a complete script that, if pasted into a Python file, would print out a bad result so that I could get to work on those specific numbers.

But, in the meantime — I was surprised to see the vectors diff_sat_target and diff_target_body as I would have expected to see diff_sat_target and diff_sat_body. Could it be that you are losing precision because diff_target_body combines a very large and a very small vector?

[steveehlert]

I don't believe I am in a position to share the full script (it isn't for public use) but I will do some more testing on my own and see if I can find another cause. And I will certainly check to see if changing the subtraction around as you propose makes a difference.

[brandon-rhodes]

Understood! I meant “a complete script” only in the sense of “a standalone Python script that reduces your test case to the minimum code necessary to demonstrate the problem” — I understand completely that you will not want to share your full script or application.

Let me know if your own investigations turn anything up, I’ll be interested to hear.

[steveehlert]

So I figured out what the problem was - it is a floating point precision issue but not related to the minimum perpendicular distance calculation but instead of figuring out if the Earth is between the satellite and target along the line of sight (the l_sat_target > l_target_body test in the code above). I wrote a new version that should be more robust to this below.

Forgive me for not having a working version with a test case just yet. This was a bit faster to put together.

def is_target_occulted(self, ut, target_pos, test_body_name, ephemfile = 'de421.bsp', dist_limit = ERAD/1000. + 200.):

 """Returns a True/False as to whether or not the 
        line of sight between the satellite position and a target position of interest
        is occulted by a planetary body. Originally written to account for 
        occultations of particular RA/DEC positions by the Earth. User defines the limiting distance for occultation
        and the ephemeris file that has the target body's ephemeris included.
        Requred inputs are the ut time of observations, the target_pos Skyfield position, and 
        the name of the (possibly) occulting body

        """

        ephemdata = apiload(ephemfile)

        test_body = ephemdata[test_body_name]
        sat_origin_pos = test_body.at(ut)
        sat_ecliptic = sat_origin_pos.ecliptic_position().au +   self.at(ut).ecliptic_position().au
        test_body_ecliptic = test_body.at(ut).ecliptic_position().au 
        target_ecliptic  = target_pos.ecliptic_position().au

        diff_body_sat = test_body_ecliptic - sat_ecliptic
        l_body_sat = sqrt(dot(diff_body_sat,diff_body_sat))

        diff_body_targ = test_body_ecliptic - target_ecliptic
        l_body_targ = sqrt(dot(diff_body_targ,diff_body_targ))

        diff_targ_sat = target_ecliptic - sat_ecliptic
        l_targ_sat = sqrt(dot(diff_targ_sat, diff_targ_sat))

        #This is the (x0-x1) X ( x0-x2) version of the min dist formula
        crossprod = cross(diff_body_sat,diff_body_targ)

        lcross = sqrt(dot(crossprod,crossprod))

        dist_perp = lcross / l_targ_sat

        #Floating point precision is an issue for this calculation since the
        # diff_body_sat << diff_body_targ. When target is at a default
        #distance of 1 Gpc
        #So we actually calculate the value
        # of the parameter t where it reaches minimum distance.

        ttt = -1.* dot( -1. * diff_body_sat, diff_targ_sat) / l_targ_sat**2

        #Note I am not sure if there is a function for a planetary body's radius like I allude to 
        #below. Note that the return is True of occultation occurs and False if the LOS is clear
        #Yhr dist_limit argument allows for distances larger than the planetary radius itself to be used
        return ( (dist_perp < (dist_limit / AU_KM ) ) and (ttt > 0) )

[brandon-rhodes]

Thanks for the update!

Do you have sample arguments I could call is_target_occulted() with to see the misbehavior in the old version and the success in the new version? I don't yet have any intuition as to why the first floating point comparison was fragile, and seeing some real numbers flow through the routine would help me develop a sense for how they're behaving.

[steveehlert]

Here is a script that runs through the Earth occultation calculation for some representative times. It agrees with my "good" results. It has the exact code as above added to sgp4lib.py

from skyfield.api import EarthSatellite, load, position_of_radec
from skyfield.nutationlib import iau2000b
import numpy as np
import tqdm

de421 = load('de421.bsp')
earth = de421['earth']
ts = load.timescale()

this_sat = EarthSatellite('1 00000U 21  0A   20146.43738583  .00000000  00000-0  00000-0 0  0000',
                          '2 00000   0.2000 000.0000 0000000   0.0000 000.0000 14.88336848949720', 'Test_600kmAlt', ts)

#100K times between 2022 and 2024

n_steps = 10000
jdrange = np.linspace(2459690.5,2460055.75,n_steps)

utrange = ts.ut1_jd(jdrange)
utrange._nutation_angles = iau2000b(utrange.tt)

#Crab Nebula
target_ra = 83.63308333
target_dec = 22.0145

#My RA is in degrees, hence the division by 15. 
target_pos = position_of_radec(target_ra / 15., target_dec)

#Currently not optimized for speed
bool_list = []
for this_ut in tqdm.tqdm(utrange):
    this_bool = this_sat.is_target_occulted(this_ut, target_pos, earth)
    bool_list.append(this_bool)

bool_mask = np.array(bool_list)

n_occulted = np.shape(np.where(bool_mask))[1]

frac_occulted = float(n_occulted)/float(n_steps)

print("A total of %.2f percent of the time steps were occulted" %(frac_occulted * 100.))

The target_occulted function got an update as well to move the ephemeris and body calculations out of the function

    def is_target_occulted(self, ut, target_pos, test_body, dist_limit = ERAD/1000. + 200.):
        """Returns a True/False as to whether or not the 
        line of sight between the satellite position and a target position of interest
        is occulted by a planetary body. Originally written to account for 
        occultations of particular RA/DEC positions by the Earth. User defines the limiting distance for occultation
        and the ephemeris file that has the target body's ephemeris included.
        Requred inputs are the ut time of observations, the target_pos Skyfield position, and 
        the name of the (possibly) occulting body

        """

        sat_origin_pos = test_body.at(ut)
        sat_ecliptic = sat_origin_pos.ecliptic_position().au +   self.at(ut).ecliptic_position().au
        test_body_ecliptic = test_body.at(ut).ecliptic_position().au 
        target_ecliptic  = target_pos.ecliptic_position().au

        diff_body_sat = test_body_ecliptic - sat_ecliptic
        l_body_sat = sqrt(dot(diff_body_sat,diff_body_sat))

        diff_body_targ = test_body_ecliptic - target_ecliptic
        l_body_targ = sqrt(dot(diff_body_targ,diff_body_targ))

        diff_targ_sat = target_ecliptic - sat_ecliptic
        l_targ_sat = sqrt(dot(diff_targ_sat, diff_targ_sat))

        #This is the (x0-x1) X ( x0-x2) version of the min dist formula
        crossprod = cross(diff_body_sat,diff_body_targ)

        lcross = sqrt(dot(crossprod,crossprod))

        dist_perp = lcross / l_targ_sat

        #Floating point precision is an issue for this calculation since the
        # diff_body_sat << diff_body_targ. When target is at a default
        #distance of 1 Gpc
        #So we actually calculate the value
        # of the parameter t where it reaches minimum distance.

        ttt = -1.* dot( -1. * diff_body_sat, diff_targ_sat) / l_targ_sat**2

        #Note I am not sure if there is a function for a planetary body's radius like I allude to 
        #below. Note that the return is True of occultation occurs and False if the LOS is clear
        #Yhr dist_limit argument allows for distances larger than the planetary radius itself to be used
        return ( (dist_perp < (dist_limit / AU_KM ) ) and (ttt > 0) )

[brandon-rhodes]

Thanks very much for the working code! Having run it several times, my first question is whether you can replace l_sat_target > l_target_body in your original code with l_body_sat > l_target_body? It would not suffer from any floating point limitation, and it would be intuitively what I would expect to type: “if the target is further away from the satellite than the body is from the satellite, then the satellite might see the body in front of the target.”

A more exact answer would involve the distance along the sat→target line to the first intersection with the sphere; if the target is more distance than that intersection, then it's hidden by the body's surface. I'll see what I can work out for a solution to that question, and we can compare the math with yours.

[brandon-rhodes]

All right, I'm back following some work with my yellow pad.

I propose that what you want is a routine that given the vector sat→target and the vector sat→Earth center and the radius of Earth, returns the distance from sat to the nearest point on the Earth sphere. If that distance is non-NaN (meaning, there is indeed an intersection) and less than the length of the sat→target vector (otherwise the target is in front of Earth), then the target is occluded.

Using a Google result:

http://paulbourke.net/geometry/circlesphere/index.html#linesphere

— and some simplification and the some NumPy-ification, I arrived at:

from skyfield.functions import length_of
import numpy as np

def intersect_line_with_sphere(line_endpoint, sphere_center, sphere_radius):
    """Compute distance to the nearest intersection of a line with a sphere.

    line_endpoint: NumPy array giving a point on a line starting at the origin.
    sphere_origin: NumPy array giving the center of the sphere.
    sphere_radius: Float giving the sphere radius.

    """
    unit = line_endpoint / length_of(line_endpoint)
    b = -2.0 * (unit * sphere_center).sum()
    c = (sphere_center * sphere_center).sum() - r * r

    discriminant = b * b - 4 * c
    return np.nan if discriminant < 0 else (-b - np.sqrt(discriminant)) / 2.0

target = np.array([1000, 0.0, 0.0])
ctr = np.array([2.0, 0.99999, 0.0])  # circle barely touching line
r = 1.0

print(intersect_line_with_sphere(target, ctr, r))

Could you verify whether this routine works for you? Since it reduces the sat→target vector immediately to a unit vector as part of the necessary math, I think you will find it immune to any scale differences between the target’s distance and the Earth’s.

[steveehlert]

Okay, I will test this soon. However, just to be absolutely certain I understand the coordinates based on what is written, line_endpoint in my case should be (target -sat) in ecliptic coordinates, and sphere_center should be (earth - sat) in ecliptic coordinates? Assuming that satellite position in barycentric ecliptic coordinates is the origin now.

[brandon-rhodes]

Yes, exactly! Or ra-dec, if you prefer to avoid the expense of the rotation into ecliptic coordinates; any x/y/z coordinate system at which the satellite is at the origin.

[steveehlert]

Okay, first test complete. I got the same result (statistically) between your code and mine once I put the coordinates in the correct system. I had to make one small change to the return statement:

return np.nan if discriminant < 0 else (max(-b - np.sqrt(discriminant), -b + np.sqrt(discriminant)) / 2.0)

to give the maximum value of the two points where it intersects. From reading the link you sent, we only care (in this case) about the t > 0 case (i.e. we intersect the Earth on the way to the target) since t < 0 occurs "behind" the satellite.

I still need to test the agreement time bin by time bin.

[brandon-rhodes]

But in the case where the ray from the satellite intersections the earth in front of the satellite twice, isn't it the min where the satellite will lose visibility? Maybe we want "min, if min is positive, else max"? I'm stepping away from the keyboard so it will only be a bit later that I can work out if I'm correct or visualizing the math wrong.

[steveehlert]

If you are interested in identifying what point on the Earth gives you the first "impact" then yes you will want to keep min instead of max. I am only interested in knowing whether either intersection point was positive for this particular project, so I kept max.

Also, I have now done the time bin by time bin comparison, which provides perfect agreement on 10,000 time bins.

[steveehlert]

Finally, I will let you know that all of the files used in this conversation are in my own fork of skyfield (https://github.com/steveehlert/python-skyfield) so feel free to borrow from there from now on.

[brandon-rhodes]

I have just released a new version of Skyfield which offers an is-in-sunlight routine powered by the above math:

https://rhodesmill.org/skyfield/earth-satellites.html#finding-when-a-satellite-is-in-sunlight

I suppose that next I'll need to generalize the example to handle the situation in this issue: whether the Earth occults an arbitrary object, rather than the Sun.

[steveehlert]

I looked over the example and it looks great! I think we have all of the math figured out so I suppose it is just a matter of putting the proper functionalities in at this point. Let me know if you need any help with that!

[brandon-rhodes]

@steveehlert — I found a few minutes this morning to try adding a method! Is the name is_behind_earth() clear enough or too obscure? The names is_blocked_by_earth() or is_blocked_by_earths_disc() or is_occluded_by_earth() or is_earth_occluded() might have been clearer, but would also take up more space on a line and would take longer to type.

Anyway, if you get the chance to try out the new method, I would be interested in hearing whether it agrees with your own code. Thanks!