Sensor Modeling

This tutorial demonstrates how to model sensor geometries.

Setup

import numpy as np
import pandas as pd

import plotly.graph_objs as go

from ostk.mathematics.geometry import Angle as MathAngle
from ostk.mathematics.geometry.d2.object import Point as Point2d
from ostk.mathematics.geometry.d2.object import LineString as LineString2d
from ostk.mathematics.geometry.d2.object import Polygon as Polygon2d
from ostk.mathematics.geometry.d3.object import Point as Point3d
from ostk.mathematics.geometry.d3.object import Polygon as Polygon3d
from ostk.mathematics.geometry.d3.object import Pyramid
from ostk.mathematics.geometry.d3.object import Cone
from ostk.mathematics.geometry.d3.transformation.rotation import Quaternion
from ostk.mathematics.geometry.d3.transformation.rotation import RotationVector

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.coordinate.spherical import LLA
from ostk.physics.coordinate.spherical import AER
from ostk.physics.coordinate import Position
from ostk.physics.coordinate import Frame
from ostk.physics.coordinate import Transform
from ostk.physics.coordinate.frame.provider import Dynamic as DynamicFrameProvider
from ostk.physics import Environment
from ostk.physics.environment.object import Geometry
from ostk.physics.environment.object import Celestial

from ostk.astrodynamics.trajectory import Orbit
from ostk.astrodynamics.trajectory.orbit.model import Kepler
from ostk.astrodynamics.trajectory.orbit.model.kepler import COE

---

Computation

Environment

environment = Environment.default()

earth = environment.access_celestial_object_with_name("Earth")

earth_geometry = earth.get_geometry_in(Frame.ITRF(), environment.get_instant())

Satellite

Let's define the orbit of a satellite in LEO. In this example, we're modeling the orbit using SGP4.

a = Length.kilometers(7000.0)
e = 0.0
i = Angle.degrees(45.0)
raan = Angle.degrees(0.0)
aop = Angle.degrees(0.0)
nu = Angle.degrees(0.0)

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

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

keplerian_model = Kepler(
    coe, epoch, earth, Kepler.PerturbationType.No, True
)  # True = COE expressed in ITRF frame

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[6], line 3
      1 epoch = Instant.date_time(DateTime(2018, 9, 5, 0, 0, 0), Scale.UTC)
----> 3 keplerian_model = Kepler(
      4     coe, epoch, earth, Kepler.PerturbationType.No, True
      5 )  # True = COE expressed in ITRF frame

RuntimeError:  0# ostk::core::error::RuntimeError::RuntimeError<>(char const*) at /usr/local/include/OpenSpaceToolkit/Core/Error/RuntimeError.hpp:33
 1# ostk::astrodynamics::trajectory::orbit::model::kepler::COE::getCartesianState(ostk::physics::unit::Derived const&, std::shared_ptr<ostk::physics::coordinate::Frame const> const&) const at /app/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Kepler/COE.cpp:331
 2# ostk::astrodynamics::trajectory::orbit::model::Kepler::InertialCoeFromFixedCoe(ostk::astrodynamics::trajectory::orbit::model::kepler::COE const&, ostk::physics::time::Instant const&, ostk::physics::environment::object::Celestial const&) at /app/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Kepler.cpp:318
 3# ostk::astrodynamics::trajectory::orbit::model::Kepler::Kepler(ostk::astrodynamics::trajectory::orbit::model::kepler::COE const&, ostk::physics::time::Instant const&, ostk::physics::environment::object::Celestial const&, ostk::astrodynamics::trajectory::orbit::model::Kepler::PerturbationType const&, bool) at /app/src/OpenSpaceToolkit/Astrodynamics/Trajectory/Orbit/Model/Kepler.cpp:57
 4# 0x00007F94DEE6657E in /usr/local/lib/python3.11/dist-packages/ostk/astrodynamics/OpenSpaceToolkitAstrodynamicsPy.cpython-311-x86_64-linux-gnu.so
 5# 0x00007F94DED8D247 in /usr/local/lib/python3.11/dist-packages/ostk/astrodynamics/OpenSpaceToolkitAstrodynamicsPy.cpython-311-x86_64-linux-gnu.so
 6# 0x00007F94DF06E68A in /usr/local/lib/python3.11/dist-packages/ostk/astrodynamics/OpenSpaceToolkitAstrodynamicsPy.cpython-311-x86_64-linux-gnu.so
 7# 0x00007F94DF032EF6 in /usr/local/lib/python3.11/dist-packages/ostk/astrodynamics/OpenSpaceToolkitAstrodynamicsPy.cpython-311-x86_64-linux-gnu.so
 8# 0x00007F94DEFC14B7 in /usr/local/lib/python3.11/dist-packages/ostk/astrodynamics/OpenSpaceToolkitAstrodynamicsPy.cpython-311-x86_64-linux-gnu.so
 9# 0x00007F94DEFC16D9 in /usr/local/lib/python3.11/dist-packages/ostk/astrodynamics/OpenSpaceToolkitAstrodynamicsPy.cpython-311-x86_64-linux-gnu.so
10# 0x00007F94DEC87C28 in /usr/local/lib/python3.11/dist-packages/ostk/astrodynamics/OpenSpaceToolkitAstrodynamicsPy.cpython-311-x86_64-linux-gnu.so
11# 0x00000000005EDDA9 in /usr/bin/python3.11
12# _PyObject_MakeTpCall in /usr/bin/python3.11
13# 0x0000000000639FFA in /usr/bin/python3.11
14# _PyObject_Call in /usr/bin/python3.11
15# 0x00000000005C819B in /usr/bin/python3.11
16# 0x00000000005D6834 in /usr/bin/python3.11
17# 0x00007F94E67C1B7D in /usr/local/lib/python3.11/dist-packages/ostk/core/OpenSpaceToolkitCorePy.cpython-311-x86_64-linux-gnu.so
18# _PyObject_MakeTpCall in /usr/bin/python3.11
19# _PyEval_EvalFrameDefault in /usr/bin/python3.11
20# 0x00000000006E142F in /usr/bin/python3.11
21# PyEval_EvalCode in /usr/bin/python3.11
22# 0x00000000006E32F3 in /usr/bin/python3.11
23# _PyEval_EvalFrameDefault in /usr/bin/python3.11
24# 0x000000000070F37C in /usr/bin/python3.11
25# PyIter_Send in /usr/bin/python3.11
26# _PyEval_EvalFrameDefault in /usr/bin/python3.11
27# 0x000000000070F37C in /usr/bin/python3.11
28# PyIter_Send in /usr/bin/python3.11
29# _PyEval_EvalFrameDefault in /usr/bin/python3.11
30# 0x000000000070F37C in /usr/bin/python3.11
31# 0x000000000070F617 in /usr/bin/python3.11
32# 0x000000000062AB87 in /usr/bin/python3.11
33# PyObject_Vectorcall in /usr/bin/python3.11
34# _PyEval_EvalFrameDefault in /usr/bin/python3.11
35# 0x000000000056567E in /usr/bin/python3.11
36# 0x0000000000639869 in /usr/bin/python3.11
37# PyObject_Call in /usr/bin/python3.11
38# _PyEval_EvalFrameDefault in /usr/bin/python3.11
39# 0x000000000070F37C in /usr/bin/python3.11
40# PyIter_Send in /usr/bin/python3.11
41# _PyEval_EvalFrameDefault in /usr/bin/python3.11
42# 0x000000000070F37C in /usr/bin/python3.11
43# PyIter_Send in /usr/bin/python3.11
44# _PyEval_EvalFrameDefault in /usr/bin/python3.11
45# 0x000000000070F37C in /usr/bin/python3.11
46# PyIter_Send in /usr/bin/python3.11
47# _PyEval_EvalFrameDefault in /usr/bin/python3.11
48# 0x000000000070F37C in /usr/bin/python3.11
49# PyIter_Send in /usr/bin/python3.11
50# _PyEval_EvalFrameDefault in /usr/bin/python3.11
51# 0x000000000070F37C in /usr/bin/python3.11
52# PyIter_Send in /usr/bin/python3.11
53# _PyEval_EvalFrameDefault in /usr/bin/python3.11
54# 0x000000000070F37C in /usr/bin/python3.11
55# PyIter_Send in /usr/bin/python3.11
56# 0x00007F952332EB7C in /usr/lib/python3.11/lib-dynload/_asyncio.cpython-311-x86_64-linux-gnu.so
57# 0x00007F952332F3A7 in /usr/lib/python3.11/lib-dynload/_asyncio.cpython-311-x86_64-linux-gnu.so
58# 0x00000000005EDF12 in /usr/bin/python3.11
59# 0x000000000049BB80 in /usr/bin/python3.11
60# 0x000000000049BBF8 in /usr/bin/python3.11
61# 0x00000000005EDE6B in /usr/bin/python3.11
62# PyObject_Call in /usr/bin/python3.11
63# _PyEval_EvalFrameDefault in /usr/bin/python3.11
64# 0x00000000006E142F in /usr/bin/python3.11
65# PyEval_EvalCode in /usr/bin/python3.11
66# 0x00000000006E32F3 in /usr/bin/python3.11
67# 0x00000000005EDE6B in /usr/bin/python3.11
68# PyObject_Vectorcall in /usr/bin/python3.11
69# _PyEval_EvalFrameDefault in /usr/bin/python3.11
70# _PyFunction_Vectorcall in /usr/bin/python3.11
71# PyObject_Call in /usr/bin/python3.11
72# 0x00000000006B3270 in /usr/bin/python3.11
73# Py_RunMain in /usr/bin/python3.11
74# Py_BytesMain in /usr/bin/python3.11
75# __libc_start_main in /lib/x86_64-linux-gnu/libc.so.6
76# _start in /usr/bin/python3.11
Frame must be Quasi Inertial

keplerian_model.get_classical_orbital_elements()

-- Classical Orbital Elements ----------------------------------------------------------------------
    Semi-major axis:                         6999999.9999999972 [m]                   
    Eccentricity:                            2.6250607414312717e-16                   
    Inclination:                             44.972897836604517 [deg]                 
    Right ascension of the ascending node:   343.91390293099067 [deg]                 
    Argument of periapsis:                   0.0 [deg]                                
    True anomaly:                            359.86030748978823 [deg]                 
----------------------------------------------------------------------------------------------------

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

satellite_orbit = Orbit(keplerian_model, earth)

start_instant = Instant.date_time(DateTime(2018, 9, 5, 0, 0, 0), Scale.UTC)
end_instant = start_instant + coe.get_orbital_period(
    earth.get_gravitational_parameter()
)

interval = Interval.closed(start_instant, end_instant);

step = Duration.seconds(300.0)

instants = interval.generate_grid(step)

states = [satellite_orbit.get_state_at(instant) for instant in instants]

states_lla = [
    LLA.cartesian(
        state.in_frame(Frame.ITRF()).get_position().get_coordinates(),
        earth.get_equatorial_radius(),
        earth.get_flattening(),
    )
    for state in states
]

states_line_string = LineString2d(
    [
        Point2d(
            state_lla.get_longitude().in_degrees(),
            state_lla.get_latitude().in_degrees(),
        )
        for state_lla in states_lla
    ]
)

ground_track_df = pd.DataFrame(
    [
        [
            float(state_lla.get_longitude().in_degrees()),
            float(state_lla.get_latitude().in_degrees()),
        ]
        for state_lla in states_lla
    ],
    columns=["Longitude", "Latitude"],
);

ground_track_df.head()

	Longitude	Latitude
0	3.356969e-12	3.376763e-12
1	1.207971e+01	1.306305e+01
2	2.559488e+01	2.535757e+01
3	4.214941e+01	3.585366e+01
4	6.306843e+01	4.302728e+01

Target

latitude = Angle.degrees(0.0)
longitude = Angle.degrees(0.0)
altitude = Length.meters(10.0)

target_lla = LLA(latitude, longitude, altitude)

target_position = Position.meters(
    target_lla.to_cartesian(earth.get_equatorial_radius(), earth.get_flattening()),
    Frame.ITRF(),
)

Sensor

orbital_frame = satellite_orbit.get_orbital_frame(Orbit.FrameType.NED)
orbital_frame = satellite_orbit.get_orbital_frame(Orbit.FrameType.VVLH)

alt_states = [
    orbital_frame.get_origin_in(Frame.GCRF(), instant) for instant in instants
]

alt_states_lla = [
    LLA.cartesian(
        orbital_frame.get_origin_in(Frame.ITRF(), instant).get_coordinates(),
        earth.get_equatorial_radius(),
        earth.get_flattening(),
    )
    for instant in instants
]

alt_ground_track_df = pd.DataFrame(
    [
        [
            float(state_lla.get_longitude().in_degrees()),
            float(state_lla.get_latitude().in_degrees()),
        ]
        for state_lla in alt_states_lla
    ],
    columns=["Longitude", "Latitude"],
);

ground_track_df = alt_ground_track_df

axes = [
    np.transpose(orbital_frame.get_axes_in(Frame.GCRF(), instant).x())
    for instant in instants
]

def calculate_attitude(state, target):
    q_ORB_GCRF = (
        Frame.GCRF()
        .get_transform_to(orbital_frame, state.get_instant())
        .get_orientation()
    )

    q_B_ORB = Quaternion.unit()
    q_B_ORB = Quaternion.rotation_vector(RotationVector.x(MathAngle.degrees(0.0)))

    q_B_GCRF = (q_B_ORB * q_ORB_GCRF).to_normalized()

    return q_B_GCRF

def body_frame_transform_generator(instant):
    state = satellite_orbit.get_state_at(instant)

    q_B_GCRF = calculate_attitude(state, target_position)

    return Transform.passive(
        instant,
        -state.get_position().get_coordinates(),
        np.array((0.0, 0.0, 0.0)),
        q_B_GCRF,
        np.array((0.0, 0.0, 0.0)),
    )


body_frame_provider = DynamicFrameProvider(body_frame_transform_generator)

if Frame.exists("Body"):
    Frame.destruct("Body")

body_frame = Frame.construct("Body", False, Frame.GCRF(), body_frame_provider)

def calculate_intersection(target, state, sensor_geometry):
    target_lla = LLA.cartesian(
        target.get_coordinates(), earth.get_equatorial_radius(), earth.get_flattening()
    )

    ned_frame = earth.get_frame_at(target_lla, Celestial.FrameType.NED)

    target_position_NED = target.in_frame(ned_frame, state.get_instant())
    satellite_position_NED = state.get_position().in_frame(
        ned_frame, state.get_instant()
    )

    aer = AER.from_position_to_position(target_position_NED, satellite_position_NED)

    observer_geometry_ITRF = sensor_geometry.in_frame(Frame.ITRF(), state.get_instant())

    intersection_ITRF = observer_geometry_ITRF.intersection_with(earth_geometry)

    if not intersection_ITRF.is_defined():
        return None

    intersection_points = [
        Point2d(lla.get_longitude().in_degrees(), lla.get_latitude().in_degrees())
        for lla in [
            LLA.cartesian(
                point_ITRF.as_vector(),
                earth.get_equatorial_radius(),
                earth.get_flattening(),
            )
            for point_ITRF in intersection_ITRF.access_composite()
            .access_object_at(0)
            .as_line_string()
        ]
    ]

    if len(intersection_points) < 3:
        return None

    intersection_polygon = Polygon2d(intersection_points)

    return intersection_ITRF

apex_B = Point3d(0.0, 0.0, 0.0)
base_B = Polygon3d(
    Polygon2d(
        [
            Point2d(-1.0, -1.0),
            Point2d(+1.0, -1.0),
            Point2d(+1.0, +1.0),
            Point2d(-1.0, +1.0),
        ]
    ),
    apex_B + np.array((0.0, 0.0, 1.0)),
    np.array((1.0, 0.0, 0.0)),
    np.array((0.0, 1.0, 0.0)),
)

pyramid_B = Pyramid(base_B, apex_B)

cone_B = Cone(apex_B, np.array((0.0, 0.0, 1.0)), MathAngle.degrees(30.0))

sensor_geometry = Geometry(pyramid_B, body_frame)

intersections_ITRF = [
    calculate_intersection(target_position, state, sensor_geometry) for state in states
]

intersections_ITRF = [
    intersection_ITRF
    for intersection_ITRF in intersections_ITRF
    if intersection_ITRF is not None
]

intersections_LLs = [
    [
        Point2d(lla.get_longitude().in_degrees(), lla.get_latitude().in_degrees())
        for lla in [
            LLA.cartesian(
                point_ITRF.as_vector(),
                earth.get_equatorial_radius(),
                earth.get_flattening(),
            )
            for point_ITRF in intersection_ITRF.access_composite()
            .access_object_at(0)
            .as_line_string()
        ]
    ]
    for intersection_ITRF in intersections_ITRF
]

intersection_dfs = [
    pd.DataFrame(
        [
            [float(intersection_point.x()), float(intersection_point.y())]
            for intersection_point in intersection_LLs
        ],
        columns=["Longitude", "Latitude"],
    )
    for intersection_LLs in intersections_LLs
];

intersection_dfs[0].head()

	Longitude	Latitude
0	-8.915019	3.339982e-12
1	-7.484489	-1.152781e+00
2	-6.373103	-2.114364e+00
3	-5.433459	-2.975470e+00
4	-4.586204	-3.792032e+00

---

Visualization

2D plot:

data = []

data.append(
    go.Scattergeo(
        lon=ground_track_df["Longitude"],
        lat=ground_track_df["Latitude"],
        mode="lines",
        line=dict(
            width=1,
            color="rgba(255, 0, 0, 0.5)",
        ),
    )
)

for intersection_df in intersection_dfs:
    data.append(
        go.Scattergeo(
            lon=intersection_df["Longitude"],
            lat=intersection_df["Latitude"],
            mode="lines",
            line=dict(
                width=1,
                color="rgba(0, 0, 255, 0.5)",
            ),
        )
    )

figure = go.Figure(
    data=data,
    layout=go.Layout(
        title=None,
        showlegend=False,
        height=600,
        width=1200,
        geo=go.layout.Geo(
            showland=True,
            landcolor="rgb(243, 243, 243)",
            countrycolor="rgb(204, 204, 204)",
            lonaxis=dict(showgrid=True, gridcolor="rgb(102, 102, 102)", gridwidth=0.1),
            lataxis=dict(showgrid=True, gridcolor="rgb(102, 102, 102)", gridwidth=0.1),
        ),
    ),
)

figure.show("svg")

3D plot:

data = []

data.append(
    go.Scattergeo(
        lon=ground_track_df["Longitude"],
        lat=ground_track_df["Latitude"],
        mode="lines",
        line=dict(width=1, color="rgba(255, 0, 0, 0.5)"),
    )
)

for intersection_df in intersection_dfs:
    data.append(
        go.Scattergeo(
            lon=intersection_df["Longitude"],
            lat=intersection_df["Latitude"],
            mode="lines",
            line=dict(
                width=2,
                color="rgba(0, 0, 255, 0.5)",
            ),
        )
    )

figure = go.Figure(
    data=data,
    layout=go.Layout(
        title=None,
        showlegend=False,
        height=1000,
        geo=go.layout.Geo(
            showland=True,
            showlakes=True,
            showcountries=False,
            showocean=True,
            countrywidth=0.0,
            landcolor="rgb(100, 100, 100)",
            lakecolor="rgb(240, 240, 240)",
            oceancolor="rgb(240, 240, 240)",
            projection=dict(type="orthographic", rotation=dict(lon=0, lat=0, roll=0)),
            lonaxis=dict(showgrid=True, gridcolor="rgb(102, 102, 102)", gridwidth=0.5),
            lataxis=dict(showgrid=True, gridcolor="rgb(102, 102, 102)", gridwidth=0.5),
        ),
    ),
)

figure.show("svg")

Satellite body frame orientation plot:

def earth_data(earth):
    theta = np.linspace(0, 2 * np.pi, 30)
    phi = np.linspace(0, np.pi, 30)

    theta_grid, phi_grid = np.meshgrid(theta, phi)

    r = float(earth.get_equatorial_radius().in_meters())

    x = r * np.cos(theta_grid) * np.sin(phi_grid)
    y = r * np.sin(theta_grid) * np.sin(phi_grid)
    z = r * np.cos(phi_grid)

    earth_figure = go.Surface(x=x, y=y, z=z, colorscale="Viridis", showscale=False)

    return earth_figure

def orbit_data(states):
    states_GCRF = [
        state.in_frame(Frame.GCRF()).get_position().get_coordinates()
        for state in states
    ]

    trace = go.Scatter3d(
        x=[state_GCRF[0] for state_GCRF in states_GCRF],
        y=[state_GCRF[1] for state_GCRF in states_GCRF],
        z=[state_GCRF[2] for state_GCRF in states_GCRF],
        mode="lines",
        marker=dict(size=0, showscale=False),
        line=dict(width=1),
    )

    return trace

def frame_data(frame, instants, scale):
    frame_origins_GCRF = [
        np.transpose(frame.get_origin_in(Frame.GCRF(), instant).get_coordinates())
        for instant in instants
    ]
    frame_axess_GCRF = [
        frame.get_axes_in(Frame.GCRF(), instant) for instant in instants
    ]

    def axis_data(origin, axis, color):
        data = go.Scatter3d(
            x=[origin[0], origin[0] + scale * axis],
            y=[origin[1], origin[1] + scale * axis],
            z=[origin[2], origin[2] + scale * axis],
            marker=dict(size=1),
            line=dict(color=color, width=6),
        )

        return data

    data = []

    for [frame_origin_GCRF, frame_axes_GCRF] in zip(
        frame_origins_GCRF, frame_axess_GCRF
    ):
        data.append(
            axis_data(frame_origin_GCRF, np.transpose(frame_axes_GCRF.x())[0], "red")
        )
        data.append(
            axis_data(frame_origin_GCRF, np.transpose(frame_axes_GCRF.y())[0], "blue")
        )
        data.append(
            axis_data(frame_origin_GCRF, np.transpose(frame_axes_GCRF.z())[0], "green")
        )

    return data

figure = go.Figure(
    data=[
        earth_data(earth),
        orbit_data(states),
        *frame_data(body_frame, instants, scale=500000),
    ],
    layout=go.Layout(
        title=None,
        showlegend=False,
        height=1000,
        width=600,
        scene=dict(
            xaxis=dict(
                gridcolor="rgb(255, 255, 255)",
                zerolinecolor="rgb(255, 255, 255)",
                showbackground=True,
                backgroundcolor="rgb(230, 230,230)",
            ),
            yaxis=dict(
                gridcolor="rgb(255, 255, 255)",
                zerolinecolor="rgb(255, 255, 255)",
                showbackground=True,
                backgroundcolor="rgb(230, 230,230)",
            ),
            zaxis=dict(
                gridcolor="rgb(255, 255, 255)",
                zerolinecolor="rgb(255, 255, 255)",
                showbackground=True,
                backgroundcolor="rgb(230, 230,230)",
            ),
            camera=dict(
                up=dict(x=0, y=0, z=1),
                eye=dict(
                    x=-1.7428,
                    y=0.5707,
                    z=0.9100,
                ),
            ),
            aspectratio=dict(x=1, y=1, z=1),
            aspectmode="manual",
        ),
    ),
)

figure.show("png")

---