interpolating    mapping    searching

The Backtrack Orbit Search Algorithm Detailed Description

Introduction:

Definitions:

The distance from the pole to the satellite as it passes closest to the pole is the declination  This is (inclination - 90).

The amount of time the satellite takes to complete one orbit is the period.

The width of the field of view of the sensor is the swath width .

The highest point achieved by a given orbit is the inflection point .  This is the point at which the orbit goes from ascending to descending.

The highest latitude achieved by an orbit, the latitude of the inflection point, is the inflection latitude.  This is (90 - declination).

The lowest latitude covered by a swath, when the satellite is at the point of inflection, is the minimum inflection latitude.  It is the latitude of the left swath edge when the satellite is at the point of inflection.

The highest latitude covered by a swath, when the satellite is at the point of inflection, is the maximum inflection latitude.  It is the latitude of the right swath edge at the point of inflection.  We make this a circular latitude - it could be greater than 90 if the swath covers the pole.  So max_inflection_lat = (inflection_lat + (inflection_lat - min_inflection_lat)).

The latitude above which all orbits cover the area is the total coverage latitude.  Not all sensors have this.  Sensors with a narrow swath do not include the pole.  So total_coverage_lat = min(90, (180 - max_inflection_lat)

The latitude above which no orbits cover the area is the maximum coverage latitude.  Not all sensor have this.  Sensors with a wide swath always include the pole.  So max_coverage_lat = min(90, max_inflection_lat)
 
 

Step by step:

0) Break the perimeter into sample points.

The algorithm works on points so given an area create a set of points along the perimeter that are less than width/2 apart.  The points need to be close enough so the resulting crossing ranges will overlap and hence be mergable into one searchable crossing range.  Only perimeter points are needed because that's where the extremes of the area are.  I. e. any swath that includes an interior point of the area must also include a perimeter point.  Steps 1-6 are repeated for each sample point then the crossing ranges are merged prior to the search step.

1) Find the orbit with nadir crossing the given point.

Let the point of interest be (lat_p, lon_p).  We seek the ascending equator crossing point of nadir at (0, lon_n )

First consider a static (non-rotating) sphere.  Given (lat_p , lon_p) we can construct a great circle.  Because we know the inclination of the orbit we know the maximum latitude of the orbit so we can construct the great circle that passes through (lat_p, lon _p) and has it's inflection point at (lat_inf, lon_inf ) .  That great circle and the longitude line of the point intersect at the point. Using those two great circles, the equator, and the longitude segment from the north pole to the inflection point (lat_inf, lon _inf ) we can create two spherical triangles.  We do not know the longitude of the inflection point but we do know the latitude because we know the inclination of the orbit.  That latitude is (90-declination) where declination = (inclination-90) So we know the length of the great circle arc along lon_inf that connects the inflection point to the north pole.

Given a spherical triangle with sides alpha, beta, and gamma and opposite interior angles Alpha, Beta, and Gamma. the Law of sines for spherical triangles is:

sin(alpha)/sin(Alpha) = sin(beta)/sin(Beta) = sin(gamma)/sin(Gamma)

Let theta be the angle between the meridian of the point and the orbital great circle.  Then by the Law of sines

sin(90 - lat_inf) / sin(theta) = sin(90 - lat_p) / sin(90)   so
sin(theta) = sin(90 - lat_inf) / sin(90 - lat_p )

Similarly

sin(lon_n - lon_p)/ sin(theta) = sin(lat_p ) / sin(lat_inf)  so
sin(lon_n - lon_p) = sin(theta) *sin(lat_p ) / sin(lat_inf)  so
lon_n = asin(sin(theta) *sin(lat_p) / sin(lat _inf ))  + lon_p

Similarly we can find the length of the arc between (0, lon_n ) and (lat_p, lon_p)

sin(L_arc) / sin(90) = sin(lat_p) / sin(lat_inf )  so
L_arc = asin(sin(lat_p) / sin(lat_inf ))

The true equatorial crossing can be found by adjusting for the rotation of the Earth during the time it takes the satellite to travel the distance of L_arc.

RealLon_n = lon_n + [ period * L_arc* ROTATION_RATE / circumference]

Eventually we'll want to adjust everything by that amount - but for now we still need to use the static Earth.
 

2) Convert to Cartesian 3-space and find the orbital plane.

These two points (lat_p, lon_p) and (lat_n , lon_n) lie on the Great Circle that is the orbit.  Combined with the origin they define a plane.

Convert to Cartesian

x = r * cos(lon) * cos(lat);
y = r * sin(lon) * cos(lat);
z = r * sin(lat);

The plane of the Great Circle is defined by:

(y_pz_n - y_nz_p)x - (x_p z_n - x_nz_p)y + (x_py _n - x_ny_p)z = 0

Let

a = (y_pz_n - y_nz_p)
b =  -(x_pz_n - x_nz_p )
c = (x_py_n - x_ny_p)

So the plane of the orbit on a static sphere is: ax + by + cz = 0

3) Find points on the swath edges:

Given the swath width find points (lat_edge, lon _edge ) on each edge by going half that distance perpendicular to the ground track.  This only works for downward looking sensors.  Forward looking and backward looking sensors can also use this method by creating avirtual satellite that is downward looking.  For side viewing sensors simply dividing the swath width in half won't work.  Instead the swath width can be broken apart into left and right distances from nadir which requires one addtional parameter in the database table that contains the orbit information.

For the left edge of the swath we need to find the point (lat_left , lon_left) that is distance (width/2) from the point of interest (lat_p , lon_p) along the heading (180 +  theta).

lat_left =  asin( (sin(lat_p)*cos(width/2)) + (cos(lat_p)*sin(width/2)*sin(180+theta)) )

lon_left = lon_p - Acos( (cos(width/2) - sin(lat _p)*sin(lat_left)) / (cos(lat_p)*cos(lat_left )) )

For the right edge of the swath we need to find the point (lat_right , lon_right) that is distance (width/2) from the point of interest (lat_p , lon_p) along the heading (theta).

lat_right =  asin( (sin(lat_p)*cos(width/2)) + (cos(lat_p)*sin(width/2)*sin(theta))
lon_right = lon_p + Acos( (cos(width/2) - sin(lat _p)*sin(lat_right)) / (cos(lat_p)*cos(lat _right )) )
 

4) Find the planes of the swath edges:

Convert the edge point to cartesian (x_edge, y_edge , z_edge).

On a static Earth the swath edge is a small circle parallel to the Great Circle defined by nadir.  In cartesian space that defines a plane parallel to the plane ax+by+cz=0

That plane is:  ax + by + cz = d  where  d = (ax_edge + by_edge + cz_edge)

and the equation for the plane of the small circle which is the other edge of the swath is:

ax + by + cz = -d

5) Find the intersects with the given latitude:

We want the point where each small circle intersects the latitude of the given point, and the sphere.  That is:
ax + by + cz = d    : swath edge
x² +y² + z² = r²     : sphere
z = z_p                    : lat = lat_p

Three equations, three unknowns.  We're keeping the latitude constant so we only care about the longitude which is arctan(y/x).  Substituting z = z_p and rearranging we get:

ax + by + cz_p -d = 0  and
x² -y² + z_p² - r² = 0

In[11]:= Solve[{(x2) + (y2) - (r2) + (z2) == 0, 
         (ax) + (by) + (cz) - d == 0}, 
         {x, y}]

Out[11]= {
 
>   {x -> (d - (b2d/(a2+b2)) - cz + (b2cz/(a2+b2)) - 
 
>         (b * Sqrt[-(a2d2) + a4r2 + a2b2r2 + 2a2cdz - a4z2  - 
 
>              a2b2z2 - a2c2z2]) / (a2 + b2)) / a, 

>     y -> (bd - bcz + Sqrt[-(a2d2) + a4r2 + a2b2r2 + 2a2cdz - 
 
>           a4z2 - a2b2z2 - a2c2z2]) / (a2 + b2)}, 
 
       
>    {x -> (d - (b2d/(a2+b2)) - cz + (b2cz/(a2+b2)) + 
 
>         (b * Sqrt[-(a2d2) + a4r2 + a2b2r2 + 2a2cdz - a4z2 - 
 
>              a2b2z2  - a2c2z2 ]) / (a2 + b2 )) / a, 
 
>     y -> (bd - bcz - Sqrt[-(a2d2) + a4r2 + a2b2r2 + 2a2cdz - 
 
>           a4z2  - a2b2z2  - a2c2z2]) / (a2 + b2 )}}

P = cz_p -d   (Constants in the Plane equation)
S = z_p² - r² (Constants in the Sphere equation)
 

In[12]:= Solve[{(x2) + (y2) + S == 0, ax + by + P == 0},{x,y}]

                    
Out[12]= {
                  b2P       b * Sqrt[-a4S - a2b2S - a2P2 ]
           -P + ------- - ----------------------------------
                 a2 + b2                 a2 + b2
    { x -> ------------------------------------------------, 
                                  a
 
           -bP + Sqrt[-(a4S) - a2b2S - a2P2 ]
>     y -> ----------------------------------------}, 
                           a2 + b2

                  b2P       b * Sqrt[-(a4S) - a2b2S - a2P2 ]
           -P + ------- + -----------------------------------
                a2 + b2                 a2 + b2
>    {x -> ------------------------------------------------, 
                                  a
 
            -bP - Sqrt[-(a4S) - a2b2S - a2P2 ]
>     y -> ----------------------------------------}}
                           a2 + b2

lon = ArcTan(y/x)

There are at most two solutions for each edge of the swath.  Of the four possible solutions we need the closest two, one east, one west.  The other two are where the planes of the swath edges cross the latitude on the other side of the sphere.
 

6) Find crossings for the extreme orbits:

Find the distance, in longitude, from the original point to the edge points (adjusted for the dateline).  Let:

D_E = lon_east - lon_p
D_W = lon_p - lon_west

Find the point D_E west of (lat_p, lon _p), which is (lat_p, lon_p - D_E), and the point D_W east of (lat_p, lon_p ), which is (lat_p, lon_p + D_W)

Note: An error is introduced into the algorithm here.  The sensor will see the points (lat_p, lon_p - D_E) and (lat_p, lon_p + D_W) some time prior to, or after, it passes over (lat_p, lon _p).  To be completely accurate we should adjust for the rotation of the Earth during the intervening time.  But the error is always small, makes less difference as it gets bigger, and adjusting for it means taking the sensor geometry into account. So we tolerate it in an effort to keep the algorithm simple.  Still, we could adjust at this point using some "generic" sensor geometry to increase the accuracy somewhat without sacrificing the generality of the algorithm.

Find the nodal crossing of the orbit (on a rotating sphere) that passes through (lat_p, lon_p - D_E) using step 1 above and call it lon _nadirWest.   Any orbit with a nodal crossing at lon_nadirWest includes the given point (lat _p , lon_p) on its east edge.

Find the nodal crossing of the orbit (on a rotating sphere) that passes through (lat_p, lon_p + D_W) using step 1 above andcall it lon _nadirEast.  Any orbit with a nodal crossing at lon_nadirEast includes the given point (lat_p, lon _p) on its west edge.

Any orbit with a nodal crossing in the range [lon_nadirWest, lon_nadirEast] must include the given point (lat_p , lon _p)

7) Merge the ranges found in steps 1-6

Each sample point created in step 0 produces a crossing range and those ranges should overlap.  Consequently they can be merged into a single crossing range that defines where the satellite must have crossed the equator on the previous ascending pass if the sensor were to see the area of interest during.

8) Query the database or other information storage system.

The crossing range can be used to create a spatial clause that searches on crossings in a database query.  For our database that clause looks like this:

where nodal_crossing between (lon_nadirWest, lon_nadirEast )

Combined with other clauses (temporal, parameter, quality flags, etc.) a single, relatively small, query is all that is required to return record locators matching the users criteria.
 

Descending Passes:

The angle theta can be found the same way - but the the point of interest is now on the opposite side of the sphere as the ascending crossing so adjustments have to be made.

sin(lon_n - (lon_p + 180)) / sin(theta) = sin(180 - lat_p ) / sin(inclination)  so
sin(lon_n - lon_p - 180) = sin(theta) *sin(lat_p ) / sin(180 - lat_inf)  so
lon_n = asin(sin(theta) *sin(lat_p) / sin(lat _inf ))  + lon_p + 180

Similarly we can find the length of the arc between (0, lon_n ) and (lat_p, lon_p)

sin(L_arc) / sin(90) = sin(180 - lat_p) / sin(180 - lat_inf )  so
L_arc = asin(sin(lat_p) / sin(lat_inf ))
 

Correcting the Longitudinal Half Swath Widths:

The distance between two points in Cartesian space is sqrt((x₁ -x₂)² + (y₁-y₂)² + (z₁-z₂)²) so the distance between left edge point and the west edge point is:
D = sqrt((x_left-x_west)² + (y_left -y_west)² + (z_left-z_west) ²)

Those two points, and the center of the small circle, create an icosceles triangle in the plane of the small circle.  We can use ordinary Euclidean trigonometry to find the vertex angle of that triangle and hence the length of the arc on the small circle.

cos(v) = (2*radius_sc² - D²)/ 2*radius _sc² = 1 - (D²/2*radius_sc² )
v = acos(1 - (D²/2*radius_sc²) )

The correction would be the distance the earth rotated in the time it took the satellite to go v degrees.

The correction for the right/east edge can be found in a similar fashion.  One must keep in mind that Earth rotates west to east, so the sensor sees the west edge point sometime after the satellite passes over the point of interest, but it sees the east edge point sometime before the satellite passes over the point of interest.
 

Special Cases:

2) The absolute value of the latitude of the point of interest is greater than the total coverage latitude.
If the point of interest is above (or below) the total coverage latitude the sensor sees that point during every orbit and the algorithm returns a crossing range of [0, 360].  Once that happens there is no need to test further sample points and really no need for a crossing range clause in the search, every orbit is a match.

3) The absolute value of the latitude of the point of interest is greater than the minimum inflection latitude.
If the point of interest is above (or below) the minmum inflection latitude it means one of the swath edges is not present at that latitude making it impossible to find the longitudinal half swath width.  Since the algorithm distinguishes between ascending and descending passes the cutoff in either direction is the longitude of the inflection point.

4) The absolute value of the latitude of the point of interest is greater than the inflection latitude.
If the point of interest is above (or below) the inflection latitude, and didn't qualify for cases 1 or 2 above, it is never the case that teh satellite passes over that point, but there is still a limited range of crossings for which the sensor will see the point.  The algorithm can still use the point itself to find the longitudinal half swath widths (with appropriate cutoffs as noted in case 3) but needs to substitute the inflection latitude for the latitude of the points in order to find the initial cross track edges points and the crossings.
 

Variations:

The algorithm as described above assumes a downward looking sensor.  As mentioned above the algorithm could be adjusted to work with forward and backward looking sensors by creating a virtual sattelite.  Technically the algorithm doesn't care where the sattelite is, it cares where the data is. If a sensor is pointed x degrees forward of the sattelite location that is functionally equivalent to a downward looking sensor on a virtual sattelite x degrees ahead.  Because the Earth rotates during the time it takes the satellite to travel x degrees the inclination of the virtual satellite may be slightly different  from the inclination of the real satellite so that adjustment to the input parameters needs to be made.

One unfortunate side effect is that orbits are often defined by where the satellite is so the orbital data from a sensor looking x degrees forward may start and end at x degrees north latitude instead of at the equator. One can adjust for this by recording the start/end circular latitude (discussed below) in the inventory and incorporating that into the search.  Or one could avoid storing those extra fields in the inventory for each and every orbit by adjusting the algorithm to use x degrees north as the reference latitude.

Side viewing sensors:

It was also mentioned above that the algorithm could easily be adjusted to work for side viewing sensors.  As described the algorithm assumes the ground track of the satellite (nadir) is in the center of the swath, so the distance to the left and right edges of the sensors field of view is just half the swath width.  To accomadate side viewing sensors we just have to be more explicit about those distances.  Instead of just the swath width as input we would need the two distances from nadir to the swath edge.

For example a sensor with a 1400 kilometer wide swath centered 200 kilometers left of nadir would have left and right distances of 900 kilometers and 500 kilometers respectively.  The same sensor at a more extreme angle, centered 800 kilometers left of nadir, would have left and right distances of 1500 kilometers and -100 kilometers respectively where the negative value of the right distance indicates the right edge of the swath is actually left of nadir. The sign is a matter of convention and one can adopt whatever convention is convenient.

Partial Orbits:

The algorithm as described assumes full orbits that start and end at the same reference latitude, but sometimes it is more practical to slice the data into partial orbits which means the start and end of the data become important.  While the orbit the data is part of may intersect the area being searched the data itself may not.  Unfortunately just the geographic location isn't sufficient, the direction of the satellite is also significant.  To combine the start/end latitude of the data with the direction of the satellite into a single number we have adopted a convention of circular latitudes.  Starting at the equator the circular latitudes are [0, 90] ascending, [90, 270] descending, and [270, 360] ascending.

Care must be taken to ensure the circular latitudes are in sync with the indexed crossing of the data.  For example an descending half orbit that starts at it's northernmost point and ends at it's southernmost point, for a satellite with an inclination of  98.78 degrees, indexed to the previous ascending crossing, would have start/end cicular latitudes of {98.78, 261.22}.  The same half orbit indexed to the next ascending crossing would have start/end circular latitudes of {-261.22, -98.78}. And an ascending half orbit indexed to the included ascending crossing would have start/end cicular latitudes of {-81.22, 81.22}.

Multiple Orbits:

It is also sometimes convenient to include more than one orbit in a data granule and circular latitudes can also be used to to compensate. For example a two orbit granule, that starts/ends on the equator during the ascending pass, indexed to inital asending crossing would have start/end circular latitudes of start/end circular latitudes of {0, 720}.  Indexed to the included ascending crossing the start/end circular latitudes would be{-360, 360}.  And indexed to the ending ascending crossing the start/end circular latitudes would be{-720, 0}

Obviously the algorithm has to compensate accordingly by computing the start/end circular latitudes of the search area for the different possibilities and creating multiple clauses for the search.  To keep the algorithm general some indication of how many orbits are in each granule, and which crossing the granules are indexed to, should be added as an input parameter.  One possible convention is to always index to the first included ascending crossing, or the previous ascending crossing if none are include (as with descending half orbits), and indicate the number of orbits in each granules with a single floating point number as an input parameter.

Lookup tables:

One advantage of the algorithm is it figures out everything mathematically rather than using lookup tables. This makes it faster, more flexible, more accurate, and less expensive.  It's faster because every database query takes a certain amount of time just to initialize which is often orders of magnitude larger than the amount of time it takes to compute the answer.  It's more flexible because changes in the orbit of a satellite can be accomadated by adjusting a few input parameters rather than recreating an entire lookup table.  This is especially true during the pre-launch timeframe when people are preparing the system.  Efforts expended creating lookup tabels prior to launch may be wasted if the sattelite doesn't achieve the expected orbit.

The backtrack algorithm is also demonstrably more flexible in adding another sensor to the system.  Instead of requiring a new lookup table for each new sensor, which are costly to create and, depending on the desire accuracy, can be quite large, the algorithm requires only the addition of a few input parameters in a single row of an existing table.  Moreover the accuracy of lookup table methods is directly related to the size of the table.  Depending on context doubling the accuracy of a lookup table can mean doubling, or quadrupling, its size, which not only increases the required storage but also increases the time it takes to query the table.  The backtrack algorithm, on the other hand, starts out as accurate as the input parameters and the computations allow.

Finally it would appear that one has to do the math in order to create the lookup tables in the first place. So actually going ahead with the creation and use of the lookup tables is extra effort.  Often that effort is deemed necessary because the orbital mechanics and spherical trigonometry involved in doing the math is quite complex, prone to error during the coding process, and computationally expensive.  But the backtrack algorithm simplifies the orbital mechanics by limiting itself to circular orbits and simplifies the math involved by using a hybred of spherical trigonometry, Cartesian solid Geometry, Euclidean planar Geometry, and simple Algebra.

Still, there are places in the algorithm where lookup tables could be used.  The longitudinal half swath widths, for example, are constant for a given latitude and direction.  The correction due to the rotation of the Earth of the longitudinal half swath widths is also constant for a given latitude and direction.  If the desired accuracy is such that the lookup table would be small a search for corrected widths based on latitude and direction might be faster than doing the math.

Alternatively the entire algorithm, or pieces of it, could be used to generate lookup tables for legacy systems that require them.  This has in fact been done for one sensor which reduced the cost of generating the table by about 90%.  And with a table generator written the costs associated with changing that table, or generating another for a different sensor, are virtually nothing.