# Closing the Gap Between SGP4 and High-Precision Propagation via
Differentiable Programming

Giacomo Acciarini
University of Surrey
<EMAIL_ADDRESS>
&Atılım Güneş Baydin
University of Oxford
<EMAIL_ADDRESS>
&Dario Izzo
European Space Agency
<EMAIL_ADDRESS>
###### Abstract

The simplified general perturbations 4 (SGP4) orbital propagation is one of
the most widely used approaches used to compute a rapid and relatively
reliable prediction of the positions and velocity of objects orbiting Earth.
Over time, SGP models have undergone refinement to enhance their efficiency
and accuracy. Yet, they still fall short of the precision offered by numerical
propagators, which can predict the positions and velocities of space objects
in low-Earth orbit with significantly smaller errors. In this study, we
introduce a novel differentiable version of SGP4, named dSGP4. By
reconfiguring SGP4 into a differentiable program using PyTorch, we unlock its
differentiability for various space-related applications, encompassing
spacecraft orbit determination, state conversion, covariance similarity
transformation, state transition matrix computation, covariance propagation,
and more. A side advantage of this PyTorch-based model is its capacity for
"embarrassingly parallel" orbital propagation across batches of TLEs. It can
harness the processing power of CPUs, GPUs, and advanced hardware alike, to be
able to predict batches of satellites at future times in a distributed
fashion. Furthermore, the inherent differentiability of dSGP4 renders it
compatible with modern machine learning methodologies. Consequently, we
propose a novel orbital propagation paradigm, ML-dSGP4. In this paradigm, the
orbital propagator is imbued with neural networks. Through stochastic gradient
descent during the learning process, the inputs, outputs, and parameters of
this combined model can be iteratively refined to achieve precision surpassing
that of SGP4. Fundamentally, the neural networks function as identity
operators when the propagator adheres to its default behavior as defined by
SGP4. However, owing to the differentiability ingrained within dSGP4, the
model can be fine-tuned with ephemeris data to learn corrections to both
inputs and outputs of SGP4. This augmentation enhances precision while
maintaining the same computational speed of dSGP4 at inference time. This
empowers satellite operators and researchers, equipping them with the ability
to train the model using their specific ephemeris or high-precision numerical
propagation data.

_K_ eywords SGP4 $\cdot$ Orbital Propagation $\cdot$ Differentiable
Programming $\cdot$ Machine Learning $\cdot$ Spacecraft Collision Avoidance
$\cdot$ Kessler $\cdot$ Kessler Syndrome $\cdot$ AI for Space $\cdot$ Applied
Machine Learning for Space

##  1 Introduction

As the number of resident space objects (RSO) increases, space situational
awareness (SSA) becomes a critical aspect of today’s satellite operations,
involving the characterization, detection, and tracking of objects in orbit
around Earth. The conspicuous space debris population poses a threat to the
increasing number of operational satellites: collisions among space objects
have catastrophic consequences for the safety of surrounding objects and for
the future of space and the accessibility to certain orbital regimes. An
essential building block to ensure safe operations in space is orbital
propagation, which plays a key role in SSA since it is used to predict and
catalog the movement of RSOs in space. These propagators take an initial
state, describing the position and velocity of a space object, and integrate
them at future times using a model of the forces that act on the object in the
near-Earth environment. Some examples of these forces can include the
gravitational force exerted by the Earth and other Solar System bodies, as
well as the atmospheric drag and the solar radiation pressure. Depending on
the perturbations that are taken into account and on the numerical scheme used
to tackle the integration, the accuracy and speed of these methods in solving
the set of ordinary differential equations describing the satellites’ motion
can substantially vary.

Typically, we distinguish four types of orbital propagators [1]. General
perturbation techniques: these produce analytical solutions to the equations
of motion. These solutions are very fast to evaluate but are only an
approximation that is valid under some assumptions [2, 3]. Special
perturbation techniques: these consist of the numerical integration of the
equations of motion with the necessary perturbing accelerations. Typically,
these methods rely on some numerical scheme for the integration (e.g. Taylor
method, Runge-Kutta methods, Adams–Bashforth-Moulton methods, etc.), which
leads to slower but more accurate solutions w.r.t. general perturbation
techniques [4]. Semianalytical techniques: these try to combine the accuracy
of numerical methods with the speed of analytical techniques. Often, they
consist of numerical schemes that solve a set of equations that are derived
from the original equations of motion under certain assumptions, which enable
faster integration [5]. Hybrid techniques: these are an emergent trend that
tries to combine the speed and accuracy of the above method by using
statistical and machine learning methods. In this paradigm, the prediction of
general perturbation methods is corrected via a machine learning and/or
statistical model, to more closely resemble the accuracy of the numerical
propagator prediction [6, 7].

During the 1960s, the US Air Force initiated an endeavor to develop orbital
propagation models that could provide reliable and rapid predictions of the
positions of resident space objects based on initial position and velocity
data. This initiative led to the creation of what are known as simplified
perturbation models. One of the most popular general perturbation techniques
is the Simplified General Perturbations 4 (SGP4) method, which belongs to the
class of the simplified perturbations models (SGP, SGP4, SDP4, SGP8, SDP8)
[8]. This consists of a mathematical model with a long history of development
(starting from the 1960s) used to predict the position and velocity of an
Earth-orbiting object [9, 10], although the first release of the refined SGP4
propagator only happened in the 1970s [11]. Originally, it was only applicable
to satellites whose orbital period is less than 225 minutes, although an
extension of this model for higher periods has also been developed, which is
known as SDP4. It owes its popularity to the fact that simplified perturbation
models are the only ones that support two-line elements (TLEs) data in their
original format: these are the only openly available ubiquitous source of
information for the position and velocity of satellites. TLE is a data format
that describes in two lines of 69 characters each, the position, velocity,
characteristics (e.g. drag term), and other information (e.g. satellite ID,
etc.) of RSOs around the Earth. United States Space Force routinely tracks
objects in space and publicly releases their TLE elements through the Space-
Track website 111https://www.space-track.org/auth/login, date of access: March
2023.222https://celestrak.org/NORAD/elements/, date of access: March 2023..
Whilst TLE data are released publicly, they are created from a mixture of
analytical theories, observations and coordinate systems that are not revealed
[10]. For example, SGP4 uses power density functions to estimate the drag
force that acts on the satellite motion [12, 13]: these require the knowledge
of a term that is related to the ballistic coefficient, called Bstar
parameter. Since a mix of simplified perturbations theory and observational
data (processed through batch least squares) is used to produce a single TLE,
the Bstar parameter often absorbs force model errors [1].

All these unrevealed underlying processes necessary to produce a TLE make the
transformation from TLE to osculating elements not possible without
errors/assumptions and their use is therefore recommended only with a
simplified perturbations propagator. Typical accuracies of these propagators
with TLE data are in the order of kilometers for predictions within a week,
compared to the ground truth states of the object [14, 15, 16, 17]. For this
reason, when osculating orbital elements data are available, special
perturbation techniques are often preferred in applications that require high
accuracies (e.g. collision avoidance). This has motivated several researchers
to improve the SGP4 prediction via statistical and machine learning methods,
to make it competitive with special perturbation techniques while retaining
its computational advantages. In general, there is a plethora of studies on
orbit prediction enhancements using Gaussian processes, state vector machines,
and other machine learning techniques [18]. In the context of SGP4, most of
the proposed techniques focus on creating datasets of errors between SGP4
predicted states and ephemeris data (or high accuracy numerical propagators),
which are then used to train machine learning models that can learn to correct
these errors [19]. Recently, some researchers have also devised a fusion
strategy between orbit prediction and Bayesian inference methods to enhance
the SGP4 predictions from TLE data [20].

Motivated by these works, and by the lack of a general framework that allows
improving SGP4 predictions integrating highly accurate data and/or
predictions, we developed a differentiable SGP4 (dSGP4) version on Python,
which we make available through an open-source repository on the European
Space Agency Github organization 333https://github.com/esa/dSGP4, date of
access: March 2023.. Our contributions consist of first re-writing the SGP4
propagator to make it compatible with modern programming techniques, which
enable automatic differentiation and "embarrassingly parallel" computations
via CPU and/or GPU. We then show that having an SGP4 propagator that is
differentiable with respect to its inputs opens the possibility of using the
Jacobian for several applications (e.g. state covariance propagation, state
transition matrix computation, gradient-based optimization). This has also
consequences in many fields related to SSA, such as orbit determination, among
others. To underscore the significance of employing dSGP4 in practical
contexts, we reference a work by Vallado on SGP4 [21], where he suggests that
while elegant operational systems typically opt for numerical integration, a
straightforward analytical propagation theory like SGP4 might justify
employing finite differencing. According to him, the effort and resources
invested in developing partial derivatives for various data types and state
forms to achieve computational efficiency are often outweighed by the
simplicity of incorporating an additional iteration or two with a finite
differencing approach. Through the use of dSGP4, we demonstrate that the
computation of partial derivatives can be seamlessly executed without
incurring any extra time or financial costs, thanks to automatic
differentiation. Additionally, this approach mitigates local truncation errors
associated with finite differencing.

Finally, we also discuss the applications of such a differentiable SGP4 model
in a new hybrid propagator paradigm, where the machine learning model is made
of an initial neural network that ingests the data, which is then passed
through dSGP4, and finally through another neural network. In this way, owing
to the differentiability of dSGP4, it is possible to take gradients of the
entire model, therefore applying standard stochastic gradient descent
techniques and training dSGP4 together with neural network models. This
framework enables backpropagating through the whole model, and learning
corrections of SGP4 inputs and outputs to match the predicted state: in
practice, this means that the neural networks, if successful, learn the
differential correction terms needed to correct the low-accuracy predictions
of SGP4, to better match either observations or high-fidelity propagations,
and retain the advantageous SGP4 speed. In this way, for the first time, both
inputs and outputs (and potentially hyperparameters of SGP4) can be learned
for the given prediction task, to maximize the accuracy of the predicted
orbital states, while maintaining a fast and parallelizable propagator.

The paper is organized as follows: in Section 2, we discuss relevant
background on the SGP4 model and differentiable programming. Then, in Section
3, we discuss the features of the propagator, including automatic
differentiation and batch mode. In Sec. 4, we discuss its machine learning
integration, by both outlining the general framework and discussing simulation
results that include the use of the machine learning dSGP4 framework to
improve the SGP4 predictions on a batch of SpaceX Starlink satellites. While
in Sec. 5, we introduce the use of dSGP4 in the context of satellite orbit
determination. Other applications including gradient-based optimization,
osculating to TLE elements conversion, covariance propagation, similarity
transformation, state transition matrix computation, and batch propagation are
included as tutorials of the open-source library and released with dSGP4.
Finally, we also conduct an experiment on a group of SpaceX Starlink
satellites where dSGP4 and ephemeris data are used to train the ML-enhanced
hybrid propagator, to improve the baseline SGP4 predictions. Finally, in
Section 6, we highlight the conclusions of our work.

##  2 Related Work

###  2.1 SGP4 Model

The Simplified General Perturbations-4 (SGP4) [10] software is an orbit
propagation model and it is the only one recommended to be used in combination
with two line elements (TLEs) orbit data. Two-line elements are a very popular
data format used for space object cataloging. They encapsulate information
associated with the satellite position and velocity (using general
perturbation mean orbital elements) at a given epoch, together with some
identifiers for the satellite (e.g. satellite number, international
designator, etc.) [1]. Six quantities in this format (the mean motion $n$, the
eccentricity $e$, the inclination $i$, the right ascension of the ascending
node $\Omega$, and the mean anomaly $M$) are associated with the position and
velocity in space of the object, while the other three (i.e., the mean motion
rate $\dot{n}$, mean motion acceleration $\ddot{n}$, drag-like parameter
$B^{*}$) are used to describe the effect of perturbations in the satellite
motion, and the time in UTC is used to indicate the epoch at which the TLE is
generated. All these two-line elements variables are directly used by the SGP4
model, which then produces ephemerides for the resident space object at the
required time, in the true equator, mean equinox (TEME) coordinate system [1].

###  2.2 Differentiable Programming

Differentiable programming is a programming paradigm in which models can be
built through a computer program and then their parameters can be
differentiated through automatic differentiation (AD). This allows users to
perform gradient-based optimization of the model outputs with respect to its
parameters. This means that users can leverage gradient-based techniques to
optimize large and complex models with many parameters, adjusting those
parameters in order to minimize the deviation of the model output with respect
to ground truth data. With the advent of artificial intelligence and neural
network models, this paradigm has become increasingly popular, and many
scientific fields have successfully used differentiable programming in several
areas: from computer science to robotics, ray tracing, image processing,
probabilistic programming, electronics, etc. [22, 23, 24, 25].

By investigating the SGP4 formulation, it is apparent that it can be described
as a set of continuous and differentiable functions that given several input
parameters, returns a 6-dimensional vector, and that this can be used in some
applications such as orbit determination [26, 4]. However, there is currently
no open-source SGP4 version formulated as a differentiable program, on which
automatic differentiation can be applied not only for direct gradient-based
optimization techniques such as orbit determination and model calibration but
also to combine deep learning with SGP4. In fact, one of the current issues of
the application of SGP4 in SSA is that while it allows performing fast orbital
prediction and conjunction assessment for the full catalog, its results are
not particularly accurate, which forces many operators to refine their
analysis with more precise (but slower) numerical integration schemes that
represent the underlying dynamical model more accurately [21, 16, 27]. Some
researchers have attempted to trade off the SGP4 speed with numerical
integration accuracy [28] or to take advantage of neural networks to correct
the output of SGP4 [7, 19]. However, to our knowledge, there have been no
attempts to combine SGP4 and neural networks within the same model to maintain
the same speed as SGP4 while correcting the model’s inputs, outputs, and
hyperparameters to better match accurate data or high-fidelity orbital
propagators. In this work, we explore this aspect and show that formulating
SGP4 as a differentiable program paves the way to easy integration with neural
networks.

##  3 dSGP4: an Open-Source Differentiable SGP4 Program

###  3.1 dSGP4 Overview: Automatic Differentiation and Batch Mode

Automatic differentiation is a technique used for computing and evaluating
derivatives of functions. It can be applied to any function that can be
represented as a sequence of elementary operations, whose derivatives are
known [29]. The basic idea of AD is to build a graph of elementary operations
with known derivatives, which is then used (either in forward or reverse mode)
to precisely compute the derivative at a given point. Before automatic
differentiation, derivatives could either be computed manually, with a time-
consuming and error-prone procedure, through symbolic differentiation, or
through numerical differentiation. This last alternative is not ideal due to
the presence of truncation and round-off errors, which might substantially
jeopardize the accuracy of the found derivatives, if not chosen carefully. On
the other hand, symbolic differentiation would allow to access the derivative
of the function for the entire domain, but it can result in the size of the
expressions growing exponentially (something known as "expression swell"), as
the number of variables and operations is increased, due to the product rule
of derivative [30]. Automatic differentiation solves these issues, by
returning the exact value (up to machine precision) of the derivatives of a
computer program evaluated at a certain point, without incurring exponentially
growing expressions [31]. The drawback, besides the fact that it is not
applicable for non-differentiable functions, is that this is only the
derivative evaluated at a given point, and one does not have access to the
symbolic form of the derivative across the entire domain. The
differentiability constraint implies that the computer program needs to comply
with certain requirements and does not admit discontinuities, singularities,
or non-differentiable features. In our case, we have implemented dSGP4 using
the Python library Pytorch [32], which is one of the most popular open-source
ML libraries that supports AD in use in deep learning research and
applications. To ensure the fidelity of dSGP4, our approach encompassed
several key verification steps. First of all, we validated the dSGP4 model’s
outputs against the "Revisiting Spacetrack Report #3" SGP4 model [10]. This
involved the design and execution of continuous integration tests (that can be
found in the open-source library) to confirm that all mirrored functionalities
faithfully reproduced the behavior of the original
implementation444https://github.com/esa/dSGP4/tree/master/tests, date of
access: February 2024.. As a part of these tests, we also verified that we
could reproduce the output of the Space Track official SGP4
release555https://www.space-track.org/documentation#/sgp4, date of access:
February 2024.: this was done by propagating several different TLEs at
different future times (between a few hours to a few days), and checking that
the results of the propagation match up to below 1 mm in position and 1 mm/s
in velocity the official release results. Our verifications also extended to
the examination of automatic differentiation outputs for the input variables.
To ensure correctness, we compared (up to local truncation errors caused by
the finite difference technique) that the results obtained from dSGP4
automatic differentiation functionality match the ones found via finite
differencing of the original SGP4.

As previously discussed, dSGP4 also offers an additional capability, as we
show in Figure 1. While retaining the traditional functionality of predicting
object states from a TLE for specific times or multiple times, dSGP4 offers a
distinctive feature called batch mode. In this mode, batches of TLEs can be
provided to the model, which in turn returns the corresponding batch of
Cartesian states. This expanded functionality has also undergone comprehensive
testing within our library.

Figure 1: Schematic illustration of batch mode supported by dSGP4 (on the
right) vs the standard SGP4 implementation (on the left).

Differentiable physics simulators have the key advantage that can be run both
in forward and backward mode. When the forward mode is active, the simulator
has the commonly known behavior of returning a certain output from the
provided inputs: however, when the backward mode is activated, the model can
also return the gradients of the outputs with respect to any input and
parameter of the model. This allows the integration of this simulator with
gradient-based methods, where the gradient information can be leveraged to
learn and control the physical system. This approach enables a variety of
inverse problems that involve orbital propagation, which we will also discuss
in the coming sections (e.g. satellite orbit determination). The interested
reader can also find in the dSGP4 website 666https://esa.github.io/dSGP4, date
of access: October 2023. a plethora of tutorials on how to apply the dSGP4
framework to common gradient-based tasks for SSA applications (e.g. state
transition matrix and sensitivity matrix computation, covariance propagation,
and gradient-based optimization). The simulator is treated as a differentiable
function, and a certain objective function is minimized, which depends on the
simulator. Mathematically this is formulated as a minimization problem of a
certain objective function: $J(\boldsymbol{\mathrm{\theta}})$, where
$\boldsymbol{\mathrm{\theta}}$ is an $n$-dimensional vector of model
parameters. The goal is to update the parameters in such a way that
$J(\boldsymbol{\mathrm{\theta}})$ reaches the minimum value: this is done in
the direction opposite of the gradient of the objective function w.r.t. the
parameters. In its simplest (i.e., vanilla gradient descent) form this can be
written as
$\boldsymbol{\mathrm{\theta}}_{i+1}=\boldsymbol{\mathrm{\theta}}_{i}-\eta\cdot\nabla_{\theta}J(\boldsymbol{\mathrm{\theta}}_{i})$,
where $\eta$ is known as learning rate [33]. Therefore, it is clear that by
formulating the orbital propagator as a differential program, its integration
with any gradient-related task, such as vanilla gradient descent is possible.
This also has the direct consequence of being able to use these gradients
within root-finding algorithms that leverage gradient information.

Beyond the advantages of parallel computing, gradient computation, and more,
the differentiability of SGP4 opens the door to a new propagator paradigm.
This paradigm involves integrating a low-fidelity propagator with neural
networks to jointly correct inputs and outputs of SGP4, in a supervised
learning prediction task enhancing model accuracy while retaining
computational efficiency. To our knowledge, this is the first time in which an
orbital propagator is backpropagated together with other machine learning
architectures, to combine the stochastic gradient descent’s capability of
reducing the prediction error, with the physics-based prediction of the model.

###  3.2 Performance Profiling

As already discussed in the previous section, an important advantage of having
a differentiable SGP4 program written in PyTorch is the possibility to
massively distribute the propagation of TLEs, which we refer to as batch-mode
propagation, leveraging either multiple CPUs or GPU acceleration (through
CUDA). In this way, we enable the possibility to have embarrassingly parallel
computations of multiple TLEs at multiple times, through hardware
parallelization. To quantitatively assess this capability, we display in Tab.
1 the performance of the simple SGP4, against dSGP4 on 1 CPU and GPU. We
compare the results using float64 numbers (this is not ideal for the GPU case,
but it might be important for orbital propagation applications). In terms of
CPU processor, we have used the MacBook M1 Pro processor, while in terms of
GPU, we have used the NVIDIA V100 Tensor Core GPU777https://www.nvidia.com/en-
us/data-center/v100/, date of access: August 2023.. As we see from the table,
already for batches of size 1,000, dSGP4 offers a significant speed-up, when
run in batch mode, on both CPU and GPU, compared to the standard SGP4 program.
In particular, as we see, the GPU implementation can reach a factor of
85$\times$ speed-up compared to the standard SGP4. This level of speed-up
factor holds promise for applications such as all-versus-all conjunction
screening, where computational speed-up can yield substantial benefits.

Table 1: Execution Times for SGP4, dSGP4 on CPU, and dSGP4 on GPU Model | Run Time (ms) per Event
---|---
Single Event | 10,000 | 100,000 | 1,000,000 | 10,000,000
SGP4 | 0.002 |  | 21.11 |  | 211.1 |  | 2,111 |  | 21,110 |
dSGP4 (1 CPU) | 0.67 |  | 4.52 |  | 44.6 |  | 460 |  | 7,920 |
dSGP4 (GPU) | 4.22 |  | 7.8 |  | 8 |  | 29 |  | 249 |

##  4 Machine Learning dSGP4

###  4.1 Methods

As detailed in Section 3.1 the propagator can be seamlessly integrated into
complex machine learning pipelines, offering a differentiable program for end-
to-end training using techniques like stochastic gradient descent. While
general perturbation propagators like SGP4 boast fast propagation times, their
low accuracy limits their applicability, especially in scenarios like
collision avoidance.

This section introduces ML-dSGP4, a paradigm where a differentiable propagator
(dSGP4) is integrated with two neural networks to correct its inputs and
outputs. In Figure 2(a) and 2(b), respectively, both standard SGP4 and dSGP4
propagations are illustrated. While in Figure 2(b) and 2(d), respectively, the
combination between ML and both SGP4 and dSGP4 are shown.

|
---|---
|

(a) SGP4 model.

(b) dSGP4 model.

(c) ML-SGP4 model.

(d) ML-dSGP4 model.

Figure 2: Schematic illustration of SGP4, ML-SGP4, dSGP4, and ML-dSGP4 models.

The ML-dSGP4 model represents a new hybrid differentiable propagator that
retains the computational advantages of SGP4 while enabling training and
improvement with accurate data. The model can be configured to return either
SGP4 predictions or activate input and output corrections using trained neural
networks.

To better discuss the ML-dSGP4 setup shown in Figure 2, let us call
$\boldsymbol{\mathrm{x}}_{0}$ the input of ML-dSGP4 (i.e., the TLE elements),
$t$ the propagation time, and $\boldsymbol{\mathrm{x}}$ its output, and let us
assume that the objective is to minimize the error between the ML-dSGP4
propagation and the propagation of a high precision propagator,
$\boldsymbol{\mathrm{x}}_{\textrm{HPOP}}$ (i.e., a numerical integrator that
accounts for all forces involved in the satellite motion). In this case, the
objective is to correct ML-dSGP4 to better match the accurate propagator,
hence, the function to minimize will be a function of the state, time, and
parameters (both of the networks and also of dSGP4). Without loss of
generality, let’s consider as an objective function the squared norm between
the target state and the dSGP4 output, then the objective is to find the
parameters (i.e., weights and biases) of the two networks, as well as the
parameters of SGP4 that minimize the loss:

$\mathcal{P}:\begin{cases}\mbox{given:}&\boldsymbol{\mathrm{x}}_{0},t\\\
\mbox{find:}&\boldsymbol{\mathrm{\theta}}_{1},\boldsymbol{\mathrm{\theta}}_{\textrm{SGP4}},\boldsymbol{\mathrm{\theta}}_{2}\\\
\mbox{s.t.
min:}&J(\boldsymbol{\mathrm{x}}(\boldsymbol{\mathrm{\theta}}_{1},\boldsymbol{\mathrm{\theta}}_{\textrm{SGP4}},\boldsymbol{\mathrm{\theta}}_{2}))=||\boldsymbol{\mathrm{x}}-\boldsymbol{\mathrm{x}}_{\textrm{HPOP}}||^{2}\end{cases}\text{,}$
(1)

where $\boldsymbol{\mathrm{\theta}}_{1}$ indicates the parameters of the input
neural network, $\boldsymbol{\mathrm{\theta}}_{2}$ of the output one, and
$\boldsymbol{\mathrm{\theta}}_{\textrm{SGP}}$ the SGP4 parameters. Owing to
the differentiability of the neural networks and dSGP4, we can use gradient
descent to find the local minimum of the loss function by iteratively updating
the parameters of the neural networks and dSGP4:

$\displaystyle\begin{split}&\boldsymbol{\mathrm{\theta}}_{1,i+1}=\boldsymbol{\mathrm{\theta}}_{1,i}-\eta_{1}\nabla_{\boldsymbol{\mathrm{\theta}}_{1}}J\\\
&\boldsymbol{\mathrm{\theta}}_{\textrm{SGP4},i+1}=\boldsymbol{\mathrm{\theta}}_{\textrm{SGP4},i}-\eta_{\textrm{SGP4}}\nabla_{\boldsymbol{\mathrm{\theta}}_{\textrm{SGP4}}}J\\\
&\boldsymbol{\mathrm{\theta}}_{2,i+1}=\boldsymbol{\mathrm{\theta}}_{2,i}-\eta_{2}\nabla_{\boldsymbol{\mathrm{\theta}}_{2}}J\\\
\end{split}\text{.}$ (2)

This setup can be easily extended for multiple propagation points, for
instance, by taking the mean squared error as a loss function, and by setting
up a stochastic gradient descent algorithm to update the parameters. The
advantage of ML-dSGP4 lies in its capacity to learn both the initial neural
network weights and SGP4 parameters, therefore offering a flexible solution.
This differs from ML-SGP4, where only the output of the SGP4 propagator is
corrected, as explored in prior research [7]. In the next section, we will
delve into some numerical experiments to showcase the use of ML-dSGP4 and to
highlight the comparison between ML-dSGP4 and the ML-SGP4.

###  4.2 Experiments

This section presents experiments conducted to assess the accuracy of the ML-
dSGP4 propagator and compare it with the uncorrected SGP4 baseline and a
neural network corrector that only adjusts SGP4 outputs [7].

The dataset comprises SpaceX Starlink Two-Line Elements (TLEs) and ephemeris
data available on Space-Track888https://www.space-track.org/, date of access:
April 2023.. The ephemeris data includes a set of 1,519 Starlink satellites,
in the period from the 18th of February 2023, at 11:38:42 UTC, to the 21st of
February 2023, at 11:38:42 UTC, with a time resolution of 60 seconds. This
data is provided by SpaceX using their numerical orbital propagator. We also
collect the TLE data for the same objects during the period of the 17th and
18th of February, before the date of the first available ephemeris. We then
use dSGP4 to propagate those TLEs at the same times as the available ephemeris
data, and we treat each TLE as a separate observation. This generates a
dataset of as many entries as the number of TLEs for each object, times the
number of time stamps in the ephemeris data, multiplied by the total number of
satellites: that is, a dataset with 6,563,599 rows and 6 columns (three for
the position vector and three for the velocity vector).

We split the training, validation, and test sets at the TLE level, making sure
that the data propagated from a certain TLE belongs to the same set and is not
shared among multiple sets: this is to avoid the model essentially learning to
interpolate the data. We use 69% of the dataset as training, 16% as
validation, and 15% as testing. Moreover, we make sure that the test set
covers ranges of altitude that are unseen in the training and validation set:
this resulted in 225 Starlink satellites below 539 km (while all the other
satellites above that altitude have been used in training and validation). We
train the model with the architecture shown on the right side of Figure 2, by
using two neural networks (each with three hidden layers of 35 neurons each),
using 30 epochs, a learning rate of 3$\times 10^{-3}$, mean squared error
(MSE) as a loss, and Adam as optimizer. Moreover, we also learn the correction
factors for all elements of the inputs and outputs of dSGP4 (this is a set of
12 parameters in total). The total number of learnable parameters was 3,454.

We compare this model with a neural network output corrector that uses
slightly more learnable parameters (i.e., 3596 parameters: to have a fair
comparison): resulting from a feed-forward neural network with four hidden
layers of size 32, plus the six correction factors for the output corrections.
In both this and the ML-dSGP4 case, the networks only learn corrections to the
given inputs and/or outputs.

In Figure 3, we show the loss of the training, validation, and test data as a
function of epochs for both the ML-dSGP4 and the neural network corrector
model. As we observe, in all three datasets, both models manage to outperform
the SGP4 baseline model, and ML-dSGP4 manages to achieve a higher accuracy
than the neural network corrector, thereby confirming the added advantages of
the differentiable setup that can also learn input corrections.

In Table 2, we show the final losses of both models, compared to the baseline
results: to produce the table, the best models at validation were chosen for
the neural network corrector and ML-dSGP4.

Figure 3: Training, validation, and test losses mean squared error as a
function of epochs.

As already mentioned, the test set was in a set of altitudes (and therefore
also dynamics) completely unseen during training or validation. Nonetheless,
the ML-dSGP4 model manages to achieve accuracies in the test set that are in
the same range of values as the validation ones: this confirms that the model
successfully starts to learn a differential correction between the low-
accuracy SGP4 and the high-precision orbital propagation model. This paves the
way for a mixed strategy where low-precision propagation with dSGP4 can be
combined with high-precision data (or orbital propagation data) to learn
correction models that then enable to perform fast and more accurate
predictions of orbital states. Besides, this is also combined with the
possibility of performing these computations in batches and leveraging
parallelization via CPU or GPU. Finally, in this experiment we did not learn
any internal dSGP4 model correction factors: however, the same framework can
also be leveraged to allow dSGP4 internal parameters to be learned, thereby
providing further prediction accuracy enhancements.

|  | State [-] | Pos [km] | Vel [km/s]
---|---|---|---|---
| Training | 5.41179$\times 10^{-5}$ | 46.6 | 0.0508
SGP4 (baseline) [11] | Validation | 0.000404 | 100.32 | 0.1097
| Test | 0.0004153 | 107.42 | 0.1176
| Training | 2.69409$\times 10^{-5}$ | 33.36 | 0.0363
ML-SGP4 [7] | Validation | 0.000323 | 81.27 | 0.0888
| Test | 0.000326 | 88.66 | 0.0969
| Training | 1.75807$\times 10^{-5}$ | 25.01 | 0.0276
ML-dSGP4 | Validation | 0.000262 | 64.02 | 0.0704
| Test | 0.000256 | 71.25 | 0.0783
Table 2: Baseline, ML-SGP4, and ML-dSGP4 RMSE values for the normalized state,
the position, and velocity after 25 epochs, for training, validation, and test
datasets.

The objective of this experiment is to show that even with only three days of
ephemeris data, the ML-dSGP4 model can be leveraged to increase the accuracy
of the standard SGP4 propagator while retaining its computational advantages.
Feeding more ephemeris data, the model can further enhance its accuracy,
potentially starting to be competitive with high-precision numerical
propagators in terms of precision, while retaining the advantages of general
perturbation techniques in terms of speed. This might also be very appealing
for operators of constellations that possess a large amount of high-precision
data.

##  5 Satellite Orbit Determination

###  5.1 Methods

Satellite orbit determination consists of the accurate determination of
spacecraft states starting from a dynamical model of the forces acting on the
satellite (e.g. SGP4) and available observations (e.g. obtained through
ground-based instruments). The objective is to combine a simulated environment
and raw data to provide a good estimate of the satellites’ position and
velocity. The orbit determination processes both in general and in relation to
SGP4 are detailed in Vallado’s publications [1, 21].

In the aforementioned works, either finite differences or analytical schemes
are proposed as a means to find the Jacobian matrix needed for solving the
least squares problem. However, in our case, we will show that this is not
necessary, since the differentiability of dSGP4 enables the computation of
those quantities via AD, which turns out to be very advantageous, as it avoids
extra computation and/or additional errors.

The detailed problem setup is that we have $N$ observations,
$\\{\hat{\boldsymbol{\mathrm{x}}}_{i}\\}_{i=1}^{N}$, at different observation
times, $\\{t_{i}\\}_{i=1}^{N}$. Our objective is to find the state at a
certain target time, $t_{T}$, that best fits the set of observations.
Effectively, this means finding the state at time $t_{T}$:
$\boldsymbol{\mathrm{x}}_{T}$, which minimizes the (weighted) squared error
between the observed and modeled states at the observation times. This problem
can be described as follows:

$\mathcal{P}:\begin{cases}\mbox{given:}&\\{\hat{\boldsymbol{\mathrm{x}}}_{i},t_{i}\\}_{i=1}^{N},\
t_{T}\\\ \mbox{find:}&\boldsymbol{\mathrm{x}}_{T}\\\ \mbox{s.t.
min:}&\sum_{i=1}^{N}|(\boldsymbol{\mathrm{x}}_{i}(\boldsymbol{\mathrm{x}}_{T},t_{T})-\hat{\boldsymbol{\mathrm{x}}}_{i})w_{i}|^{2}\end{cases}\text{,}$
(3)

where $w_{i}$ are the observation weights. The modeled states
$\boldsymbol{\mathrm{x}}_{i}$ at the observed times $t_{i}$, can be found via
the nonlinear propagation of $\boldsymbol{\mathrm{x}}_{T}$ at $t_{T}$, until
$t_{i}$, using dSGP4:

$\boldsymbol{\mathrm{x}}_{i}=\textrm{dSGP4}(\boldsymbol{\mathrm{x}}_{T},t_{i}-t_{T})$
(4)

It is straightforward to observe that the above problem is equivalent to a
weighted least-squares problem, whose minimum can be sought starting from an
initial guess and setting up an update rule that depends on the Jacobian of
the modeled states at observation times, with respect to the state at the
target time [4].

Hence, starting from an initial guess $\boldsymbol{\mathrm{x}}_{T}^{0}$, we
can setup an iteration scheme:

$\boldsymbol{\mathrm{x}}_{T}^{k+1}=\boldsymbol{\mathrm{x}}_{T}^{k}+\delta\boldsymbol{\mathrm{x}}\text{,}$
(5)

where:

$\delta\boldsymbol{\mathrm{x}}=-(A^{T}WA)^{-1}A^{T}W\boldsymbol{\mathrm{b}}\text{,}$
(6)

with $W$ being the weight matrix,
$\boldsymbol{\mathrm{b}}_{i}=(\boldsymbol{\mathrm{x}}_{i}^{k}-\hat{\boldsymbol{\mathrm{x}}}_{i})$
and $A_{i}$ is the Jacobian of the $i$-th modeled state with respect to the
target state:

$A_{i}=\dfrac{\partial\boldsymbol{\mathrm{x}}_{i}}{\partial\boldsymbol{\mathrm{x}}_{T}}\bigg{|}_{\boldsymbol{\mathrm{x}}_{T}^{k}}\text{.}$
(7)

Owing to the differentiability of dSGP4, the partial derivative terms in Eq.
(7) can now be exactly computed via automatic differentiation. In case more
than one observations are available, then the $A$, $W$ matrices and
$\boldsymbol{\mathrm{b}}$ vector can be found by stacking each observation:
$A^{T}WA=[A_{0}^{T}WA_{0},\dots,A_{N}^{T}WA_{N}]^{T}$,
$A^{T}W\boldsymbol{\mathrm{b}}=[A_{0}^{T}W\boldsymbol{\mathrm{b}},\dots,A_{N}^{T}W\boldsymbol{\mathrm{b}}_{N}]$.
In case no weights are assigned for the observations, then the weight matrix
corresponds to the identity matrix.

###  5.2 Experiments

Using the procedure highlighted in Section 5, we setup an experiment to show
the usage of dSGP4 in an orbit determination problem. In our analysis, we
assume to have a set of TLEs as observations, and we use SGP4 as a nonlinear
orbital propagation model: our goal is to set up the problem as a nonlinear
least squares problem, in which we want to find a TLE for a given target epoch
that best fits the observations, using the SGP4 model. This experiment is also
documented in one of the tutorials of the open-source package. As we already
mentioned, three essential components of satellite orbit determination are a
set of observations, a target time at which we desire to extract the
satellites’ state, and an initial guess to setup the iterative procedure. We
download the set of TLEs corresponding to Sentinel-2A satellite in the period
from 1st of June to 8th of June 2022, from the Space-
Track999https://www.space-track.org/, date of access: April 2023. public
catalog. This corresponds to a set of 26 TLEs, with a maximum time difference
of about 10 hours, and a minimum time difference of about 3 hours: we use
these TLEs as observations. Then, we choose the target epoch, at which we
would like to extract the best fit for the TLE, one day after the final
observation, which corresponds to the 8th of June 2022, at 21:25:48.337248
UTC. The initial guess for the TLE elements is found from the average state of
the 26 TLEs propagated at the target epoch via dSGP4. This corresponds to the
following initial state: $\boldsymbol{\mathrm{r}}=[4.5995\times
10^{2},2.4107\times 10^{3},6.7259\times 10^{3}]^{T}$ km,
$\boldsymbol{\mathrm{v}}=[4.3698,5.5925,-2.2983]^{T}$ km/s, which results in
the nominal target TLE (i.e., $\boldsymbol{\mathrm{y}}_{T}^{0}$) shown in Tab.
3. By propagating the nominal TLE at the observation times and comparing the
norm of the residuals w.r.t. the actual observation, we get a total residual
norm in the state of 2915.6622, with residual norms ranging from a maximum of
900 km in position and 1 km/s in velocity, to a minimum of 100 km in position
and 0.1 km/s in velocity for the closest in time. We use this nominal TLE as
an initial guess for the differential corrector. After convergence, the TLE in
Tab. 4 is found, with a total residual norm in the state (norm of the sum of
all residuals) of 1.2112, which corresponds to residuals that range from a
minimum of 12.6 m in position and 0.5 m/s in velocity, and to a maximum of
513.6 m in position and 1.4 m/s in velocity. The difference between this
example and previous satellite orbit determination applications shown in
literature is that in this case, we used dSGP4 automatic differentiation to
precisely compute the partial derivatives needed to build the $A$ matrix of
Eq. (7), therefore avoiding the local truncation errors due to finite
differencing [21].

1 40697U 15028A 22159.89292057 .00000111 00000-0 59112-4 0 9998
---
2 40697 98.5685 234.7917 0000491 348.2227 119.9277 14.31085333362518
Table 3: initial nominal TLE. 1 40697U 15028A 22159.89292057 .00000111 00000-0
-19814-4 0 9998
---
2 40697 98.5716 234.7932 0000011 348.1175 120.0022 14.30817171362512
Table 4: TLE after convergence.

##  6 Conclusions

In this work, we introduce dSGP4, a differentiable version of the SGP4 model
implemented in PyTorch, representing a significant advancement in space
situational awareness. Our research work yields several considerations.

First of all, the differentiability features of the model enable the accurate
computation of partial derivatives through automatic differentiation. This
finds applications in gradient-based optimization, orbit determination, state
transition matrix computation, sensitivity matrix computation, covariance
propagation, and covariance transformation.

Moreover, we discussed the seamless integration of dSGP4 with state-of-the-art
machine learning backpropagation algorithms. This opens new avenues for
orbital propagation models. We therefore introduce ML-dSGP4, a machine
learning SGP4 model to enhance SGP4 propagations to better match ephemeris
and/or high-precision propagation data. During training, both dSGP4 parameters
and neural network weights can be jointly optimized backpropagating through
both the dSGP4 model and the neural networks, enhancing prediction accuracy to
align more closely with observational data, such as ephemeris data. In this
way, one can maintain the same SGP4 model efficiency and speed, while
improving its accuracy and enabling parallelization over CPU and GPU. We show
how the new propagation model is able to outperform both the SGP4 baseline and
black-box neural network correctors on a sample of around 1,500 Starlink
satellites.

Furthermore, the PyTorch implementation of dSGP4 facilitates straightforward
parallelization on both CPU and GPUs, thereby allowing the processing of
batches of TLE data at the same time. This acceleration holds promise for
applications that display bottlenecks in terms of computational costs, such as
many-to-many conjunction screening in an increasingly congested low-Earth
orbit.

Throughout our paper, we present several experiments that delve into each of
these applications, providing quantitative results that underscore the
practical significance and versatility of dSGP4 in advancing space situational
awareness and orbital prediction. Furthermore, we also release the dSGP4 code
and associated tests and tutorials (which also go through the experiments in
this paper) as an open-source Github repository under the ESA organization.

## References

  * [1] D. A. Vallado, Fundamentals of astrodynamics and applications, Vol. 12, Springer Science & Business Media, 2001.
  * [2] D. Brouwer, Solution of the problem of artificial satellite theory without drag, Astronomical Journal 64 (1959) 378.  doi:10.1086/107958.
  * [3] R. Lyddane, Small eccentricities or inclinations in the brouwer theory of the artificial satellite, The Astronomical Journal 68 (1963) 555.
  * [4] O. Montenbruck, E. Gill, F. Lutze, Satellite orbits: models, methods, and applications, Appl. Mech. Rev. 55 (2) (2002) B27–B28.
  * [5] J. Liu, R. Alford, Semianalytic theory for a close-earth artificial satellite, Journal of Guidance and Control 3 (4) (1980) 304–311.  doi:10.2514/3.55994.
  * [6] I. Pérez, J. F. San-Juan, M. San-Martín, L. M. López-Ochoa, Application of computational intelligence in order to develop hybrid orbit propagation methods, Mathematical Problems in Engineering 2013.  doi:10.1155/2013/631628.
  * [7] J. F. San-Juan, I. Pérez, E. Vergara, M. San Martın, R. López, A. Wittig, D. Izzo, Hybrid SGP4 propagator based on machine-learning techniques applied to GALILEO-type orbits, in: 69th International Astronautical Congress, Bremen, Germany, 2018, pp. 1–11.
  * [8] Y. Kozai, The motion of a close earth satellite, The Astronomical Journal 64 (1959) 367.  doi:10.1086/107957.
  * [9] M. H. Lane, The development of an artificial satellite theory using a power-law atmospheric density representation, in: 2nd Aerospace Sciences Meeting, 1965, p. 35.  doi:10.2514/6.1965-35.
  * [10] D. Vallado, P. Crawford, R. Hujsak, T. Kelso, Revisiting spacetrack report# 3, in: AIAA/AAS Astrodynamics Specialist Conference and Exhibit, 2006, p. 6753.  doi:10.2514/6.2006-6753.
  * [11] F. R. Hoots, R. L. Roehrich, Spacetrack report number 3: Models for propagation of NORAD element sets (1980).  doi:10.21236/ADA093554.
  * [12] K. Cranford, An improved analytical drag theory for the artificial satellite problem, in: Astrodynamics Conference, 1969, p. 925.  doi:10.2514/6.1969-925.
  * [13] M. H. Lane, F. R. Hoots, General perturbations theories derived from the 1965 lane drag theory, Tech. rep., Aerospace Defense Command Peterson AFB CO Office of Astrodynamics (1979).  doi:10.21236/ada081264.
  * [14] T. Kelso, Validation of SGP4 and IS-GPS-200D against GPS precision ephemerides, in: 17th AAS/AIAA Space Flight Mechanics Conference, Sedona, AZ, Vol. 28, 2007, pp. 1,15.
  * [15] T. Flohrer, H. Krag, H. Klinkrad, Assessment and categorization of TLE orbit errors for the us ssn catalogue, risk 8 (9) (2008) 10–11.
  * [16] W. Dong, Z. Chang-yin, An accuracy analysis of the SGP4/SDP4 model, Chinese Astronomy and Astrophysics 34 (1) (2010) 69–76.  doi:10.1016/j.chinastron.2009.12.009.
  * [17] D. A. Vallado, P. J. Cefola, Two-line element sets-practice and use, in: 63rd International Astronautical Congress, Naples, Italy, 2012, pp. 1–14.
  * [18] H. Peng, X. Bai, Machine learning approach to improve satellite orbit prediction accuracy using publicly available data, The Journal of the astronautical sciences 67 (2) (2020) 762–793.  doi:10.1007/s40295-019-00158-3.
  * [19] H. Peng, X. Bai, A medium-scale study of using machine learning fusion to improve TLE prediction precision without external information, Acta Astronautica 204 (2023) 477–491.  doi:10.1016/j.actaastro.2022.06.023.
  * [20] H. Peng, X. Bai, Fusion of a machine learning approach and classical orbit predictions, Acta astronautica 184 (2021) 222–240.  doi:10.1016/j.actaastro.2021.04.017.
  * [21] D. Vallado, P. Crawford, SGP4 orbit determination, in: AIAA/AAS Astrodynamics Specialist Conference and Exhibit, 2008, p. 6770.  doi:10.2514/6.2008-6770.
  * [22] J. Degrave, M. Hermans, J. Dambre, et al., A differentiable physics engine for deep learning in robotics, Frontiers in neurorobotics (2019) 6doi:10.3389/fnbot.2019.00006.
  * [23] T.-M. Li, M. Aittala, F. Durand, J. Lehtinen, Differentiable Monte Carlo ray tracing through edge sampling, ACM Transactions on Graphics (TOG) 37 (6) (2018) 1–11.  doi:10.1145/3272127.3275109.
  * [24] Y. Hu, L. Anderson, T.-M. Li, Q. Sun, N. Carr, J. Ragan-Kelley, F. Durand, Difftaichi: Differentiable programming for physical simulation (2020).  arXiv:1910.00935.
  * [25] D. Izzo, F. Biscani, A. Mereta, Differentiable genetic programming, in: Genetic Programming: 20th European Conference, EuroGP 2017, Amsterdam, The Netherlands, April 19-21, 2017, Proceedings 20, Springer, 2017, pp. 35–51.
  * [26] O. Montenbruck, E. Gill, Real-time estimation of sgp4 orbital elements from gps navigation data, in: International Symposium Space Flight Dynamics, Biarritz, France, 2000, pp. 26–30.
  * [27] D. A. Vallado, B. B. Virgili, T. Flohrer, Improved SSA through orbit determination of two-line element sets, in: ESA Space Debris Conference, 2013, pp. 1–7.
  * [28] J. F. San-Juan, I. Pérez, M. San-Martín, E. P. Vergara, Hybrid SGP4 orbit propagator, Acta Astronautica 137 (2017) 254–260.  doi:10.1016/j.actaastro.2017.04.015.
  * [29] A. Griewank, A. Walther, Evaluating derivatives: principles and techniques of algorithmic differentiation, SIAM, 2008.  doi:10.1137/1.9780898717761.
  * [30] G. F. Corliss, Applications of differentiation arithmetic, in: Reliability in computing, Elsevier, 1988, pp. 127–148.  doi:10.5555/60181.60188.
  * [31] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, J. M. Siskind, Automatic differentiation in machine learning: a survey, Journal of Marchine Learning Research 18 (2018) 1–43.  doi:10.5555/3122009.3242010.
  * [32] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, et al., Pytorch: An imperative style, high-performance deep learning library, Advances in neural information processing systems 32.  doi:10.5555/3454287.3455008.
  * [33] S. Ruder, An overview of gradient descent optimization algorithms (2017).  arXiv:1609.04747.


# Ocean-DC: An analysis ready data cube framework for environmental and
climate change monitoring over the port areas

Ioannis Kavouras
School of Rural, Surveying and Geoinformatics Engineering
National Technical University of Athens
9, Iroon Polytechniou st, 15772, Zografou
<EMAIL_ADDRESS>
Ioannis Rallis
School of Rural, Surveying and Geoinformatics Engineering
National Technical University of Athens
9, Iroon Polytechniou st, 15772, Zografou
<EMAIL_ADDRESS>
Nikolaos Doulamis
School of Rural, Surveying and Geoinformatics Engineering
National Technical University of Athens
9, Iroon Polytechniou st, 15772, Zografou
<EMAIL_ADDRESS>
Anastasios Doulamis
School of Rural, Surveying and Geoinformatics Engineering
National Technical University of Athens
9, Iroon Polytechniou st, 15772, Zografou
<EMAIL_ADDRESS>

###### Abstract

The environmental hazards and climate change effects causes serious problems
in land and coastal areas. A solution to this problem can be the periodic
monitoring over critical areas, like coastal region with heavy industrial
activity (i.e., ship-buildings) or areas where a disaster (i.e., oil-spill)
has occurred. Today there are several Earth and non-Earth Observation data
available from several data providers. These data are huge in size and usually
it is needed to combine several data from multiple sources (i.e., data with
format differences) for a more effective evaluation. For addressing these
issues, this work proposes the Ocean-DC framework as a solution in data
harmonization and homogenization. A strong advantage of this Data Cube
implementation is the generation of a single NetCDF product that contains
Earth Observation data of several data types (i.e., Landsat-8 and Sentinel-2).
To evaluate the effectiveness and efficiency of the Ocean-DC implementation,
it is examined a case study of an oil-spill in Saronic gulf in September of
2017. The generated 4D Data Cube considers both Landsat-8,9 and Sentinel-2
products for a time-series analysis, before, during, and after the oil-spill
event. The Ocean-DC framework successfully generated a NetCDF product,
containing all the necessary remote sensing products for monitoring the oil-
spill disaster in the Saronic gulf.

_Keywords_ Data Cubes, Climate Change, Coastal Areas

##  1 Introduction

The modern urban environments have been developed with the purpose of
interconnecting people, infrastructures, activities and several other
resources, making them susceptible to environmental disasters and climate
change effects Mi et al. (2019); Leal Filho et al. (2019); Olabi and
Abdelkareem (2022). Daily industrial activities, especially near the port
areas, have a negative impact on the local environment and climate endangering
both the human health and well-being, as well as the local nature species.

All the human activities result in the following issues: (a) lack of urban
greenness Temenos et al. (2022); (b) flooding Pour et al. (2020); (c) non-
healthy atmospheric quality Ravindra et al. (2019); (d) increment of
greenhouse gases emissions Shen et al. (2020); and (e) the phenomenon of heat
island effect Parker (2010). In most of these cases, the impact can be
evaluated by continuously monitoring the activities over the area of interest
using Earth Observation (EO) data provided by space-borne technologies (i.e.,
satellites) and aerial or field climate measurements. In other cases, the
impact can be instantaneous, causing a huge disaster over a radius, like in an
oil-spill situation.

Today, there is available a large number of EO data and it expands each week
Gomes et al. (2020). Moreover, the available data are produced from different
data providers (i.e., Copernicus, USGS, ESA, NASA, etc.), thus there are
differences in distribution format, spatial resolution, data collection, etc.
In some cases, however, two missions from different provides (i.e., Landsat-8
and Sentinel-2) have similarities, thus such data can be used in combination
for thickening a time-series dataset. For a successful data harmonization and
homogenization a Data Cube (DC) implementation can be used Kopp et al. (2019).

Figure 1: An example of a Data Cube’s dimension expanding.

A DC Baumann et al. (2019); Harinarayan et al. (1996) is a multidimensional
array containing all the necessary information (i.e., data, metadata,
coordinates, etc.) describing a phenomenon. To better understand a basic Earth
Observation DC framework, let us denote as 2D array ($B_{i}(x,y)$)
corresponding an spectral band (i.e., image) of a satellite product. Thus,
each coordinate set (x, y) corresponds to a radiometric reflectance value for
this band. However, each satellite retrieves data on several wavelengths,
which means that they generate several images for different wavelengths of the
spectrum (i.e., Coastal Aerosol, Blue, Green, Red, Near Infrared, etc.). The
combination of these images can generate several remote sensing products Zhang
et al. (2010); Wu and Li (2009) for indicating a phenomenon. Usually, the raw
data are captured in different spatial resolutions (i.e., 10m/px, 15m/px,
30m/px, 60m/px, etc), thus any combinations and/or calculations can not be
achieved in the original sate of data. For addressing this issue, let us now
denote a 3D array ($DC_{i}(b,x,y)$) for storing all the spectral bands $N$
($DC_{i}(b,x,y)=[B_{1}(x,y),B_{2}(x,y),...,B_{N}(x,y)]$) in a common spatial
resolution format and coordination reference system. Finally, the monitoring
of a phenomenon introduces the parameter of time ($t$). Thus, for effectively
monitoring a periodic phenomenon (i.e., climate) over an are and detecting
abnormalities it is necessary to generate a DC containing data from various
sources. Thus, a final DC implementation can be a 4D-DC (Hypercube) defined as
$HDC_{i}(t,b,x,y)=[DC_{1}(b,x,y),DC_{2}(b,x,y),...DC_{N}(b,x,y)]$, containing
$N$ 3D-DC of different time periods (i.e., before and after an event). Figure
1 illustrates the DC dimension expanding as described in this paragraph.

For the successful implementation of a 4D-DC it is necessary to harmonize and
homogenize all of the available 3D-DCs. If the original data are provided from
the same data provider, then the implementation can follow the same
processing, however additional processes are necessary if the data are
provided from more than one sources. For addressing this challenge, in this
work, we propose the Oceanids Data Cubes (Ocean-DCs) as a novel approach for
continuously and quickly monitoring an area of interest using several kinds of
data from the same or different data sources. Moreover, the Ocean-DC provides
data harmonization and homogenization by applying the common template for
Landsat-8 and Sentinel-2 data products, as well as it offers a complete ready
to use analysis products by calculating the most common used remote sensing
indices. To be more precise, Ocean-DCs’ products can be used as is in a
decision support system platform alongside in synergy with other data of the
platform for indicating environmental and climate change hazards in many
applications, emphasizing to port, marine and coastal applications (i.e.,
monitoring, forecasting, etc.).

For demonstrating the Ocean-DC implementation a monitoring case study is
provided in this manuscript. For the case study it is selected the oil-spill
disaster that took place on the 10th of September 2017 in the Saronikos gulf
(Greece) causing huge environmental disaster to the surrounded shores and
sea/water area. This case study aims to provide an application example of the
proposed Ocean-DC for periodically monitoring a phenomenon based on time-
series analysis. However, Ocean-DC products can be used differently depended
on the application and the related phenomenon.

The rest of this paper is organized as follow: (a) Section 2 presents the
related literature review; (b) Section 3 describes the proposed Ocean-DC
framework; (c) Section 4 provides a brief demonstration of environmental
monitoring based on an Ocean-DC product; (d) Section 5 discuss the outcomes of
this work and (e) Section 6 summarizes this work.

##  2 Related Work

A Data Cube is a multidimensional array capable or a data structure that
allows the efficient analysis of big data and complex data structures along
the multiple dimensions. A DC framework is mostly used for applications
regarding big data storage and online analytical processing (OLAP) Sarawagi et
al. (1998) based on complex queries Gray et al. (1997) and analysis Francia et
al. (2022). Major functionalities of a DC framework are: (a) multidimensional
representation Leprince et al. (2021); (b) cube structure for big data
visualizations Song et al. (2018); (c) data aggregation Giuliani et al.
(2020a); (d) fast query performance Miranda et al. (2017); (e) OLAP operations
Tardío et al. (2020); and (d) smart big data applications Majeed et al.
(2021).

The current DC’s implementations are emphasizing on organizing the data on
multiple dimensions Han et al. (2005) like time, geography, product, or
customer. Each dimension represents a different aspect of the data. Thus, DCs
structure is often visualized as cube with multiple dimensions. Each cell in
the cube’s representation is a unique combination of dimensional values. For
example, in a remote sensing DC an indicative combination can be the (time x
spectral bands x geographical coordinates). Moreover, in other DC examples,
like climatic and geographical DC, the data can be aggregated using
traditional aggregation functions like mean, count, average, minimum, and
maximum.

However, it is important to note that the numerical values alone are not
enough to describe a phenomenon. Another strong characteristic of the DC’s
architecture is the inclusion of metadata information that describes the
numerical values of the data Boulil et al. (2015). Such information includes
the metric system of the data, the origin of the data, geo-spatial
information, etc. This information can be used for achieving fast query
performance, which is one of the main advantages of the utilization of DCs,
especially when dealing with large volumes of data (i.e., addressing big data
storage performance issues).

Several Data Cube implementations have been suggested the recent years. An
indicative example can be the work of Ferreira et al. Ferreira et al. (2020),
where they proposed an effectual way for generating Earth Observation Data
Cubes (EODCs) over the Brazil region. They additionally are exploiting the
benefits of the Open Data Cubes (ODC) platform, contributing by conducting a
analysis-ready data and multidimensional data cubes from remote sensing images
to effectively mapping the land use and land cover based on employing remote
sensing imagery time series analysis and machine learning techniques.

The Australian Geoscience Data Cubes (ADGCs) is one of the first data cubes
implementation, proposed by Lewis et al. Lewis et al. (2017) and manages the
challenges of the Big Data, such us the storage volume, velocity, and
polymorphism or variety, limiting the earth observation data usefulness. Some
of the application, where the ADGCs have been applied are: (a) water
observations; (b) coastal erosion; (c) agriculture monitoring; (d) forest
cover changes; and (e) biodiversity.

In the work of Poussin et al. Poussin et al. (2019) the advantages and
limitations of the Earth Observation Swiss Data Cubes (EOSDCs) are exploited
by monitoring the snow cover change in Gran Paradiso national Park, Alps.
EOSDCs architecture is based on the Open Data Cubes software Chatenoux et al.
(2021), which is a geospatial data management and analysis software based on
open-source technologies for processing remote sensing data by providing a
framework for accessing, storing, managing, and analyzing grid-aligned
satellite Earth Observation data collections in huge quantities.

Other works utilizing the EOSDCs implementations are proposed by Giuliani et
al. Giuliani et al. (2020a, b). In the first work Giuliani et al. (2020a) they
propose a methodology for generating accurate and consistent land degradation
products based on remote sensing imagery datasets. They uses a data cube and
several sun-indicators for assessing land productivity and cover, as well as
soil organic carbon changes for evaluating land degradation. The second work
Giuliani et al. (2020b) suggests the usage of a data cube on demand (DCoD),
which is a script-based tool, for automating the generation of data cubes
under predetermined by the user requirements. Such requirements can be the
area of interest, a time period, the type of the sensor, or the satellite
product family. Their implementation has been successfully applied to two
regions (Bolivia and Democratic Republic of Congo) for environmental
monitoring.

Finally, Temenos et al. Temenos et al. (2024) propose a context-aware adaptive
data cube (C2A-DC) framework, which utilize Earth Observation data for
environmental monitoring and mitigating the climate change effects. Their
proposed methodology combines the data cubes formation with calculation of
remote sensing operations, including deep and machine learning algorithms for
data classification and harmonization. Thus, their outcome is a context-aware
adaptive framework, which successfully allows the environmental monitoring and
climate change effects based on user predetermined preferences.

###  2.1 Contribution

Considering the current literature is is noted that several DC implementations
have been proposed for different applications over specific regions. However,
these implementations are focusing only on the harmonization and
homogenization of data provided from the same source, even if they provide DCs
of several data sources. As an example, the current DC implementations can
provide different data like Landsat-8, Sentinel-2, Sentinel-3, etc., however
each cube contains data only from one source. In addition, the current DC
implementations provide additional products (i.e., NDVI, NDWI, etc.) on demand
based on some query selection.

This work approach the problem of handling the big data from a different
scope. The proposed Ocean-DC implementation provide data harmonization and
homogenization in a 4D-DC scenario, including data from several sources with
commonalities in the same cube. For example, in the case study of this
manuscript are used Landsat-8 and Sentinel-2 data for the time series
analysis, that are included in the same 4D-DC. In addition, the Ocean-DC pre-
calculates the most common remote sensing products (e.g., NDVI, NDWI, Land
Surface Temperature, etc.) providing to the end-users a ready to use product
for addressing several environmental and climate change challenges.

##  3 Proposed Ocean-DC framework

The main objective of the Ocean-DC implementation is to provide harmonization
and homogenization between the several available remote sensing and climate
data products (i.e., Earth Observation [EO] products), that can be retrieved
from multiple data providers (i.e., USGS, Copernicus, NASA). In addition, non-
EO data can be used, like the boundaries of the area of interest in the
processes. All these data are collected in the first step of the whole
framework architecture, as illustrated in Figure 2.

Figure 2: The workflow architecture of the Ocean-DC implementation.

To be more precise, the Data Collection considers several EO and non-EO data
including geolocated raster remote sensing and climate products and the non-EO
products are usually vector type. Currently, the Ocean-DC implementation
supports products from USGS USGS (2024) (i.e., Landsat-8,9), Copernicus
Services Agency (2024) (i.e., Sentinel-2,3,5P and CERRA climate data) and NASA
NASA (2024) (i.e., MODIS) EO data. The supported non-EO data are polygon
shapefile (SHP) geometries. These data are fed to the next step (i.e., Data
Harmonization/Homogenization) for generating the Ocean-DCs.

The Data Harmonization/Homogenization step is the core tool Oce (2024) of the
workflow architecture. In this step, the EO products are reprojected to a
common Coordinate Reference System (CRS) to get the same coordinates per
pixel. Usually, this coordinate system is the same of the SHP geometry, but
the user can select any CRS or EPSG. This step is necessary because the EO
data are in different CRS. For instance Landsat-8 and Sentinel-2 raw data are
provided in different CRS. Thus, this step is necessary before applying any
other process.

After the data reprojection, the framework clips the products to a box
geometry based on the user input (i.e., SHP geometry). The box geometry is
calculated by the minimum and maximum (x,y) coordinates of this SHP geometry.
Image clipping is a critical step, because it significantly minimizes the size
of the original images and as a natural result the size of the final Ocean-DC
product. After successfully clipping the raster EO products, they are in the
same coordinate reference system and visualize the area of interest.

The next following two steps are the image resizing/resampling and the
calculation of the most used remote sensing product, in the case of Landsat-8
and Sentinel-2 EO products. The image resizing/resampling process corrects the
spatial resolution of the images and ensures that all images have the same
width and height pixels. Note that it is not added new information to the
original data. This means that an image with spatial resolution 60m/px will be
resampled to 10m/px by splitting each pixel to a 6x6 grid with the same value
(i.e., pseudo-resolution). This step increases the size of the final product,
but it is very important, because the data cube cannot be constructed with
different image resolutions. In addition, it is impossible to perform raster
calculation with images in different resolutions. Thus, the image
resizing/resampling is a mandatory step to proceed in the next steps.

Finally, the calculation of the most common remote sensing products is
performed when it is necessary. To be more precise, the remote sensing indices
are calculated when specific band information is known (i.e., optical bands,
near-infrared, SWIR, etc.). Currently, this step is applied for Landsat-8,9
and Sentinel-2 products. Table 1 summarizes the satellite calculated products,
which are containing in an Ocean-DC implementation. The processing of the
other products is terminated in the image resizing/resampling step.

The implementation of the Ocean-DC is performed using Python and the state-of-
the-art and open-source DC handling library Xarray Xarray (2024). This library
is empowered with the ready to use methods of rasterio Ras (2024), which is
another python library for handling EO data. The final outcome of the data
harmonization/homogenization framework is a NetCDF file, which is a state-of-
the-art format for geodata and is supported by several GIS software like QGIS
Community (2024).

Finally, the Ocean-DC products can be combined with other non-EO data like
images or drone footage to produce advanced cartographic products or monitor
several phenomena. Some examples of post-processing applications can be the
monitoring of a disaster (i.e., flooding, forest fire, oil-spill, etc.) or the
generation of thematic maps that describes a phenomenon. Other application
could be the monitoring of air quality or climate factors (i.e., temperature,
precipitation).

Table 1: Summary of the products, where containing in an Ocean-DC file.  BAND ID | Description | BAND ID | Description
---|---|---|---
1 | COASTAL AEROSOL | 23 | LSWI-1
2 | BLUE | 24 | LSWI-2
3 | GREEN | 25 | ARVI
4 | RED | 26 | MSAVI2
5 | NIR | 27 | MTVI2
6 | SWIR-1 | 28 | VARI
7 | SWIR-2 | 29 | TGI
8 | CIRRUS | 30 | LST-1
9 | VRE-1 | 31 | LST-2
10 | VRE-2 | 32 | LST-CELSIUS
11 | VRE-3 | 33 | LST-FHRENHEIT
12 | VRE-4 | 34 | VCI
13 | WATER VAROUR | 35 | MNDWI-1
14 | PANCHROMATICS | 36 | MNDW-2
15 | TIRS-1 | 37 | WRI-1
16 | TIRS-2 | 38 | WRI-2
17 | NDVI | 39 | NDTI
18 | NDWI | 40 | AWEI
19 | NDVI (VRE-1) | 41 | OSI
20 | NDVI (VRE-2) | 42 | NBR-1
21 | NDVI (VRE-3) | 43 | NBR-2
22 | NDVI (VRE-4) |  |

##  4 Case Study Analysis of the oil-spill in Saronikos gulf at Septempber
2017

For the demonstration of the Ocean-DC framework is used the oil-spill disaster
that took place in Saronic gulf in 10th of September 2017. This disaster has
been also investigated in the work of Kolokoussis and Karathanassi Kolokoussis
and Karathanassi (2018) where they use the NDWI among other indices to
visually enhance the critical areas of the oil-spill. In a similar way, this
work investigates the oil-spill disaster emphasizing mostly on the
possibilities that Ocean-DC framework provides to the end-users. For this
reason, in this analysis it is used a dataset composed of Landsat-8 and
Sentinel-2 products in a time-period before the oil-spill (2016 to 9th of
September 2017), some days and weeks after the the oil-spill (10th of
September to October 2017) and recent images (2023-2024). Thus, the
implemented Ocean-DC product is consisted of 16 satellite imagery products (5
Landsat-8 and 11 Sentinel-2). For visualizing the oil-spill effectively it
used the NDWI, WRI-2, and OSI Zakzouk et al. (2024) indices for monitoring the
water quality of Saronic gulf during the aforementioned periods.

The results of the generated Ocean-DC product are illustrated in Figure 3. The
visualization of the generated product performed inside QGIS software and
considers the natural (RGB) color composition, the NDWI, WRI, and OSI indices
illustrated using pseudo-palette colors. The NDWI product is calculated using
the Equation 1 and the output’s image values are in range $[-1.0,1.0]$ as
follows: $[-1.0,-0.3)\rightarrow$ Drought, Non-Aqueous Surfaces;
$[-0.3,0.0)\rightarrow$ Moderate drought, non-aqueous surfaces;
$[0.0,0.2)\rightarrow$ Flooding, humidity; $[0.2,1.0]\rightarrow$ Water
Surface NH (2022). In most of cases, NDWI values in range $[0.0,0.2)$,
indicating the humidity over the Saronicos gulf. Near the coastal areas the
values are higher. The only exception is observed in the image of 13th of
September, 2017, where it is observed the abnormality in of the oil-spill in
the middle of the Saronic Gulf.

$NDWI=\frac{NIR-SWIR-1}{NIR+SWIR-1}$ (1)

The WRI index is calculated by the Equation 2 and indicates the water ratio in
range $[-1.0,1.0]$ as follows: $[-1.0,-0.75)\rightarrow$ Critical water areas;
$[-0.75,-0.25)\rightarrow$ Normal water quality $[-0.25,0.00)\rightarrow$ Wet
ground / Vegetation; $[-0.25,0.00)\rightarrow$ Wet ground / Vegetation;
$[0.00,1.00]\rightarrow$ Ground / Infrastructures. Based on the observations
of the WRI in Figure 3 is is clearly observed a periodic phenomenon of
critical water areas among the Salamina coast line. Moreover, in the image of
13th of September, 2017 critical values are observed in the area of the oil-
spill (green box in Figure 3), as well as similar areas are indicated in the
following weeks after the oil-spill. Another observation of the analysis, is
that Landsat-8 products were unable to indicate the oil-spill due to low
spatial resolution (i.e., 30m/px).

$WRI-2=\frac{GREEN-RED}{NIR+SWIR-2}$ (2)

The OSI (Oil-Spill Index) is referred in the work of Zakzouk et al. (2024) and
is calculated by the Equation 3. A first observation is the difference of the
products’ values between Landsat-8,9 and Sentinel-2 data. In Landsat-8,9
products the index effectively distinguish the area of interest into water
bodies ($values>2.5$), vegetation regions ($1.9<values<2.5$), and
building/infrastructures ($values<1.9$). In the Sentinel-2 products can be
divided to the following ranges: $[0.75,1.00)\rightarrow$ vegetation regions;
$[1.00,1.9)\rightarrow$ building/infrastructures; and $values>1.9\rightarrow$
water bodies. In the water bodies are additionally observed two trends: (a)
$[1.9,2.5)\rightarrow$ indicate areas emitting strongly in blue spectral
radiance; (b) $>2.5\rightarrow$ indicate areas where water absorbs fully the
blue spectral radiance. Considering the aforementioned observations in the
product corresponding to the 13th of September, 2017, oil-spill patterns are
indicated.

$OSI=\frac{RED+GREEN}{BLUE}$ (3)

Figure 3: Time series monitoring of the water quality of Saronikos gulf before
and after the oil-spill using the RGB, NDWI, WRI, and OSI products of the
Ocean-DC implementation file. The green box indicates the area of the oil-
spill.

Combining all the information that is observed by all these indexes, the
critical area is bounding by the green box as illustrated in Figure 3. In
addition, it is necessary to mention that this time-series analysis indicates
other information as well. For instance, among Salamina’s coast line is is
observed a periodic and recursive pattern of dark pixel values in NDWI and WDI
products. This phenomenon indicates additional areas of interest that may be
useful for specific use cases related to climate change. Moreover, due to the
repeatability of this phenomenon, it is safely to assume that is not
correlated with the oil-spill disaster, because it is also observed in images
before the oil-spill.

##  5 Discussion

The main objective of this work is to introduce a new Data Cube implementation
for coastal region monitoring named Ocean-DC. The Ocean-DC framework handles
several kinds of data combining both EO and non-EO products for generating a
NetCDF file that contains all the necessary information over a predetermined
area of interest. In addition, during the cube’s generation it calculates the
most common remote sensing products to help more the end-users by providing
them with a ready to use product. As an example, in this paper it was
generated a 4D Ocean-DC for monitoring the oil-spill disaster that took place
at 10th of September, 2017. The generated Ocean-DC product successfully
applied data from different sources (i.e., Landsat-8,9 and Sentinel-2) to a
common template and provided a strong ready to use dataset for time series
analysis. Thus, the advantages of the Ocean-DC framework are the following:

  1. 1.

It utilizes several EO and non-EO data for generating a ready to use product.

  2. 2.

It supports several EO data from different sources and harmonize/homogenize
them to a common. This feature significantly accelerates the post-processing
analysis by allowing the end-users to visualize several data using a single
file.

  3. 3.

It calculates the state-of-the-art remote sensing data providing a ready to
use product. This feature also accelerates the post-processing analysis by
allowing the end-users to emphasize only on the interpretation of the
generated products without spending time to calculate the products.

  4. 4.

The aforementioned characteristics allowing the usage of the Ocean-DC products
by both experts and non-experts, because after the generation of the DC, the
end-user have to visualize and interpreter the data.

  5. 5.

The outcome of the Ocean-DC framework is a NetCDF file, which is readable by
the most common GIS software. Thus, it is supported by both free and open-
source and commercial GIS software, allowing to end-users to use their
preferred software.

Summarizing the discussion, it is worth mentioning that the Ocean-DC framework
is not finalized yet. This means, that the framework is updated frequently by
expanding this framework with additional features. Moreover, further testing
is necessary for indicating further the efficiency of the Ocean-DC framework.
Thus, possible future works will combine the Ocean-DC framework with Machine
Learning and Deep Learning algorithms for generating more products like cloud
and shadow masks, expanding more the applicability of the Ocean-DC products
for coastal monitoring. Moreover, other application like in-land vegetation,
deforestation, forest fires, or flooding monitoring will be examined as well.

##  6 Conclusions

In conclusion, this work introduces a new DC framework that handles several EO
and non-EO data from different sources for periodic monitoring the
environmental hazards and climate change effects over the coastal areas. The
framework may still be under development, however, it can be effectively
utilized for time-series analysis and the visualization of disasters over the
coastal areas, as indicated in the presented case study. As a result, the
Ocean-DC framework can effectively utilized in monitoring phenomena
applications. However, the Ocean-DC as an emerging and currently in
development implementation will be updated frequently to support more
applications by utilizing novel and smart technologies.

##  7 Acknowledgments

This work is supported by the European Union Funded project OCEANIDS ”User-
driven applications and tools for Climate-Informed Maritime Spatial Planning
and integrated seascape management, towards a resilient & inclusive Blue
Economy”, under the EU HE research and innovation program under GA No.
101112919.

## References

  * Mi et al. [2019] Zhifu Mi, Dabo Guan, Zhu Liu, Jingru Liu, Vincent Viguié, Neil Fromer, and Yutao Wang.  Cities: The core of climate change mitigation.  _Journal of Cleaner Production_ , 207:582–589, 2019.
  * Leal Filho et al. [2019] Walter Leal Filho, Abdul-Lateef Balogun, Olawale Emmanuel Olayide, Ulisses M Azeiteiro, Desalegn Y Ayal, Pastor David Chavez Muñoz, Gustavo J Nagy, Paulette Bynoe, Otienoh Oguge, N Yannick Toamukum, et al.  Assessing the impacts of climate change in cities and their adaptive capacity: Towards transformative approaches to climate change adaptation and poverty reduction in urban areas in a set of developing countries.  _Science of the Total Environment_ , 692:1175–1190, 2019.
  * Olabi and Abdelkareem [2022] AG Olabi and Mohammad Ali Abdelkareem.  Renewable energy and climate change.  _Renewable and Sustainable Energy Reviews_ , 158:112111, 2022.
  * Temenos et al. [2022] Anastasios Temenos, Ioannis N Tzortzis, Maria Kaselimi, Ioannis Rallis, Anastasios Doulamis, and Nikolaos Doulamis.  Novel insights in spatial epidemiology utilizing explainable ai (xai) and remote sensing.  _Remote Sensing_ , 14(13):3074, 2022.
  * Pour et al. [2020] Sahar Hadi Pour, Ahmad Khairi Abd Wahab, Shamsuddin Shahid, Md Asaduzzaman, and Ashraf Dewan.  Low impact development techniques to mitigate the impacts of climate-change-induced urban floods: Current trends, issues and challenges.  _Sustainable Cities and Society_ , 62:102373, 2020.
  * Ravindra et al. [2019] Khaiwal Ravindra, Preety Rattan, Suman Mor, and Ashutosh Nath Aggarwal.  Generalized additive models: Building evidence of air pollution, climate change and human health.  _Environment international_ , 132:104987, 2019.
  * Shen et al. [2020] Maocai Shen, Wei Huang, Ming Chen, Biao Song, Guangming Zeng, and Yaxin Zhang.  (micro) plastic crisis: un-ignorable contribution to global greenhouse gas emissions and climate change.  _Journal of Cleaner Production_ , 254:120138, 2020.
  * Parker [2010] David E Parker.  Urban heat island effects on estimates of observed climate change.  _Wiley Interdisciplinary Reviews: Climate Change_ , 1(1):123–133, 2010.
  * Gomes et al. [2020] Vitor CF Gomes, Gilberto R Queiroz, and Karine R Ferreira.  An overview of platforms for big earth observation data management and analysis.  _Remote Sensing_ , 12(8):1253, 2020.
  * Kopp et al. [2019] Steve Kopp, Peter Becker, Abhijit Doshi, Dawn J Wright, Kaixi Zhang, and Hong Xu.  Achieving the full vision of earth observation data cubes.  _Data_ , 4(3):94, 2019.
  * Baumann et al. [2019] Peter Baumann, Dimitar Misev, Vlad Merticariu, and Bang Pham Huu.  Datacubes: Towards space/time analysis-ready data.  _Service-Oriented Mapping: Changing Paradigm in Map Production and Geoinformation Management_ , pages 269–299, 2019.
  * Harinarayan et al. [1996] Venky Harinarayan, Anand Rajaraman, and Jeffrey D Ullman.  Implementing data cubes efficiently.  _Acm Sigmod Record_ , 25(2):205–216, 1996.
  * Zhang et al. [2010] RenHua Zhang, Jing Tian, ZhaoLiang Li, HongBo Su, ShaoHui Chen, and XinZhai Tang.  Principles and methods for the validation of quantitative remote sensing products.  _Science China Earth Sciences_ , 53:741–751, 2010.
  * Wu and Li [2009] Hua Wu and Zhao-Liang Li.  Scale issues in remote sensing: A review on analysis, processing and modeling.  _Sensors_ , 9(3):1768–1793, 2009.
  * Sarawagi et al. [1998] Sunita Sarawagi, Rakesh Agrawal, and Nimrod Megiddo.  Discovery-driven exploration of olap data cubes.  In _Advances in Database Technology—EDBT’98: 6th International Conference on Extending Database Technology Valencia, Spain, March 23–27, 1998 Proceedings 6_ , pages 168–182. Springer, 1998.
  * Gray et al. [1997] Jim Gray, Surajit Chaudhuri, Adam Bosworth, Andrew Layman, Don Reichart, Murali Venkatrao, Frank Pellow, and Hamid Pirahesh.  Data cube: A relational aggregation operator generalizing group-by, cross-tab, and sub-totals.  _Data mining and knowledge discovery_ , 1:29–53, 1997.
  * Francia et al. [2022] Matteo Francia, Patrick Marcel, Verónika Peralta, and Stefano Rizzi.  Enhancing cubes with models to describe multidimensional data.  _Information Systems Frontiers_ , 24(1):31–48, 2022.
  * Leprince et al. [2021] Julien Leprince, Clayton Miller, and Wim Zeiler.  Data mining cubes for buildings, a generic framework for multidimensional analytics of building performance data.  _Energy and Buildings_ , 248:111195, 2021.
  * Song et al. [2018] Ying Song, Yingling Fan, Xin Li, and Yanjie Ji.  Multidimensional visualization of transit smartcard data using space–time plots and data cubes.  _Transportation_ , 45:311–333, 2018.
  * Giuliani et al. [2020a] Gregory Giuliani, Bruno Chatenoux, Antonio Benvenuti, Pierre Lacroix, Mattia Santoro, and Paolo Mazzetti.  Monitoring land degradation at national level using satellite earth observation time-series data to support sdg15–exploring the potential of data cube.  _Big Earth Data_ , 4(1):3–22, 2020a.
  * Miranda et al. [2017] Fabio Miranda, Lauro Lins, James T Klosowski, and Claudio T Silva.  Topkube: A rank-aware data cube for real-time exploration of spatiotemporal data.  _IEEE Transactions on visualization and computer graphics_ , 24(3):1394–1407, 2017.
  * Tardío et al. [2020] Roberto Tardío, Alejandro Maté, and Juan Trujillo.  A new big data benchmark for olap cube design using data pre-aggregation techniques.  _Applied Sciences_ , 10(23):8674, 2020.
  * Majeed et al. [2021] Arfan Majeed, Yingfeng Zhang, Shan Ren, Jingxiang Lv, Tao Peng, Saad Waqar, and Enhuai Yin.  A big data-driven framework for sustainable and smart additive manufacturing.  _Robotics and Computer-Integrated Manufacturing_ , 67:102026, 2021.
  * Han et al. [2005] Jiawei Han, Yixin Chen, Guozhu Dong, Jian Pei, Benjamin W Wah, Jianyong Wang, and Y Dora Cai.  Stream cube: An architecture for multi-dimensional analysis of data streams.  _Distributed and Parallel Databases_ , 18:173–197, 2005.
  * Boulil et al. [2015] Kamal Boulil, Sandro Bimonte, and Francois Pinet.  Conceptual model for spatial data cubes: A uml profile and its automatic implementation.  _Computer Standards & Interfaces_, 38:113–132, 2015.
  * Ferreira et al. [2020] Karine R. Ferreira, Gilberto R. Queiroz, Lubia Vinhas, Rennan F. B. Marujo, Rolf E. O. Simoes, Michelle C. A. Picoli, Gilberto Camara, Ricardo Cartaxo, Vitor C. F. Gomes, Lorena A. Santos, Alber H. Sanchez, Jeferson S. Arcanjo, José Guilherme Fronza, Carlos Alberto Noronha, Raphael W. Costa, Matheus C. Zaglia, Fabiana Zioti, Thales S. Korting, Anderson R. Soares, Michel E. D. Chaves, and Leila M. G. Fonseca.  Earth observation data cubes for brazil: Requirements, methodology and products.  _Remote Sensing_ , 12(24), 2020.  ISSN 2072-4292.  doi:10.3390/rs12244033.  URL https://www.mdpi.com/2072-4292/12/24/4033.
  * Lewis et al. [2017] Adam Lewis, Simon Oliver, Leo Lymburner, Ben Evans, Lesley Wyborn, Norman Mueller, Gregory Raevksi, Jeremy Hooke, Rob Woodcock, Joshua Sixsmith, Wenjun Wu, Peter Tan, Fuqin Li, Brian Killough, Stuart Minchin, Dale Roberts, Damien Ayers, Biswajit Bala, John Dwyer, Arnold Dekker, Trevor Dhu, Andrew Hicks, Alex Ip, Matt Purss, Clare Richards, Stephen Sagar, Claire Trenham, Peter Wang, and Lan-Wei Wang.  The australian geoscience data cube — foundations and lessons learned.  _Remote Sensing of Environment_ , 202:276–292, 2017.  ISSN 0034-4257.  doi:https://doi.org/10.1016/j.rse.2017.03.015.  URL https://www.sciencedirect.com/science/article/pii/S0034425717301086.
  * Poussin et al. [2019] Charlotte Poussin, Yaniss Guigoz, Elisa Palazzi, Silvia Terzago, Bruno Chatenoux, and Gregory Giuliani.  Snow cover evolution in the gran paradiso national park, italian alps, using the earth observation data cube.  _Data_ , 4(4), 2019.  ISSN 2306-5729.  doi:10.3390/data4040138.  URL https://www.mdpi.com/2306-5729/4/4/138.
  * Chatenoux et al. [2021] Bruno Chatenoux, Jean-Philippe Richard, David Small, Claudia Roeoesli, Vladimir Wingate, Charlotte Poussin, Denisa Rodila, Pascal Peduzzi, Charlotte Steinmeier, Christian Ginzler, et al.  The swiss data cube, analysis ready data archive using earth observations of switzerland.  _Scientific data_ , 8(1):295, 2021.
  * Giuliani et al. [2020b] Gregory Giuliani, Bruno Chatenoux, Thomas Piller, Frédéric Moser, and Pierre Lacroix.  Data cube on demand (dcod): Generating an earth observation data cube anywhere in the world.  _International Journal of Applied Earth Observation and Geoinformation_ , 87:102035, 2020b.
  * Temenos et al. [2024] Anastasios Temenos, Nikos Temenos, Ioannis N. Tzortzis, Ioannis Rallis, Anastasios Doulamis, and Nikolaos Doulamis.  C2a-dc: A context-aware adaptive data cube framework for environmental monitoring and climate change crisis management.  _Remote Sensing Applications: Society and Environment_ , 34:101171, 2024.  ISSN 2352-9385.  doi:https://doi.org/10.1016/j.rsase.2024.101171.  URL https://www.sciencedirect.com/science/article/pii/S2352938524000351.
  * USGS [2024] USGS.  Earth explorer, 2024.  URL https://earthexplorer.usgs.gov/.  (Accessed on 29 Feb 2024).
  * Agency [2024] European Space Agency.  Copernicus hub, 2024.  URL https://dataspace.copernicus.eu/.  (Accessed on 29 Feb 2024).
  * NASA [2024] NASA.  Modis data hub, 2024.  URL https://modis.gsfc.nasa.gov/data/.  (Accessed on 29 Feb 2024).
  * Oce [2024] Ocean-dc framework, 2024.  URL https://gitlab.com/JohnCrabs/oceanides_data_cubes.  (Accessed on 29 Feb 2024).
  * Xarray [2024] Xarray.  Xarray documentation, 2024.  URL https://docs.xarray.dev/en/stable/.  (Accessed on 29 Feb 2024).
  * Ras [2024] Rasterio documentation, 2024.  URL https://rasterio.readthedocs.io/en/stable/.  (Accessed on 29 Feb 2024).
  * Community [2024] QGIS Community.  Quantum gis, 2024.  URL https://qgis.org/en/site/.  (Accessed on 29 Feb 2024).
  * Kolokoussis and Karathanassi [2018] Polychronis Kolokoussis and Vassilia Karathanassi.  Oil spill detection and mapping using sentinel 2 imagery.  _Journal of Marine Science and Engineering_ , 6(1), 2018.  ISSN 2077-1312.  doi:10.3390/jmse6010004.  URL https://www.mdpi.com/2077-1312/6/1/4.
  * Zakzouk et al. [2024] Mohamed Zakzouk, Islam Abou El-Magd, Elham M Ali, Abdulaziz M Abdulaziz, Amjad Rehman, and Tanzila Saba.  Novel oil spill indices for sentinel-2 imagery: A case study of natural seepage in qaruh island, kuwait.  _MethodsX_ , 12:102520, 2024.  ISSN 2215-0161.  doi:https://doi.org/10.1016/j.mex.2023.102520.  URL https://www.sciencedirect.com/science/article/pii/S2215016123005162.
  * NH [2022] Riyaz Khan NH.  Study of fluctuations in surface area of lake haramaya using ndwi and mndwi methods.  _JGISE: Journal of Geospatial Information Science and Engineering_ , 5(1):36–41, 2022.