Lenapy for spherical harmonics and gravity data

Lenapy includes functions to read and process spherical harmonics and gravity data.

These function are accessible to be applied to a dataset using xarray accessor. With a xr.Dataset instance ds, one can call the function to be applied onto ds with ds.lnharmo.func().

The corresponding functions are located in the package lenapy.utils (for lnharmo it corresponds to function mainly from lenapy.utils.harmo and lenapy.utils.gravity).

Imports

Lenapy is an overlayer of xarray. For using lenapy fonctionnalities you need to import both xarray and lenapy.

[32]:

import xarray as xr
import lenapy

import os
import numpy as np

# for plotting:
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

Apply L2 -> L2B coefficient corrections

For the following demonstration, we use COST-G products (GRACE RL01 + GRACE-FO operational RL02) between 2002 and December 2022.

[50]:

# Dataset with COST-G months of GRACE and GRACE-FO
ds_path = '../../data/COSTG_n12_2002_2022.nc'

ds = xr.open_dataset(ds_path)

Apply simple spatial filtering

GRACE products need to be filtered to remove north-south striping. It can be done with a gaussian filter in lenapy.

[51]:

from lenapy.utils.gravity import gauss_weights

# radius in meters
ds_gauss_weights = gauss_weights(radius=300000, lmax=12)

ds_filtered = ds.lnharmo * ds_gauss_weights

Read TN14

Both TN14 and TN13 dataset can be download here : https://podaac.jpl.nasa.gov/gravity/grace-documentation#TechnicalNotes

\(C_{2,0}\) coefficient has a large uncertainty for GRACE products compared to other determination with SLR satellites. It is the same for \(C_{3,0}\) coefficient after August 2016 and the accelerometer failure of one GRACE and GRACE-FO satellite.

To improve the quality of GRACE products, supplementary dataset, TN14, that can be used to corrected these coefficients is delivered by the processing centers.

The function read the time-series with or without the mean value using the argument rmmean= (False by default).

[52]:

from lenapy.readers.gravi_reader import read_tn14, read_tn13

ds_C20_C30 = read_tn14('../../data/TN-14_C30_C20_GSFC_SLR.txt', rmmean=False)

To be really precise, you need to change the Earth radius associated to the coefficients used for the correction so that it correspond to the constant associated the the GRACE product used (COST-G here) and then you can replace the coefficients in the product.

The function change_reference() read the ‘radius’ and ‘earth_gravity_constant’ initial and final variables from the dataset.attrs or they can be given with optionnal attributes of the function.

Note that the modification for coefficient :math:`C_{2,0}` and :math:`C_{3,0}` is negligeable

[53]:

from lenapy.utils.gravity import change_reference

ds_C20_C30 = change_reference(ds_C20_C30, new_radius=ds.attrs['radius'], new_earth_gravity_constant=ds.attrs['earth_gravity_constant'])

ds.clm.sel(l=2, m=0).plot(label='COST-G $C_{2,0}$')
# Correct C20
ds.clm.loc[{'l':2, 'm':0}] = ds_C20_C30.clm.sel(l=2, m=0, time=ds.time)
ds.clm.sel(l=2, m=0).plot(label='TN14 $C_{2,0}$')
plt.legend()

plt.figure()
ds.clm.sel(l=3, m=0).plot(label='COST-G $C_{3,0}$')
# Correct C30 after July 2016 for GRACE and GRACE-FO (start on August 2016)
ds.clm.loc[{'l':3, 'm':0, 'time':slice('2016-08-01', None)}] = ds_C20_C30.clm.sel(l=3, m=0, time=slice('2016-08-01', str(ds.time.values[-1])[:10]))
ds.clm.sel(l=3, m=0).plot(label='TN14 $C_{3,0}$')
plt.xlim(np.datetime64('2016-01-01'), np.datetime64('2023-01-01'))
plt.legend()

[53]:

<matplotlib.legend.Legend at 0x71d50c0ffe10>

../_images/tutorials_Tutorial_lnharmo_14_1.png

../_images/tutorials_Tutorial_lnharmo_14_2.png

Read TN13

GRACE products does not contain degree 1 coefficients and are then referenced with an origin corresponding to the center of mass of the Earth. To add degree 1 coefficient, processing centers are delevering TN13, degree 1 time-series product. Each SDS centers has is own TN13 file that is generated with center products using Sun et al. 2016 and Swenson et al. 2008.

To read and apply the correction

ds_deg1 = lenapy.readers.gravi_reader.read_tn13(‘TN-13_GEOC_CSR_RL06.2.txt’)

ds.clm.loc[{‘l’:1, ‘m’:0}] = ds_deg1.clm.sel(l=1, m=0, time=ds.time)

ds.clm.loc[{‘l’:1, ‘m’:1}] = ds_deg1.clm.sel(l=1, m=1, time=ds.time)

ds.slm.loc[{‘l’:1, ‘m’:1}] = ds_deg1.slm.sel(l=1, m=1, time=ds.time)

Manipulating Spherical harmonics datasets

Operation of spherical harmonics dataset

You may have seen that I can apply operation like power (2) over ds.lnharmo object. In fact several methods have been overwriten to be applied only over ds.clm** and ds.slm variables in the dataset.

By this way you can add, subtract, multiply dataset with spherical harmonics variables by applying the operation onto ds.lnharmo object.

Note that you can not apply the divide operation (however, you can multiply by \(\frac{1}{x}\)).

[59]:

ds_sum = ds.isel(time=slice(0)).lnharmo + ds.isel(time=slice(1))

ds_mean = (ds.isel(time=slice(0)).lnharmo + ds.isel(time=slice(1))).lnharmo * 0.5

Plot based on spherical harmonics

Time-serie plot

You can just the time-series corresponding to one coefficient with ds.sel(l=l, m=m).plot()

If you give argument to the .plot() function they are interpreted as plot argument like ‘label’ or ‘color’.

After the plot, you can modify the matplotlib figure as usual.

[60]:

# time series of one coefficient
ds.isel(l=2, m=0).clm.plot()
plt.title("$C_{2,0}$")

[60]:

Text(0.5, 1.0, '$C_{2,0}$')

../_images/tutorials_Tutorial_lnharmo_30_1.png

Triangular plot

With a spherical harmonics dataset, we might try to look at it in terms of SH coefficients amplitude with the function ds.lnharmo.plot_hs() that call the function located in lenapy.plots.plotting.

You first need to create a dataset without extra dimension (reduce the time with a .sel(time=…) or a .mean(‘time’) or other). Then you can plot the dataset.

[61]:

import matplotlib.ticker as ticker

# reduce the dataset to one date
ds_square = ds.isel(time=0).lnharmo ** 2

# triangular plot
ax = ds_square.lnharmo.plot_hs(norm='log', lmax=12)

../_images/tutorials_Tutorial_lnharmo_32_0.png

HS power plot

You can plot the degree power for a certain date with the function plot_power_hs() located in lenapy.plots.plotting.

You first need to create a dataset without extra dimension (reduce the time with a .sel(time=…) or other). Then you can plot the dataset.

[62]:

from lenapy.plots.plotting import plot_power

# plot the degree power
plot_power(ds.isel(time=0), lmin=3, label='COST-G first date')
plot_power(ds.isel(time=1).lnharmo*2, lmin=3, label=r'COST-G second date $\times$ 2')
plt.legend()
plt.title("Degree power spectrum")

[62]:

Text(0.5, 1.0, 'Degree power spectrum')

../_images/tutorials_Tutorial_lnharmo_34_1.png

Writing .gfc file

You can write a .gfc file with ds.lnharmo.to_gfc(). The function dataset_to_gfc located in lenapy.writers.gravi_writer allows to write .gfc ascii file from a spherical harmonics dataset.

The writer will inspect ds.attrs to create the header of the gfc file. If values are not contained in ds.attrs or in the function parameters, default values are used.

First, you need to create a dataset without extra dimensions by reducing the ‘time’ dimension using methods like .sel(time=…) or .mean(‘time’).

[15]:

example_path = 'tmp/test.gfc'

# Write the spherical harmonics data of the first time slice to a .gfc file
ds.isel(time=0).lnharmo.to_gfc(example_path)

Lenapy for spherical harmonics and gravity data

Imports

Reading

ascii file with .gfc format reading

gfct temporal files

Multifiles

ascii file with Level 2 GRACE format

Apply L2 -> L2B coefficient corrections

Apply simple spatial filtering

Read TN14

Read TN13

Projecting spherical harmonics into spatial grid

Details of the projection

EWH (default unit)

Other units

Function for the conversion

Using grid variable

Access residuals and geographic plot

Going back to Spherical Harmonics