Sentinel-1b propagation

Data¶

This notebook compares Sentinel 1b GNSS position data against propagated position data in the GCRF frame. The Sentinel 1b data is retrieved automatically from http://aux.sentinel1.eo.esa.int/POEORB/ for the hard coded time frames 2019-06-26T21:40:41.000Z to 2019-07-02T00:59:42.000Z. Maneuvers are regarded and extracted from https://sentinel.esa.int/documents/247904/3376743/s1b.man.

Comparing the GNSS data against the propagated data¶

We are starting with the results being plotted over time, showing the

• radial (R)
• along-track (T) and
• cross-track (N) errors.

In [7]:

plt.plot(range(len(propagated_gcrf_state_vectors)), radial_differences, label="R")
plt.plot(range(len(propagated_gcrf_state_vectors)), tangential_differences, label="T")
plt.plot(range(len(propagated_gcrf_state_vectors)), normal_differences, label="N")
plt.legend()
plt.xlabel('Time / 100 s')
plt.ylabel('Position difference / km')
plt.grid(True)
plt.title('Sentinel 1b GNSS / NEPTUNE Propagation Result')
plt.show()

plt.plot(range(len(propagated_gcrf_state_vectors)), radial_differences, label="R") plt.plot(range(len(propagated_gcrf_state_vectors)), tangential_differences, label="T") plt.plot(range(len(propagated_gcrf_state_vectors)), normal_differences, label="N") plt.legend() plt.xlabel('Time / 100 s') plt.ylabel('Position difference / km') plt.grid(True) plt.title('Sentinel 1b GNSS / NEPTUNE Propagation Result') plt.show()

Discussion¶

Over the 7 days interval of the downloaded validation data an errors builds up in radial, along-track and cross-track components in the propagation results. The greatest error can be found in the along-track component. After 7 days, the error is at a magnitude of 500 m. The radial and cross-track errors are 1 magnitude smaller. The reason for the increased build up of the along-track error can be found in the radial error, which in turn leads to a deviation of the velocity vector. Even a small error in the velocity vector leads to a a growing along-track error. Also the graph shows a small offset around the time a maneuver was applied, which leads to an increase of the velocity error and thus again a very visible error on the along-track component.

Below you find the code of the whole validation process, starting with the inclusion of the needed modules.

In [1]:

from okapi_pkg.okapi_init import okapi_init
from okapi_pkg.okapi_send_request import okapi_send_request
from okapi_pkg.okapi_get_result import okapi_get_result
import json
import numpy as np
from  astro_conversions import itrf_to_gcrf, get_rtn_matrix, gcrf_to_itrf
import os.path
from os import path, mkdir
import wget
import xml.dom.minidom
from tqdm import tqdm
from datetime import datetime
import matplotlib.pyplot as plt

from okapi_pkg.okapi_init import okapi_init from okapi_pkg.okapi_send_request import okapi_send_request from okapi_pkg.okapi_get_result import okapi_get_result import json import numpy as np from astro_conversions import itrf_to_gcrf, get_rtn_matrix, gcrf_to_itrf import os.path from os import path, mkdir import wget import xml.dom.minidom from tqdm import tqdm from datetime import datetime import matplotlib.pyplot as plt

Calling the NEPTUNE propagator¶

The OKAPI-Python-Connector is used to dispatch the propagation and retrieve the result. The initial state vector is extrated from the sentinel data and rotated from the ITRF to GCRF. The assumed mass is 2151.964 kg. The drag area is 6.2. An output state is created every 100 s.

In [2]:

# Make sure we only propagate once, when the data is already there, we do not repeat the request.
if not path.exists('sentinel-1b_neptune_gcrf.json'):

    okapi_login, error = okapi_init("https://api.okapiorbits.com/", "<username>", "<password>")

    if error.get('status', '') == 'FATAL':
        print(error)
        exit('Error during authentification.')
    #print(okapi_login)

    propagate_simple_request_body = {
        "orbit": {
            "type": "opm.json",
            "content": {
                "OPM_HEADER": {
                    "CCSDS_OPM_VERS": 2
                },
                "OPM_META_DATA": {
                    "OBJECT_ID": "41456",
                    "OBJECT_NAME": "Sentinel-1b",
                    "CENTER_NAME": "EARTH",
                    "REF_FRAME": "GCRF",
                    "REF_FRAME_EPOCH": "2000-01-01T00:00:00Z",
                    "TIME_SYSTEM": "UTC"
                },
                "OPM_DATA": {
                    "EPOCH": "2019-06-24T22:59:42.000Z",
                    "X":  5123.592195,
                    "Y": -462.889158,
                    "Z": -4865.970976,
                    "X_DOT": -5.16151308,
                    "Y_DOT": -1.00719510,
                    "Z_DOT": -5.34379188,
                    "MASS": 2151.964,
                    "SOLAR_RAD_COEFF": 1.3,
                    "DRAG_AREA": 6.2,
                    "DRAG_COEFF": 2.2,
                    "MANEUVERS": [
                        {
                            "MAN_EPOCH_IGNITION": "2019-06-26T21:40:41.000Z",
                            "MAN_DURATION": 23.875,
                            "MAN_REF_FRAME": "RTN",
                            "USER_DEFINED_MAN_A_1": -0.21367128E-08,
                            "USER_DEFINED_MAN_A_2": 0.34196316E-06,
                            "USER_DEFINED_MAN_A_3": 0.24668746E-07
                        },
                        {
                            "MAN_EPOCH_IGNITION": "2019-06-26T22:29:45.000Z",
                            "MAN_DURATION": 14.125,
                            "MAN_REF_FRAME": "RTN",
                            "USER_DEFINED_MAN_A_1": -0.14149311E-08,
                            "USER_DEFINED_MAN_A_2": -0.37489161E-06,
                            "USER_DEFINED_MAN_A_3": 0.26093520E-07
                        }
                    ],
                },
            },
        },
        "settings": {
            "type": "prop_settings.json",
            "content": {
                "propagation_end_epoch": "2019-07-02T00:59:42.000Z",
                "more": {
                    "output_step_size": 100,
                    "stormer_cowell_rel_tolerance": 1e-12,
                    "stormer_cowell_abs_tolerance": 1e-13,
                    "geopotential_degree_order": 70,
                    "atmospheric_drag": 1,
                    "horizontal_wind": 1,
                    "sun_gravity": 1,
                    "moon_gravity": 1,
                    "solar_radiation_pressure": 1,
                    "solid_tides": 1,
                    "ocean_tides": 0,
                    "earth_albedo": 1,
                    "geopotential_model": "EIGEN-GL04C",
                    "atmospheric_model": "NRLMSISE-00",
                    "shadow_model": "conical",
                    "shadow_boundary_correction": 0,
                    "covariance_propagation": False
                }
            }
        }
    }

    request_neptune_simple, error = okapi_send_request(
        okapi_login,
        propagate_simple_request_body,
        'propagate-orbit/neptune/requests')
    if error['status'] == 'FATAL':
        print(error)
        exit()
    elif error['status'] == 'WARNING':
        print(error)
    print(str(request_neptune_simple))
    print("Getting OEM result from id ")
    error['web_status'] = 202
    counter = 0
    result_oem = {}
    while ((error['web_status'] == 202) and counter < 15):
        result_oem, error = okapi_get_result(okapi_login, request_neptune_simple,
        'propagate-orbit/neptune/results/{}/oem')
        if error['status'] == 'FATAL':
            print(error)
            exit()

    with open('sentinel-1b_neptune_gcrf.json', 'w') as outfile:
        json.dump(result_oem, outfile)
    print("Storing propagated data")

# Make sure we only propagate once, when the data is already there, we do not repeat the request. if not path.exists('sentinel-1b_neptune_gcrf.json'): okapi_login, error = okapi_init("https://api.okapiorbits.com/", "", "") if error.get('status', '') == 'FATAL': print(error) exit('Error during authentification.') #print(okapi_login) propagate_simple_request_body = { "orbit": { "type": "opm.json", "content": { "OPM_HEADER": { "CCSDS_OPM_VERS": 2 }, "OPM_META_DATA": { "OBJECT_ID": "41456", "OBJECT_NAME": "Sentinel-1b", "CENTER_NAME": "EARTH", "REF_FRAME": "GCRF", "REF_FRAME_EPOCH": "2000-01-01T00:00:00Z", "TIME_SYSTEM": "UTC" }, "OPM_DATA": { "EPOCH": "2019-06-24T22:59:42.000Z", "X": 5123.592195, "Y": -462.889158, "Z": -4865.970976, "X_DOT": -5.16151308, "Y_DOT": -1.00719510, "Z_DOT": -5.34379188, "MASS": 2151.964, "SOLAR_RAD_COEFF": 1.3, "DRAG_AREA": 6.2, "DRAG_COEFF": 2.2, "MANEUVERS": [ { "MAN_EPOCH_IGNITION": "2019-06-26T21:40:41.000Z", "MAN_DURATION": 23.875, "MAN_REF_FRAME": "RTN", "USER_DEFINED_MAN_A_1": -0.21367128E-08, "USER_DEFINED_MAN_A_2": 0.34196316E-06, "USER_DEFINED_MAN_A_3": 0.24668746E-07 }, { "MAN_EPOCH_IGNITION": "2019-06-26T22:29:45.000Z", "MAN_DURATION": 14.125, "MAN_REF_FRAME": "RTN", "USER_DEFINED_MAN_A_1": -0.14149311E-08, "USER_DEFINED_MAN_A_2": -0.37489161E-06, "USER_DEFINED_MAN_A_3": 0.26093520E-07 } ], }, }, }, "settings": { "type": "prop_settings.json", "content": { "propagation_end_epoch": "2019-07-02T00:59:42.000Z", "more": { "output_step_size": 100, "stormer_cowell_rel_tolerance": 1e-12, "stormer_cowell_abs_tolerance": 1e-13, "geopotential_degree_order": 70, "atmospheric_drag": 1, "horizontal_wind": 1, "sun_gravity": 1, "moon_gravity": 1, "solar_radiation_pressure": 1, "solid_tides": 1, "ocean_tides": 0, "earth_albedo": 1, "geopotential_model": "EIGEN-GL04C", "atmospheric_model": "NRLMSISE-00", "shadow_model": "conical", "shadow_boundary_correction": 0, "covariance_propagation": False } } } } request_neptune_simple, error = okapi_send_request( okapi_login, propagate_simple_request_body, 'propagate-orbit/neptune/requests') if error['status'] == 'FATAL': print(error) exit() elif error['status'] == 'WARNING': print(error) print(str(request_neptune_simple)) print("Getting OEM result from id ") error['web_status'] = 202 counter = 0 result_oem = {} while ((error['web_status'] == 202) and counter < 15): result_oem, error = okapi_get_result(okapi_login, request_neptune_simple, 'propagate-orbit/neptune/results/{}/oem') if error['status'] == 'FATAL': print(error) exit() with open('sentinel-1b_neptune_gcrf.json', 'w') as outfile: json.dump(result_oem, outfile) print("Storing propagated data")

The propagation results have been stored as a json document in sentinel-1b_gcrf.json. The file is read and the state vectors (GCRF) are stored in a list.

In [3]:

oem_body = {}
with open('sentinel-1b_neptune_gcrf.json', 'r') as file:
    oem_body = json.load(file)

propagated_gcrf_state_vectors = list()
for oem_data in oem_body['orbit']['OEM_META_DATA_DATA'][0]['DATA']:
    gcrf_state_vector = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
    gcrf_state_vector[0] = oem_data['X']
    gcrf_state_vector[1] = oem_data['Y']
    gcrf_state_vector[2] = oem_data['Z']
    gcrf_state_vector[3] = oem_data['X_DOT']
    gcrf_state_vector[4] = oem_data['Y_DOT']
    gcrf_state_vector[5] = oem_data['Z_DOT']
    propagated_gcrf_state_vectors.append(gcrf_state_vector)

file.close()

oem_body = {} with open('sentinel-1b_neptune_gcrf.json', 'r') as file: oem_body = json.load(file) propagated_gcrf_state_vectors = list() for oem_data in oem_body['orbit']['OEM_META_DATA_DATA'][0]['DATA']: gcrf_state_vector = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0]) gcrf_state_vector[0] = oem_data['X'] gcrf_state_vector[1] = oem_data['Y'] gcrf_state_vector[2] = oem_data['Z'] gcrf_state_vector[3] = oem_data['X_DOT'] gcrf_state_vector[4] = oem_data['Y_DOT'] gcrf_state_vector[5] = oem_data['Z_DOT'] propagated_gcrf_state_vectors.append(gcrf_state_vector) file.close()

Retrieving the Sentinel-1b validation data¶

Next the Sentinel 1b GNSS data is retrieved from ESA's service. The XML files are stored in the poeData directory. They are also only retrieved once.

In [4]:

url = 'http://aux.sentinel1.eo.esa.int/POEORB/'
dates = [
    '2019/07/15',
    '2019/07/16',
    '2019/07/17',
    '2019/07/18',
    '2019/07/19',
    '2019/07/20',
    '2019/07/21']
filenames = [
    'S1B_OPER_AUX_POEORB_OPOD_20190715T110507_V20190624T225942_20190626T005942.EOF',
    'S1B_OPER_AUX_POEORB_OPOD_20190716T110556_V20190625T225942_20190627T005942.EOF',
    'S1B_OPER_AUX_POEORB_OPOD_20190717T110525_V20190626T225942_20190628T005942.EOF',
    'S1B_OPER_AUX_POEORB_OPOD_20190718T110707_V20190627T225942_20190629T005942.EOF',
    'S1B_OPER_AUX_POEORB_OPOD_20190719T110400_V20190628T225942_20190630T005942.EOF',
    'S1B_OPER_AUX_POEORB_OPOD_20190720T110629_V20190629T225942_20190701T005942.EOF',
    'S1B_OPER_AUX_POEORB_OPOD_20190721T110515_V20190630T225942_20190702T005942.EOF']
dest = 'poeData'

if not os.path.exists(f'{dest}'):
    os.mkdir(f'{dest}')

i = 0
for filename in filenames:
    if not os.path.exists(f'{dest}/{filename}'):
        wget.download(f'{url}/{dates[i]}/{filename}', out=f'{dest}/{filename}')
    else:
        print('Nothing to download.')
    i=i+1

url = 'http://aux.sentinel1.eo.esa.int/POEORB/' dates = [ '2019/07/15', '2019/07/16', '2019/07/17', '2019/07/18', '2019/07/19', '2019/07/20', '2019/07/21'] filenames = [ 'S1B_OPER_AUX_POEORB_OPOD_20190715T110507_V20190624T225942_20190626T005942.EOF', 'S1B_OPER_AUX_POEORB_OPOD_20190716T110556_V20190625T225942_20190627T005942.EOF', 'S1B_OPER_AUX_POEORB_OPOD_20190717T110525_V20190626T225942_20190628T005942.EOF', 'S1B_OPER_AUX_POEORB_OPOD_20190718T110707_V20190627T225942_20190629T005942.EOF', 'S1B_OPER_AUX_POEORB_OPOD_20190719T110400_V20190628T225942_20190630T005942.EOF', 'S1B_OPER_AUX_POEORB_OPOD_20190720T110629_V20190629T225942_20190701T005942.EOF', 'S1B_OPER_AUX_POEORB_OPOD_20190721T110515_V20190630T225942_20190702T005942.EOF'] dest = 'poeData' if not os.path.exists(f'{dest}'): os.mkdir(f'{dest}') i = 0 for filename in filenames: if not os.path.exists(f'{dest}/{filename}'): wget.download(f'{url}/{dates[i]}/{filename}', out=f'{dest}/{filename}') else: print('Nothing to download.') i=i+1

Nothing to download.
Nothing to download.
Nothing to download.
Nothing to download.
Nothing to download.
Nothing to download.
Nothing to download.

Next the XML files are parsed and a state vector is extrated every 100 s (in line with the propagtion settings). Also the overlap between the files is ignored. Every state vector is transformed from ITRFto GCRF.

In [5]:

iso_time = None
last_iso_time = None
timediff = None
sentinel_gcrf_vectors = list()
i=0
for filename in tqdm(filenames):
    if os.path.exists(f'{dest}/{filename}'):
        dom = xml.dom.minidom.parse(f'{dest}/{filename}')
        itrf_vectors = dom.getElementsByTagName("OSV")
        for itrf_vector in itrf_vectors:
            state_vector = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
            for node in itrf_vector.childNodes:
                if (node.nodeName == "X"):
                    state_vector[0] = float(node.childNodes[0].nodeValue)/1000.0
                elif (node.nodeName == "Y"):
                    state_vector[1] = float(node.childNodes[0].nodeValue)/1000.0
                elif (node.nodeName == "Z"):
                    state_vector[2] = float(node.childNodes[0].nodeValue)/1000.0
                elif (node.nodeName == "VX"):
                    state_vector[3] = float(node.childNodes[0].nodeValue)/1000.0
                elif (node.nodeName == "VY"):
                    state_vector[4] = float(node.childNodes[0].nodeValue)/1000.0
                elif (node.nodeName == "VZ"):
                    state_vector[5] = float(node.childNodes[0].nodeValue)/1000.0
                elif (node.nodeName == "UTC"):
                    iso_time = datetime.fromisoformat(node.childNodes[0].nodeValue[4:])
                    if (iso_time != None and last_iso_time != None):
                        timediff = iso_time - last_iso_time
            if (timediff == None or timediff.total_seconds() >= 100.0):
                last_iso_time = iso_time
                gcrf_state_vector = itrf_to_gcrf(iso_time, state_vector)
                sentinel_gcrf_vectors.append(gcrf_state_vector)

print("Done reading files")

iso_time = None last_iso_time = None timediff = None sentinel_gcrf_vectors = list() i=0 for filename in tqdm(filenames): if os.path.exists(f'{dest}/{filename}'): dom = xml.dom.minidom.parse(f'{dest}/{filename}') itrf_vectors = dom.getElementsByTagName("OSV") for itrf_vector in itrf_vectors: state_vector = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0]) for node in itrf_vector.childNodes: if (node.nodeName == "X"): state_vector[0] = float(node.childNodes[0].nodeValue)/1000.0 elif (node.nodeName == "Y"): state_vector[1] = float(node.childNodes[0].nodeValue)/1000.0 elif (node.nodeName == "Z"): state_vector[2] = float(node.childNodes[0].nodeValue)/1000.0 elif (node.nodeName == "VX"): state_vector[3] = float(node.childNodes[0].nodeValue)/1000.0 elif (node.nodeName == "VY"): state_vector[4] = float(node.childNodes[0].nodeValue)/1000.0 elif (node.nodeName == "VZ"): state_vector[5] = float(node.childNodes[0].nodeValue)/1000.0 elif (node.nodeName == "UTC"): iso_time = datetime.fromisoformat(node.childNodes[0].nodeValue[4:]) if (iso_time != None and last_iso_time != None): timediff = iso_time - last_iso_time if (timediff == None or timediff.total_seconds() >= 100.0): last_iso_time = iso_time gcrf_state_vector = itrf_to_gcrf(iso_time, state_vector) sentinel_gcrf_vectors.append(gcrf_state_vector) print("Done reading files")

100%|██████████| 7/7 [00:26<00:00,  3.81s/it]Done reading files

Comparing the results¶

The state vectors (Sentinel 1b's GNSS' and the propagated ones) of each time step are transformed from GCRF to the satellite-centered coordinate system (RTN). The differences are then expressed in that coordinate system as

• radial,
• along-track (tangential) and
• cross-track (normal).

In [6]:

index = 0
rtn_vectors = list()
radial_differences = list()
tangential_differences = list()
normal_differences = list()
rtn_diff = np.array([0.0, 0.0, 0.0])
r_propagated_rtn = np.array([0.0, 0.0, 0.0])
r_sentinel_rtn = np.array([0.0, 0.0, 0.0])
for propagated_state_vector in tqdm(propagated_gcrf_state_vectors):
    if (index >= len(sentinel_gcrf_vectors)): break
    r_propagated = propagated_state_vector[0:3]
    r_sentinel = sentinel_gcrf_vectors[index][0:3]
    vecRTN = get_rtn_matrix(propagated_state_vector)
    r_propagated_rtn = np.dot(vecRTN, r_propagated)
    r_sentinel_rtn = np.dot(vecRTN, r_sentinel)
    rtn_diff = r_sentinel_rtn - r_propagated_rtn
    rtn_vectors.append(rtn_diff)
    radial_differences.append(r_sentinel_rtn[0] - r_propagated_rtn[0])
    tangential_differences.append(r_sentinel_rtn[1] - r_propagated_rtn[1])
    normal_differences.append(r_sentinel_rtn[2] - r_propagated_rtn[2])
    index += 1

index = 0 rtn_vectors = list() radial_differences = list() tangential_differences = list() normal_differences = list() rtn_diff = np.array([0.0, 0.0, 0.0]) r_propagated_rtn = np.array([0.0, 0.0, 0.0]) r_sentinel_rtn = np.array([0.0, 0.0, 0.0]) for propagated_state_vector in tqdm(propagated_gcrf_state_vectors): if (index >= len(sentinel_gcrf_vectors)): break r_propagated = propagated_state_vector[0:3] r_sentinel = sentinel_gcrf_vectors[index][0:3] vecRTN = get_rtn_matrix(propagated_state_vector) r_propagated_rtn = np.dot(vecRTN, r_propagated) r_sentinel_rtn = np.dot(vecRTN, r_sentinel) rtn_diff = r_sentinel_rtn - r_propagated_rtn rtn_vectors.append(rtn_diff) radial_differences.append(r_sentinel_rtn[0] - r_propagated_rtn[0]) tangential_differences.append(r_sentinel_rtn[1] - r_propagated_rtn[1]) normal_differences.append(r_sentinel_rtn[2] - r_propagated_rtn[2]) index += 1

100%|██████████| 6121/6121 [00:00<00:00, 12736.32it/s]