How can I perform rotation of my satellite sensor in relation to nadir?

Hey folks,

My satellite has field of regard (FOR) of 90° (-45° to +45°). Its sensor is capable of imaging a 12km width area. Also, it has nadir pointing at 0 degrees.

I am struggling a bit to construct the LOS + the rotation in relation to the nadir.

When I put LOS {0,0,1} I successfully plot the point in the same projected path of the satellite.

EDIT: the blue dot near the satellite’s path, that’s my geolocated point with LOS {001}.

However, my control image falls in continental Portugal, which means there’s a rotation of the instrument in relation to the nadir pointing.

How do I apply this rotation to my LOS so I can correctly geolocate my points?

So far, I’ve trying to use the quaternions, but the output shows me unexpected locations.

I thank you all in advance for any help.

Fields of view are defined according to the spacecraft frame. The attitude of the satellite is what is used to compute the orientation of the spacecraft frame with respect to the inertial frame.
What attitude provider did you register to the orbit propagator? You should probably set up something like a LofOffset provider configured with one of the LOFType constants, or a NadirPointing provider, eventually wrapped within a YawCompensation provider if needed.

Hey, Luc, thanks for the reply! And sorry for taking some time to get back to you. I’ve been trying to gather enough information to reply to you.

Well, from the image above, I was using the following setting of the propagator, with attitude provider as NadirPointing:

#Creating the initial orbit based on PV, Frame and Date.
initial_orbit = KeplerianOrbit(pv_gcrf, gcrf, date_i, Constants.WGS84_EARTH_MU)

#Setting the Propagator.
propagator = EcksteinHechlerPropagator(initial_orbit, #Initial Orbit
                                    NadirPointing(gcrf, earth), #Attitude Provider
                                    Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS, #
                                    Constants.EIGEN5C_EARTH_MU, 
                                    Constants.EIGEN5C_EARTH_C20,
                                    Constants.EIGEN5C_EARTH_C30, 
                                    Constants.EIGEN5C_EARTH_C40,
                                    Constants.EIGEN5C_EARTH_C50, 
                                    Constants.EIGEN5C_EARTH_C60)

But I was not being able to get the rotation.

Now, seeking simplification, I am trying to geolocate only one pixel.

attitude_provider = LofOffset(gcrf,  
                              LOFType.VNC #Constant for Velocity - Normal - Co-normal frame 
                                        #(X axis aligned with velocity, 
                                        # Y axis aligned with orbital momentum).)

#Setting the Propagator.
propagator = EcksteinHechlerPropagator(initial_orbit,
                                    attitude_provider, #NadirPointing(gcrf, earth),
                                    Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
                                    Constants.EIGEN5C_EARTH_MU, 
                                    Constants.EIGEN5C_EARTH_C20,
                                    Constants.EIGEN5C_EARTH_C30, 
                                    Constants.EIGEN5C_EARTH_C40,
                                    Constants.EIGEN5C_EARTH_C50, 
                                    Constants.EIGEN5C_EARTH_C60)

Now I wonder, to get the angle applied to the LOS 001, should I use the spacecraft attitude rotation and applyTo() in the LOS Vector?

Or would be necessary to convert to ECI frame?

Hello @rodreras,

If my understanding is right, you are struggling to create your FieldOfView with respect to the satellite attitude (in this case NadirPointing) ?

First thing to do is to understand what will be the reference frame of the spacecraft when specifying its attitude provider. When using NadirPointing, the spacecraft (0,0,1) vector will point towards the ground (as you could see on the image of your first post). Moreover, the (0,1,0) vector is the vector resulting from the cross product of the nadir vector and the orbital velocity i.e. opposite to the orbital momentum vector (this is what is performed at code level, javadoc should be clearer about that).

Knowing this, and depending on how you create your FieldOfView sensor, you should succeed in using the right orientation. You can use quaternions to build your rotation but if you want to do some testing first, i believe that the following Rotation constructor is more “natural” :
Rotation(Vector3D u1, Vector3D u2, Vector3D v1, Vector3D v2)

For example, in this case, u1 and u2 would be (0,1,0) and (0,0,1) as you would build your rotation starting from the spacecraft attitude.

Hope this helps ! Frames and rotations can be quite tricky to visualize and explains.

Cheers,
Vincent

2 Likes

Ok so i re-used another tutorial that i made to make a simple example of what can be done to visualize your field of view :slight_smile:.

It uses pandas and plotly to visualize the results.

BEWARE: change ‘PATH/TO/OREKIT/DATA’

Code

import orekit
import pandas as pd
import plotly.graph_objects as go
from java.util import List
from orekit import JArray_double
from org.hipparchus.geometry.euclidean.threed import Vector3D
from org.hipparchus.ode.nonstiff import DormandPrince853Integrator
from org.hipparchus.util import FastMath
from org.orekit.attitudes import LofOffset, NadirPointing
from org.orekit.bodies import CelestialBodyFactory, GeodeticPoint
from org.orekit.frames import FramesFactory, Transform
from org.orekit.geometry.fov import DoubleDihedraFieldOfView, FieldOfView
from org.orekit.models.earth import ReferenceEllipsoid
from org.orekit.orbits import EquinoctialOrbit
from org.orekit.orbits import OrbitType
from org.orekit.orbits import PositionAngleType
from org.orekit.propagation import SpacecraftState, Propagator
from org.orekit.propagation.numerical import NumericalPropagator
from org.orekit.propagation.sampling import PythonOrekitFixedStepHandler
from org.orekit.time import AbsoluteDate
from org.orekit.time import TimeScalesFactory
from org.orekit.utils import Constants
from org.orekit.utils import IERSConventions
from pandas import DataFrame

vm = orekit.initVM()

from orekit.pyhelpers import setup_orekit_curdir

setup_orekit_curdir('PATH/TO/OREKIT/DATA')

### CONSTANTS
UTC = TimeScalesFactory.getUTC()
INERTIAL = FramesFactory.getEME2000()
ITRF = FramesFactory.getITRF(IERSConventions.IERS_2010, False)
EARTH_BODYSHAPE = ReferenceEllipsoid.getIers2010(ITRF)
EARTH_RADIUS = Constants.WGS84_EARTH_EQUATORIAL_RADIUS
INITIAL_DATE = AbsoluteDate(2024, 8, 30, 15, 0, 00.000, UTC)


### CLASSES
class FOVHandler(PythonOrekitFixedStepHandler):

    def init(self, spacecraftState, absoluteDate, double):
        pass

    def finish(self, spacecraftState):
        pass

    def __init__(self, fov: FieldOfView, angular_step: float):
        super().__init__()
        self.fov = fov
        self.angular_step = angular_step
        self.footprint = []

    def handleStep(self, state):
        inertToBody = state.getFrame().getTransformTo(EARTH_BODYSHAPE.getBodyFrame(), state.getDate())
        fovToBody = Transform(state.getDate(), state.toTransform().getInverse(), inertToBody)
        self.footprint.append(self.fov.getFootprint(fovToBody, EARTH_BODYSHAPE, self.angular_step))

    def get_footprint(self):
        return self.footprint


### METHODS
def get_numerical_propagator(initial_orbit, attitude_provider):
    min_step = 0.001
    maxstep = 1000.0
    init_step = 60.0
    position_tolerance = 1.0

    tolerances = NumericalPropagator.tolerances(position_tolerance,
                                                initial_orbit,
                                                initial_orbit.getType())

    integrator = DormandPrince853Integrator(min_step, maxstep,
                                            JArray_double.cast_(tolerances[0]),
                                            # Double array of doubles needs to be casted in Python
                                            JArray_double.cast_(tolerances[1]))
    integrator.setInitialStepSize(init_step)

    return NumericalPropagator(integrator, attitude_provider)

def get_orbit():
    # Define equatorial orbit
    a = 2000. * 1000. + EARTH_RADIUS
    return EquinoctialOrbit(a, 0., 0., 0., 0.,
                            FastMath.toRadians(30.), PositionAngleType.TRUE, INERTIAL,
                            INITIAL_DATE, Constants.WGS84_EARTH_MU)


# Method to convert footprint to usable data (comes from another solution on the forum https://forum.orekit.org/t/access-getfootprint-data/2485)
def get_footprint_list(handler):
    footprint_list = handler.get_footprint()
    converted_footprint_list = []
    for footprint in footprint_list:
        footprint = list(footprint)
        footprint = [list(List.cast_(points)) for points in footprint]
        converted_footprint = []
        for points in footprint:
            converted_points = []
            for point in points:
                converted_points.append(GeodeticPoint.cast_(point))
            converted_footprint.append(converted_points)
        converted_footprint_list.append(converted_footprint)
    return converted_footprint_list


def plot_footprint(footprint_list):
    # Create empty figure to fill with footprints
    fig = go.Figure()

    # Fill figure
    for footprint in footprint_list:
        df = DataFrame()
        for points in footprint:
            points = [[point.getLatitude(), point.getLongitude()] for point in points]
            df = pd.concat([df, DataFrame(points, columns=['Latitude', 'Longitude'])])
            df.loc[-1] = points[0]  # adding first point as last point to close figure on plot

        fig.add_scattergeo(lat=df['Latitude'], lon=df['Longitude'], mode='lines')


    # Display figure
    fig.show()


### MAIN
if __name__ == '__main__':
    # Define nadir pointing (by default, vector (0,0,1) is nadir and vector (0,1,0) is the normal to the plane defined
    # by nadir and the orbital velocity)
    sat_attitude_provider = NadirPointing(INERTIAL, EARTH_BODYSHAPE)

    # Get initial orbit
    initial_orbit = get_orbit()

    # Get Numerical propagator
    propagator_num = get_numerical_propagator(initial_orbit, sat_attitude_provider)

    # Create initial state
    initial_state = SpacecraftState(initial_orbit)

    # Set orbit type & initial state
    propagator_num.setOrbitType(OrbitType.CARTESIAN)
    propagator_num.setInitialState(initial_state)

    # Define sensor FOV (nadir is defined using vector (0,0,1), roll around (1,0,0) and pitch around (0,1,0)
    sensor_fov = DoubleDihedraFieldOfView(Vector3D.PLUS_K,
                                          Vector3D.PLUS_I, FastMath.toRadians(45.),
                                          Vector3D.PLUS_J, FastMath.toRadians(10.),
                                          0.)

    # Create custom FOV handler
    angular_step = FastMath.toRadians(90.)
    fov_handler = FOVHandler(sensor_fov, angular_step)

    # Add handler to register FOV footprint
    casted_propagator = Propagator.cast_(propagator_num) # Cast is necessary to use method from Java interface Propagator
    casted_propagator.setStepHandler(300., fov_handler)

    # Propagate for one hour
    propagator_num.propagate(initial_orbit.getDate().shiftedBy(3600.))

    # Get footprint list
    footprint_list = get_footprint_list(fov_handler)

    # Plot on map
    plot_footprint(footprint_list)

Plots


Cheers,
Vincent

1 Like

Another (cleaner and more complete) code example is available here : FOV and grountrack with FootprintOverlapDetector

Cheers,
Vincent CUCCHIETTI

Hello, @Vincent

Thanks a lot for the explanations, examples and sources. They are going to be very useful to me.

It turns out that I was capable of applying the rotation to my satellite sensor based on the angle generated from the Quaternions. I learned that the convention for rotating with Orekit is the opposite from what one would expect.

As a result, I could geolocate my pixels. For now, it works for me, but I could not geolocate using Rugged and obtain the coordinate elevation, which makes things more precise.

I will put here the workflow for geolocation without DEM, so other people can profit from it:

1. - Get satellite position
2.- Transform vectors from ECEF to ECI
3.- Convert time to UTC AbsoluteDate
4. - Define LOS. I used los_body_nadir = Vector3D(0.,0.,1.), since my satellite points to nadir.
5. - Generate rotation from the Quaternions
6. - Rotate the LOS based on the rotation
7. - Transform LOS from ECI to ECEF, and then convert to StaticTransform seeking the transform of the LOS vector.
8. - Perform intersection with Earth by building a line

#Setting Earth ellipsoid
earth = OneAxisEllipsoid(
                        Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
                        Constants.WGS84_EARTH_FLATTENING,
                        itrf
)

#Building the line for intersection
line = Line(pos_ecef,
            pos_ecef.add(los_body_rotated_itrf),
            1e-10
) 

#Intersecting with Earth
intersection = earth.getIntersectionPoint(
                                        line, #the line
                                        pos_ecef, #sat pos
                                        itrf, #frame
                                        date #date
)

#Coordinates
latitude = math.degrees(intersection.getLatitude())
longitude = math.degrees(intersection.getLongitude())
altitude = math.degrees(intersection.getAltitude()) 

With that, I could geolocate the 1st pixel of my image, falling exactly where it should be.

Make sure that, when you apply the rotation to your vectors, always check whether you should use applyTo() or applyInverseTo().

I think that, by finding this solution, I can see that my original post was not clear enough. But still I appreciate the help from you all! It was not in vain, since I’m still going to use it. For next time, I’ll try to be more accurate.

All the best,

Rodrigo.