Access Computation¶

This tutorial demonstrates how to compute access.

Setup¶

import pandas as pd

import plotly.graph_objs as go

from ostk.mathematics.object import RealInterval

from ostk.physics import Environment
from ostk.physics.unit import Length
from ostk.physics.unit import Angle
from ostk.physics.time import Scale
from ostk.physics.time import Instant
from ostk.physics.time import Duration
from ostk.physics.time import Interval
from ostk.physics.time import DateTime
from ostk.physics.time import Time
from ostk.physics.coordinate.spherical import LLA
from ostk.astrodynamics.trajectory import Orbit
from ostk.astrodynamics.trajectory.orbit.model import Kepler
from ostk.astrodynamics.trajectory.orbit.model.kepler import COE
from ostk.astrodynamics.trajectory.orbit.model import SGP4
from ostk.astrodynamics.trajectory.orbit.model.sgp4 import TLE
from ostk.astrodynamics.access import Generator as AccessGenerator
from ostk.astrodynamics.access import AccessTarget
from ostk.astrodynamics.access import VisibilityCriterion
from ostk.astrodynamics.utilities import compute_trajectory_geometry
from ostk.astrodynamics.utilities import compute_time_lla_aer_coordinates

Access¶

An access represents an object-to-object visibility period.

In this example, let's compute accesses between a fixed position on the ground and a satellite in LEO.

Environment¶

Let's setup an environment (which describes where planets are, etc...):

environment = Environment.default(True)  # Set global environment instance as the default instance 
earth = environment.access_celestial_object_with_name("Earth")

Origin¶

Visibility Criterion¶

Let's define a visibility criterion, there are currently 4 types of Visibility Criterion

AERInterval : Defined by Azimuth, Elevation and Range intervals (typically used for accesses with Ground Stations)
ElevationInterval : Defined solely by an Elevation interval (typically used for imaging windows)
AERMask : Defined by an Azimuth-Elevation mask, and a Range interval (typically used for accesses with Ground Stations)
LineOfSight : Defined by a pure line of sight constraint, occluded by the bodies in the provided Environment (typically used for Satellite <-> Satellite accesses)

visibility_criterion = VisibilityCriterion.from_aer_interval(
    azimuth_interval=RealInterval.closed(0.0, 360.0),  # [deg]
    elevation_interval=RealInterval.closed(0.0, 90.0),  # [deg]
    range_interval=RealInterval.closed(0.0, 10000e3),  # [m]
)

# visibility_criterion = VisibilityCriterion.from_aer_mask(
#     azimuth_elevation_mask={0.0: 0.0, 90.0: 45.0, 270.0: 30.0},
#     range_interval=RealInterval.closed(0.0, 10000e3),
# )

# visibility_criterion = VisibilityCriterion.from_elevation_interval(
#     elevation_interval=RealInterval.closed(80.0, 90.0),
# )

# visibility_criterion = VisibilityCriterion.from_line_of_sight(
#     environment=environment,
# )

Access Target¶

Now that we have defined a visibility criterion, we can create an access target to compute accesses with. Let's define a fixed ground position, using its geographic coordinates. We can create a list of multiple access targets and compute a list of accesses corresponding to each access target

access_target = AccessTarget.from_lla(
    visibility_criterion=visibility_criterion,
    lla=LLA(Angle.degrees(50.0), Angle.degrees(20.0), Length.meters(30.0)),
    celestial=earth,
)
another_access_target = AccessTarget.from_lla(
    visibility_criterion=visibility_criterion,
    lla=LLA(Angle.degrees(-50.0), Angle.degrees(10.0), Length.meters(30.0)),
    celestial=earth,
)

access_targets = [access_target, another_access_target]

Target¶

Let's consider a satellite in Low-Earth Orbit.

We can define its orbit with Classical Orbital Elements:

a = earth.get_equatorial_radius() + Length.kilometers(500.0)
e = 0.000
i = Angle.degrees(97.8893)
raan = Angle.degrees(100.372)
aop = Angle.degrees(0.0)
nu = Angle.degrees(0.0201851)

coe = COE(a, e, i, raan, aop, nu)

... and by using a Keplerian orbital model:

epoch = Instant.date_time(DateTime(2018, 1, 1, 0, 0, 0), Scale.UTC)

keplerian_model = Kepler(coe, epoch, earth, Kepler.PerturbationType.J2)

Or with a Two-Line Element (TLE) set:

tle = TLE(
    "ISS (ZARYA)",
    "1 25544U 98067A   18268.86272795  .00002184  00000-0  40781-4 0  9990",
    "2 25544  51.6405 237.0010 0003980 205.4375 242.3358 15.53733046134172",
)

... along with its associated SGP4 orbital model:

sgp4_model = SGP4(tle)

Locking local repository [/var/cache/open-space-toolkit-data/data/manifest]...
Fetching Data Manifest from [https://github.com/open-space-collective/open-space-toolkit-data/raw/v1/data/manifest.json]...
Unlocking local repository [/var/cache/open-space-toolkit-data/data/manifest]...
Data Manifest [/var/cache/open-space-toolkit-data/data/manifest/manifest.json] has been successfully fetched from [https://github.com/open-space-collective/open-space-toolkit-data/raw/v1/data/manifest.json].
Fetching latest Bulletin A...
Locking local repository [/var/cache/open-space-toolkit-data/data/coordinate/frame/provider/iers]...
Creating temporary directory [/var/cache/open-space-toolkit-data/data/coordinate/frame/provider/iers/bulletin-A/tmp]...
Fetching Bulletin A from [https://github.com/open-space-collective/open-space-toolkit-data/raw/v1/data//coordinate/frame/provider/iers/bulletin-A/ser7.dat]...
Unlocking local repository [/var/cache/open-space-toolkit-data/data/coordinate/frame/provider/iers]...
Bulletin A [/var/cache/open-space-toolkit-data/data/coordinate/frame/provider/iers/bulletin-A/ser7.dat] has been successfully fetched from [https://github.com/open-space-collective/open-space-toolkit-data/raw/v1/data//coordinate/frame/provider/iers/bulletin-A/ser7.dat].

Below, we select which orbital model to use:

orbital_model = keplerian_model
# orbital_model = sgp4_model

We then obtain the satellite orbit (which is a Trajectory object):

satellite_orbit = Orbit(orbital_model, earth)

Alternatively, the Orbit class can provide some useful shortcuts (for usual orbit types):

epoch = Instant.date_time(DateTime(2018, 1, 1, 0, 0, 0), Scale.UTC)

satellite_orbit = Orbit.sun_synchronous(
    epoch, Length.kilometers(500.0), Time(12, 0, 0), earth
)

Access¶

Now that the origin and the target trajectories are well defined, we can compute the Access.

Let's first define an analysis interval:

start_instant = Instant.date_time(DateTime.parse("2018-01-01 00:00:00"), Scale.UTC)
end_instant = Instant.date_time(DateTime.parse("2018-01-10 00:00:00"), Scale.UTC)

interval = Interval.closed(start_instant, end_instant)

Then, using an Access Generator, we can compute the accesses within the intervals of interest:

access_generator = AccessGenerator(environment=environment)

accesses = access_generator.compute_accesses(
    interval=interval,
    access_targets=access_targets,
    to_trajectory=satellite_orbit,
)

assert len(accesses) == len(access_targets)  # a list of accesses per target

first_target_accesses = accesses[0]
second_target_accesses = accesses[1]
all_accesses = first_target_accesses + second_target_accesses

accesses_data = [
    (
        str(access.get_type()),
        repr(access.get_acquisition_of_signal()),
        repr(access.get_time_of_closest_approach()),
        repr(access.get_loss_of_signal()),
        float(access.get_duration().in_seconds()),
        "first_target",
    ) for access in first_target_accesses
]

accesses_data += [
    (
        str(access.get_type()),
        repr(access.get_acquisition_of_signal()),
        repr(access.get_time_of_closest_approach()),
        repr(access.get_loss_of_signal()),
        float(access.get_duration().in_seconds()),
        "second_target",
    ) for access in second_target_accesses
]

And format the output using a dataframe:

accesses_df = pd.DataFrame(
    data=accesses_data,
    columns=["Type", "AOS", "TCA", "LOS", "Duration", "Target"],
)

Output¶

Print accesses:

accesses_df

	Type	AOS	TCA	LOS	Duration	Target
0	Type.Complete	2018-01-01 00:10:48.691.727.513 [UTC]	2018-01-01 00:13:07.422.341.866 [UTC]	2018-01-01 00:15:32.871.969.978 [UTC]	284.180242	first_target
1	Type.Complete	2018-01-01 09:57:49.624.946.040 [UTC]	2018-01-01 10:02:55.818.813.623 [UTC]	2018-01-01 10:07:55.028.955.753 [UTC]	605.404010	first_target
2	Type.Complete	2018-01-01 11:31:16.376.710.963 [UTC]	2018-01-01 11:37:02.533.622.949 [UTC]	2018-01-01 11:42:43.005.386.542 [UTC]	686.628676	first_target
3	Type.Complete	2018-01-01 13:06:09.698.392.313 [UTC]	2018-01-01 13:09:58.212.688.283 [UTC]	2018-01-01 13:13:42.224.919.802 [UTC]	452.526527	first_target
4	Type.Complete	2018-01-01 20:41:15.275.981.816 [UTC]	2018-01-01 20:46:05.882.558.137 [UTC]	2018-01-01 20:51:01.381.094.552 [UTC]	586.105113	first_target
...	...	...	...	...	...	...
108	Type.Complete	2018-01-09 01:53:59.269.923.527 [UTC]	2018-01-09 01:57:28.651.787.612 [UTC]	2018-01-09 02:00:53.534.163.413 [UTC]	414.264240	second_target
109	Type.Complete	2018-01-09 09:28:38.419.034.034 [UTC]	2018-01-09 09:33:39.084.591.238 [UTC]	2018-01-09 09:38:44.723.391.949 [UTC]	606.304358	second_target
110	Type.Complete	2018-01-09 11:01:26.277.987.717 [UTC]	2018-01-09 11:07:08.135.307.687 [UTC]	2018-01-09 11:12:56.397.738.412 [UTC]	690.119751	second_target
111	Type.Complete	2018-01-09 12:38:38.602.425.206 [UTC]	2018-01-09 12:41:50.444.699.315 [UTC]	2018-01-09 12:45:09.170.303.183 [UTC]	390.567878	second_target
112	Type.Complete	2018-01-09 22:26:55.632.669.418 [UTC]	2018-01-09 22:31:43.647.618.185 [UTC]	2018-01-09 22:36:24.604.249.470 [UTC]	568.971580	second_target

113 rows × 6 columns

Let's calculate the geographic coordinate of the satellite, during access:

def compute_access_geometry(access):
    return [
        compute_time_lla_aer_coordinates(state, access_target.get_position(), environment)
        for state in satellite_orbit.get_states_at(
            access.get_interval().generate_grid(Duration.seconds(1.0))
        )
    ]

satellite_orbit_geometry_df = pd.DataFrame(
    [lla.to_vector() for lla in compute_trajectory_geometry(satellite_orbit, interval)],
    columns=["Latitude", "Longitude", "Altitude"],
)

satellite_orbit_geometry_df.head()

	Latitude	Longitude	Altitude
0	-0.020152	-0.001105	500000.002625
1	3.772321	-0.734928	500091.839402
2	7.564114	-1.472946	500367.684689
3	11.354544	-2.219521	500822.628803
4	15.142919	-2.979323	501448.578113

access_geometry_dfs = [
    pd.DataFrame(
        compute_access_geometry(access),
        columns=[
            "Time",
            "Latitude",
            "Longitude",
            "Altitude",
            "Azimuth",
            "Elevation",
            "Range",
        ],
    )
    for access in all_accesses
]

def get_max_elevation(df):
    return df.loc[df["Elevation"].idxmax()]["Elevation"]

And plot the geometries onto a map:

data = []

# Target geometry

for access_target in access_targets:
    data.append(
        dict(
            type="scattergeo",
            lon=[float(access_target.get_lla(earth).in_degrees())],
            lat=[float(access_target.get_lla(earth).in_degrees())],
            mode="markers",
            marker=dict(size=10, color="orange"),
        )
    )

# Orbit geometry

data.append(
    dict(
        type="scattergeo",
        lon=satellite_orbit_geometry_df["Longitude"],
        lat=satellite_orbit_geometry_df["Latitude"],
        mode="lines",
        line=dict(
            width=1,
            color="rgba(0, 0, 0, 0.1)",
        ),
    )
)

# Access geometry

for access_geometry_df in access_geometry_dfs:
    data.append(
        dict(
            type="scattergeo",
            lon=access_geometry_df["Longitude"],
            lat=access_geometry_df["Latitude"],
            mode="lines",
            line=dict(
                width=1,
                color="red",
            ),
        )
    )

layout = dict(
    title=None,
    showlegend=False,
    width=1200,
    height=600,
    geo=dict(
        showland=True,
        landcolor="rgb(243, 243, 243)",
        countrycolor="rgb(204, 204, 204)",
    ),
)

figure = go.Figure(data=data, layout=layout)

#figure.show("svg")
figure.show()

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[25], line 9
# Target geometry
for access_target in access_targets:
   data.append(
       dict(
           type="scattergeo",
----> 9             lon=[float(access_target.get_lla(earth).in_degrees())],
           lat=[float(access_target.get_lla(earth).in_degrees())],
           mode="markers",
           marker=dict(size=10, color="orange"),
       )
   )
# Orbit geometry
data.append(
   dict(
       type="scattergeo",
   (...)
   )
)

AttributeError: 'ostk.physics.coordinate.spherical.LLA' object has no attribute 'in_degrees'