Projecting the surface

Since we have extracted the numerical surface in spherical coordinates, we can also find the points that correspond to the tangent direction and project them into the FOV of the instrument. Let’s suppose that for a particular image, the satellite is located at the position \(\vec{r}_{sat}=(x_{sat}, y_{sat}, z_{sat})\) in the GSE system. The tangent points of the surface for this location will satisfy the condition:

(7)\[ \vec{n} \cdot (\vec{r}_{point}-\vec{r}_{sat} )= 0 \]

where \(\vec{r}_{point}\) is the position of each point of the numerical surface and \(\vec{n}\) is the normal vector of the surface for each point, in the GSE system. We can get \(\vec{r}_{point}\) by simply applying the transformation from the spherical system to the cartesian as defined in Coordinate system.

Fig. 14 Tangent points in GSE.

To get the normal vectors, we calculate the gradient of the surface with respect to \(\theta\) and \(\phi\) and normalize the corresponding vectors. We can calculate the normals for each point through their cross product:

\[ \vec{n} = \frac{\partial \vec{r}}{\partial \theta} \times \frac{\partial \vec{r}}{\partial \phi} \]

To get the tangent curve, we simply solve Equation (7) numerically, to get the points that satisfy it. In our case, we applied a simple mask with an arbitrary threshold of \(0.05^\circ\) tolerance, as a first approximate and fast solution.

This is defined in the function tangent_points of the TangentFitting.py script, where we find the normal vector for each point of the grid of the surface, calculate the vectors from the satellite to the grid points, and return the ones that have a dot product lower than the desired tolerance.

def tangent_points(R_numerical, Theta, Phi, satellite_pos):
    x = R_numerical * np.cos(Theta)
    y = R_numerical * np.sin(Theta) * np.cos(Phi)
    z = R_numerical * np.sin(Theta) * np.sin(Phi)

    dx_dtheta = np.gradient(x, Theta[:,0], axis=0)
    dy_dtheta = np.gradient(y, Theta[:,0], axis=0)
    dz_dtheta = np.gradient(z, Theta[:,0], axis=0)

    dx_dphi = np.gradient(x, Phi[0], axis=1)
    dy_dphi = np.gradient(y, Phi[0], axis=1)
    dz_dphi = np.gradient(z, Phi[0], axis=1)

    dtheta = np.stack((dx_dtheta,dy_dtheta,dz_dtheta), axis=-1)
    dphi = np.stack((dx_dphi, dy_dphi, dz_dphi), axis=-1)

    normals = np.cross(dtheta, dphi)
    norms = np.linalg.norm(normals, axis=-1, keepdims=True)
    normals /= norms    

    x_sat = x - satellite_pos[0]
    y_sat = y - satellite_pos[1]
    z_sat = z - satellite_pos[2]
    vec_to_sat = np.stack((x_sat, y_sat, z_sat), axis=-1)  # shape: (n_theta, n_phi, 3)
    norm_vec = np.linalg.norm(vec_to_sat, axis=-1, keepdims=True)
    vec_to_sat_unit = vec_to_sat / norm_vec  # shape: (n_theta, n_phi, 3)

    dot_product_num = np.sum(normals * vec_to_sat_unit, axis=-1)  # shape: (n_theta, n_phi)
    # dot_product_num = normals[..., 0] * x_sat + normals[..., 1] * y_sat + normals[..., 2] * z_sat

    tolerance = np.deg2rad(0.05)  # e.g. 5 degrees grazing
    mask = np.abs(dot_product_num) < np.cos(np.pi/2 - tolerance)
    # mask = np.abs(dot_product_num) < 0.1
    x_sol = x[mask]
    y_sol = y[mask]
    z_sol = z[mask]

    return x_sol, y_sol, z_sol

Fig. 15 Extraction of tangent point to satellite for a numerical surface.

The orientation of the satellite is described by three vectors, giving each base vector of the satellite coordinate system, expressed in the \ac{GSE} frame. To move from one system to another, it is sufficient to define the rotation matrix through these vectors and apply a translation according to the position of the satellite, meaning the origin of the satellite’s coordinate system. We construct the rotation matrix:

\[\begin{split} R = \begin{bmatrix} | & | & | \\ \mathbf{e}_x^{sat} & \mathbf{e}_y^{sat} & \mathbf{e}_z^{sat} \\ | & | & | \end{bmatrix} \end{split}\]

The relation between a vector in the GSE frame and the same vector in the satellite frame is:

\[ \vec{r}_{GSE} = R \cdot \vec{r}_{SF} + \vec{t} \]

where \(r_{SF}\) is the vector of a point in the satellite frame and \(\vec{t}\) is the translation vector. Therefore, the transformation of a point in GSE to the satellite frame is performed by the following operation:

\[ \vec{v}_{sat} = R^\top \cdot (\vec{v}_{GSE} - \vec{t}) \]

def GSE_to_Sat(x_GSE, y_GSE, z_GSE, satellite_pos, vx, vy, vz):
    """
    Convert GSE coordinates to satellite coordinates
    """
    R = np.array([vx, vy, vz])
    GSE_vector = np.array([x_GSE, y_GSE, z_GSE]) - satellite_pos
    sat_vector = np.dot(R, GSE_vector.T )  # Perform matrix multiplication
    return sat_vector[0], sat_vector[1], sat_vector[2]

In order to project this in the SXI coordinates we need to compute the angle of each point from the z axis (boresight) along the x and y direction of the satellite frame:

\[ azimuth = \arctan{\frac{x}{z}} \]

\[ elevation = \arctan{\frac{y}{z}} \]

And limit this to the FOV of the imager, with \(azimuth \in [-7.8,7.7]^\circ\) and \(elevation\in[-13.2,13.2]^\circ\).

def Sat_to_SXI(x_sat, y_sat, z_sat):
    """
    Convert satellite coordinates to SXI FOV angles.
    """
    azim = np.arctan2(x_sat, z_sat)  # X-Z plane: up/down
    elev = np.arctan2(-y_sat, z_sat)  # Y-Z plane: left/right

    # FOV limits
    azim_mask = (azim >= np.radians(-7.8)) & (azim<= np.radians(7.7))
    elev_mask = (elev >= np.radians(-13.2)) & (elev <= np.radians(13.2))
    # Combine masks: keep values inside either azim or elev FOV
    fov_mask = azim_mask & elev_mask  # <- this means inside both FOV limits
    # Apply mask
    azim = azim[fov_mask]
    elev = elev[fov_mask]

    return azim, elev

We can also define a coordinate transformation from the SXI degree coordinates to the indexes of the image, by mapping these angles to grid indices \(i_{\text{az}}\) and \(i_{\text{el}}\):

\[ i_{\text{az}} = \text{round} \left( \frac{\theta_{\text{az}} - \text{Az}{\min}}{\text{Az}{\max} - \text{Az}{\min}} \cdot (N{\text{az}} - 1) \right) \]

\[ i_{\text{el}} = \text{round} \left( \frac{\theta_{\text{el}} - \text{El}{\min}}{\text{El}{\max} - \text{El}{\min}} \cdot (N{\text{el}} - 1) \right) \]

where \(N_{\text{az}}\) and \(N_{\text{el}}\) are the number of azimuth and elevation grid points, respectively. To ensure the indices stay within valid bounds we also clip to the minimum of the grid:

\[i_{\text{az}} = \min\left(\max\left(i_{\text{az}}, 0\right), N_{\text{az}} - 1\right)\]

\[i_{\text{el}} = \min\left(\max\left(i_{\text{el}}, 0\right), N_{\text{el}} - 1\right)\]

The complete algorithm for the projection of the surface to the FOV is encapsulated in the funtion surface_tangent_SXI.

def surface_tangent_SXI(R_numerical,Theta,Phi, Az, El, satellite_pos, vx, vy, vz):
    # Solve the tangent equation in GSE
    x_tangent_GSE, y_tangent_GSE, z_tangent_GSE = tangent_points(R_numerical, Theta, Phi, satellite_pos)
    # Convert GSE coordinates to satellite coordinates
    x_tangent_sat , y_tangent_sat, z_tangent_sat = np.zeros((len(x_tangent_GSE))), np.zeros((len(x_tangent_GSE))), np.zeros((len(x_tangent_GSE)))
    for i in range(len(x_tangent_GSE)):
        x_tangent_sat[i], y_tangent_sat[i], z_tangent_sat[i] = GSE_to_Sat(x_tangent_GSE[i], y_tangent_GSE[i], z_tangent_GSE[i], satellite_pos, vx, vy, vz)
    # Convert satellite coordinates to SXI FOV angles
    azim_SXI, elev_SXI = Sat_to_SXI(x_tangent_sat, y_tangent_sat, z_tangent_sat)
    # interpolate the azimuth and elevation angles to the image grid
    az_ind, el_ind = interpolate(azim_SXI, elev_SXI, Az, El)
    return azim_SXI, elev_SXI, z_tangent_sat