Handling processed AQG data¶

Once the raw data is processed and converted into an AQGDataset-class, all further visualization and offered methods are carried out on this variable.

Read processed data from an .nc-file¶

read_aqg_dataset() uses xarray.open_dataset() to create an xarray.Dataset instance, and then wraps it as an AQGDataset.

In [1]:

from gravitools.aqg import read_aqg_dataset
dataset = read_aqg_dataset("../data/test.nc")
dataset

from gravitools.aqg import read_aqg_dataset dataset = read_aqg_dataset("../data/test.nc") dataset

Out[1]:

<AQGDataset: 20240620_163341>

Save to file¶

AQGDataset.to_nc() exports the dataset as NETCDF4, a binary data format. The resulting file contains the full signal of all data columns and all relevant metadata. It is also transparently compressed.

In [2]:

dataset.to_nc("../data/test.nc")

dataset.to_nc("../data/test.nc")

Accessing dataset properties and metadata¶

AQGDataset is a wrapper to an xarray.Dataset, which provides advantages for reading the data from file since individual columns can be read on-demand. The AQGDataset.ds attribute can also be accessed directly, though this is usually not needed. All metadata is contained in the ds dataset and can be accessed via AQGDataset.ds.attrs, which is also accessible with the AQGDataset.metadata property.

In [3]:

dataset.ds.attrs["software"]

dataset.ds.attrs["software"]

Out[3]:

'1.7.4'

In [4]:

dataset.metadata["software"]

dataset.metadata["software"]

Out[4]:

'1.7.4'

Viewing the gravity time series¶

Usually it is not necessary or meaningfull to consider individual drops because the variance is very high and, due to the nature of the measurement principle, successive drops are correlated. For this reason, AQGDataset provides two methods for accessing the gravity signal:

AQGDataset.full_g(): The full gravity signal, showing all drops at 0.54 s intervals. Note that this does not follow a white noise distribution.

AQGDataset.g(): The gravity signal downsampled to 10 second interval means. If there are no interruptions to the measurement, a 10 second interval contains 18 or 19 drops. Importantly, one this time scale, the signal does show a white noise distribution.

In [5]:

# Configure Pandas display float precision
import pandas as pd
pd.options.display.float_format = '{:.1f}'.format

# Configure Pandas display float precision import pandas as pd pd.options.display.float_format = '{:.1f}'.format

In [6]:

dataset.full_g()

dataset.full_g()

Out[6]:

GravityTimeSeries: 20240620_163341
{'vgg': Gradient(value=-3287.0, error=0), 'height': 1.25}
                                         y
utc                                       
2024-06-20 16:33:41.790999889 9808369416.9
2024-06-20 16:33:42.332000017 9808369400.0
2024-06-20 16:33:42.871999979 9808369547.9
2024-06-20 16:33:43.411000013 9808369492.0
2024-06-20 16:33:43.950999975 9808369485.9
...                                    ...
2024-06-20 23:59:57.690000057 9808369449.7
2024-06-20 23:59:58.230000019 9808369411.1
2024-06-20 23:59:58.769999981 9808369495.1
2024-06-20 23:59:59.309999943 9808369625.1
2024-06-20 23:59:59.848999977 9808369506.5

[48690 rows x 1 columns]

In [7]:

dataset.g()

dataset.g()

Out[7]:

GravityTimeSeries: 20240620_163341
{'vgg': Gradient(value=-3287.0, error=0), 'height': 1.25}
                               y  y_err  y_count
utc                                             
2024-06-20 16:33:40          NaN   22.9        6
2024-06-20 16:33:50 9808369288.0   40.9       19
2024-06-20 16:34:00 9808369631.0   25.0       18
2024-06-20 16:34:10 9808369500.5   52.7       19
2024-06-20 16:34:20 9808369449.0   70.4       19
...                          ...    ...      ...
2024-06-20 23:59:20 9808369057.1   16.3       19
2024-06-20 23:59:30 9808369167.6   28.3       18
2024-06-20 23:59:40 9808369313.7   37.6       19
2024-06-20 23:59:50 9808369158.4   46.0       19
2024-06-21 00:00:00          NaN   86.3        9

[2679 rows x 3 columns]

You can obtain the gravity signal at a different transfer height by transferring this result (dataset.g().transfer(h=0)) or by providing a height argument (dataset.g(h=0)).

When further downsampling to 10 minute intervals, for instance, the errors (standard error of the mean) should accurately represent the statistical uncertainty.

In [8]:

dataset.g().resample("10min")

dataset.g().resample("10min")

Out[8]:

GravityTimeSeries: 20240620_163341
{'vgg': Gradient(value=-3287.0, error=0), 'height': 1.25}
                               y  y_err  y_count
utc                                             
2024-06-20 16:30:00 9808369398.0   51.7        7
2024-06-20 16:40:00 9808369346.9   19.5       56
2024-06-20 16:50:00 9808369316.3   22.0       60
2024-06-20 17:00:00 9808369347.0   16.7       60
2024-06-20 17:10:00 9808369356.3   18.3       60
2024-06-20 17:20:00 9808369379.1   17.2       50
2024-06-20 17:30:00 9808369337.6   21.9       60
2024-06-20 17:40:00 9808369329.7   17.1       60
2024-06-20 17:50:00 9808369329.3   14.9       60
2024-06-20 18:00:00 9808369303.5   17.6       60
2024-06-20 18:10:00 9808369309.5   18.3       60
2024-06-20 18:20:00 9808369285.4   19.7       60
2024-06-20 18:30:00 9808369316.5   18.8       60
2024-06-20 18:40:00 9808369376.4   18.7       60
2024-06-20 18:50:00 9808369278.9   16.0       55
2024-06-20 19:00:00 9808369361.0   16.4       60
2024-06-20 19:10:00 9808369298.0   19.3       60
2024-06-20 19:20:00 9808369293.2   17.1       60
2024-06-20 19:30:00 9808369309.8   18.9       60
2024-06-20 19:40:00 9808369339.9   17.4       58
2024-06-20 19:50:00 9808369384.7   17.8       60
2024-06-20 20:00:00 9808369371.9   14.1       60
2024-06-20 20:10:00 9808369353.2   15.6       60
2024-06-20 20:20:00 9808369317.6   17.8       50
2024-06-20 20:30:00 9808369331.5   16.7       60
2024-06-20 20:40:00 9808369405.9   17.7       60
2024-06-20 20:50:00 9808369384.1   18.6       60
2024-06-20 21:00:00 9808369328.3   17.5       58
2024-06-20 21:10:00 9808369311.4   19.0       60
2024-06-20 21:20:00 9808369304.4   18.4       60
2024-06-20 21:30:00 9808369318.1   23.2       60
2024-06-20 21:40:00 9808369326.4   17.2       60
2024-06-20 21:50:00 9808369318.9   16.6       55
2024-06-20 22:00:00 9808369347.4   18.0       60
2024-06-20 22:10:00 9808369314.2   20.9       60
2024-06-20 22:20:00 9808369315.2   18.3       60
2024-06-20 22:30:00 9808369345.2   19.0       60
2024-06-20 22:40:00 9808369377.6   18.2       60
2024-06-20 22:50:00 9808369364.8   15.1       60
2024-06-20 23:00:00 9808369325.0   16.9       60
2024-06-20 23:10:00 9808369324.7   19.1       60
2024-06-20 23:20:00 9808369376.0   18.7       50
2024-06-20 23:30:00 9808369376.8   17.6       60
2024-06-20 23:40:00 9808369285.0   16.1       60
2024-06-20 23:50:00 9808369275.2   17.9       60
2024-06-21 00:00:00 9808369241.0   26.0       30

Naturally, some intervals will not contain as many drops and the statistical error can therefore be much higher. You can exclude these intervals by requiring at least 50% of the max. number of drops within an interval, for instance.

In [9]:

interval = "10min"
g = dataset.g().resample(interval, min_count=0.5)
g

interval = "10min" g = dataset.g().resample(interval, min_count=0.5) g

Out[9]:

GravityTimeSeries: 20240620_163341
{'vgg': Gradient(value=-3287.0, error=0), 'height': 1.25}
                               y  y_err  y_count
utc                                             
2024-06-20 16:30:00          NaN   51.7        7
2024-06-20 16:40:00 9808369346.9   19.5       56
2024-06-20 16:50:00 9808369316.3   22.0       60
2024-06-20 17:00:00 9808369347.0   16.7       60
2024-06-20 17:10:00 9808369356.3   18.3       60
2024-06-20 17:20:00 9808369379.1   17.2       50
2024-06-20 17:30:00 9808369337.6   21.9       60
2024-06-20 17:40:00 9808369329.7   17.1       60
2024-06-20 17:50:00 9808369329.3   14.9       60
2024-06-20 18:00:00 9808369303.5   17.6       60
2024-06-20 18:10:00 9808369309.5   18.3       60
2024-06-20 18:20:00 9808369285.4   19.7       60
2024-06-20 18:30:00 9808369316.5   18.8       60
2024-06-20 18:40:00 9808369376.4   18.7       60
2024-06-20 18:50:00 9808369278.9   16.0       55
2024-06-20 19:00:00 9808369361.0   16.4       60
2024-06-20 19:10:00 9808369298.0   19.3       60
2024-06-20 19:20:00 9808369293.2   17.1       60
2024-06-20 19:30:00 9808369309.8   18.9       60
2024-06-20 19:40:00 9808369339.9   17.4       58
2024-06-20 19:50:00 9808369384.7   17.8       60
2024-06-20 20:00:00 9808369371.9   14.1       60
2024-06-20 20:10:00 9808369353.2   15.6       60
2024-06-20 20:20:00 9808369317.6   17.8       50
2024-06-20 20:30:00 9808369331.5   16.7       60
2024-06-20 20:40:00 9808369405.9   17.7       60
2024-06-20 20:50:00 9808369384.1   18.6       60
2024-06-20 21:00:00 9808369328.3   17.5       58
2024-06-20 21:10:00 9808369311.4   19.0       60
2024-06-20 21:20:00 9808369304.4   18.4       60
2024-06-20 21:30:00 9808369318.1   23.2       60
2024-06-20 21:40:00 9808369326.4   17.2       60
2024-06-20 21:50:00 9808369318.9   16.6       55
2024-06-20 22:00:00 9808369347.4   18.0       60
2024-06-20 22:10:00 9808369314.2   20.9       60
2024-06-20 22:20:00 9808369315.2   18.3       60
2024-06-20 22:30:00 9808369345.2   19.0       60
2024-06-20 22:40:00 9808369377.6   18.2       60
2024-06-20 22:50:00 9808369364.8   15.1       60
2024-06-20 23:00:00 9808369325.0   16.9       60
2024-06-20 23:10:00 9808369324.7   19.1       60
2024-06-20 23:20:00 9808369376.0   18.7       50
2024-06-20 23:30:00 9808369376.8   17.6       60
2024-06-20 23:40:00 9808369285.0   16.1       60
2024-06-20 23:50:00 9808369275.2   17.9       60
2024-06-21 00:00:00 9808369241.0   26.0       30

Obtain mean gravity value¶

Calculate the dataset's mean gravity value (with a statistical error) at the standard transfer height (gravitools.utils.STANDARD_HEIGHT).

In [10]:

g_mean = dataset.mean()
g_mean

g_mean = dataset.mean() g_mean

Out[10]:

AbsGValue(value=np.float64(9808369333.151157), error=np.float64(19.447983338980816), height=1.25, vgg=Gradient(value=-3287.0, error=0))

In [11]:

str(g_mean)

str(g_mean)

Out[11]:

'9,808,369,333.2 ± 19.4 nm/s² (1.25 m)'

The mean gravity at a different height can be obtained by transferring this value or by supplying a custom height to the mean() method.

In [12]:

str(g_mean.transfer(height=0))

str(g_mean.transfer(height=0))

Out[12]:

'9,808,373,441.9 ± 19.4 nm/s² (0 m)'

In [13]:

str(dataset.mean(h=0))

str(dataset.mean(h=0))

Out[13]:

'9,808,373,441.9 ± 19.4 nm/s² (0 m)'

Visualizing the data¶

Set optional Matplotlib styling parameters

In [14]:

import matplotlib.pyplot as plt

plt.rcParams["figure.figsize"] = 8, 4
plt.rcParams["date.autoformatter.hour"] = "%H:%M\n%d %b"
plt.rcParams["axes.ymargin"] = 0.2

import matplotlib.pyplot as plt plt.rcParams["figure.figsize"] = 8, 4 plt.rcParams["date.autoformatter.hour"] = "%H:%M\n%d %b" plt.rcParams["axes.ymargin"] = 0.2

The GravityTimeSeries class provides a plot() method, which uses pyplot.errorbar() to create the graph.

In [15]:

g.plot()

g.plot()

Out[15]:

<ErrorbarContainer object of 3 artists>

You can add more elements to the plot and customize it, just like any other matplotlib plot.

In [16]:

plt.figure()
g.plot("x", color="green", label=f"{interval} means", capsize=3)
plt.axhline(float(g_mean) - g.y0, label="Mean", color="k", linestyle="--")
plt.grid(color="lightgrey")
plt.title(dataset.name)
plt.legend(ncols=2)
plt.show()

plt.figure() g.plot("x", color="green", label=f"{interval} means", capsize=3) plt.axhline(float(g_mean) - g.y0, label="Mean", color="k", linestyle="--") plt.grid(color="lightgrey") plt.title(dataset.name) plt.legend(ncols=2) plt.show()

Plot the gravity signal relative to the mean gravity

In [17]:

g.plot(y0=float(g_mean))
plt.show()

g.plot(y0=float(g_mean)) plt.show()

Create a stacked plot of the gravity signal and another data column¶

In [18]:

fig, axs = plt.subplots(nrows=2, sharex=True)
fig.subplots_adjust(hspace=0)
axs[0].set_title(dataset.name)
g.plot(ax=axs[0], y0=float(g_mean), label=f"{interval} means")
axs[1].plot(dataset.get("dg_pressure"), label="Pressure correction")
for ax in axs:
    ax.set_ylabel(r"$\Delta g$ [nm/s²]")
    ax.legend(loc="upper left")
plt.show()

fig, axs = plt.subplots(nrows=2, sharex=True) fig.subplots_adjust(hspace=0) axs[0].set_title(dataset.name) g.plot(ax=axs[0], y0=float(g_mean), label=f"{interval} means") axs[1].plot(dataset.get("dg_pressure"), label="Pressure correction") for ax in axs: ax.set_ylabel(r"$\Delta g$ [nm/s²]") ax.legend(loc="upper left") plt.show()

Create a standardized measurement report as PDF file¶

The AQGDataset-class provides the functionality to generate a standardized PDF report of the measurements. This report intents to include all relevant information as well as standard plots of many important parameters.

In [19]:

dataset.save_report("../data/AQG_measurementReport.pdf")

dataset.save_report("../data/AQG_measurementReport.pdf")

<Figure size 1360x600 with 0 Axes>

<Figure size 1360x600 with 0 Axes>

Reading multiple processed AQG datasets as one¶

still to come..

Comparison to a reference function¶

still to come..